<%BANNER%>

Modeling Influences of Canal Stage Raises on Groundwater and Soil Water in the C-111 Basin of South Florida

MISSING IMAGE

Material Information

Title:
Modeling Influences of Canal Stage Raises on Groundwater and Soil Water in the C-111 Basin of South Florida
Physical Description:
1 online resource (235 p.)
Language:
english
Creator:
Kisekka, Isaya
Publisher:
University of Florida
Place of Publication:
Gainesville, Fla.
Publication Date:

Thesis/Dissertation Information

Degree:
Doctorate ( Ph.D.)
Degree Grantor:
University of Florida
Degree Disciplines:
Agricultural and Biological Engineering
Committee Chair:
Migliaccio, Kati White
Committee Co-Chair:
Munoz-Carpena, Rafael
Committee Members:
Li, Yuncong
Schaffer, Bruce A
Boyer, Treavor H

Subjects

Subjects / Keywords:
canal -- modeling -- soil -- water
Agricultural and Biological Engineering -- Dissertations, Academic -- UF
Genre:
Agricultural and Biological Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract:
The study purpose was to develop modeling tools to investigate how raises in canal stage proposed under the C-111 spreader canal project could impact groundwater and soil water content in the C-111 basin.  The objectives were to: 1) evaluate the impact of surface water management on soil and limestone water content using Dynamic Factor Analysis (DFA), 2) estimate hydraulic parameters governing canal-aquifer interaction using analytical modeling, 3) simulate water table response to the proposed incremental raises in canal stage and 4) simulate soil and bedrock water content dynamics in response to proposed changes in canal stage. Canal stage significantly (t > 2) drives temporal variation in soil and bedrock water contents, followed by net surface recharge.The effect of water table evaporation was not significant at all sites. DFA Model performance was better at sites with smaller depths to water table ( The developed MODFLOW based model was able to reproduce measured water table elevation, with average Nash-Sutcliffe (NSE) > 0.9 and Root Mean Square Error (RMSE) -3. Soil water content before and after the incremental raises in canal stage were predicted to be significantly different (p<0.001), sitesvthat had surface elevation less than 2.0 m NGVD29 were predicted to experience root zone saturation.
General Note:
In the series University of Florida Digital Collections.
General Note:
Includes vita.
Bibliography:
Includes bibliographical references.
Source of Description:
Description based on online resource; title from PDF title page.
Source of Description:
This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility:
by Isaya Kisekka.
Thesis:
Thesis (Ph.D.)--University of Florida, 2013.
Local:
Adviser: Migliaccio, Kati White.
Local:
Co-adviser: Munoz-Carpena, Rafael.

Record Information

Source Institution:
UFRGP
Rights Management:
Applicable rights reserved.
Classification:
lcc - LD1780 2013
System ID:
UFE0045647:00001

MISSING IMAGE

Material Information

Title:
Modeling Influences of Canal Stage Raises on Groundwater and Soil Water in the C-111 Basin of South Florida
Physical Description:
1 online resource (235 p.)
Language:
english
Creator:
Kisekka, Isaya
Publisher:
University of Florida
Place of Publication:
Gainesville, Fla.
Publication Date:

Thesis/Dissertation Information

Degree:
Doctorate ( Ph.D.)
Degree Grantor:
University of Florida
Degree Disciplines:
Agricultural and Biological Engineering
Committee Chair:
Migliaccio, Kati White
Committee Co-Chair:
Munoz-Carpena, Rafael
Committee Members:
Li, Yuncong
Schaffer, Bruce A
Boyer, Treavor H

Subjects

Subjects / Keywords:
canal -- modeling -- soil -- water
Agricultural and Biological Engineering -- Dissertations, Academic -- UF
Genre:
Agricultural and Biological Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract:
The study purpose was to develop modeling tools to investigate how raises in canal stage proposed under the C-111 spreader canal project could impact groundwater and soil water content in the C-111 basin.  The objectives were to: 1) evaluate the impact of surface water management on soil and limestone water content using Dynamic Factor Analysis (DFA), 2) estimate hydraulic parameters governing canal-aquifer interaction using analytical modeling, 3) simulate water table response to the proposed incremental raises in canal stage and 4) simulate soil and bedrock water content dynamics in response to proposed changes in canal stage. Canal stage significantly (t > 2) drives temporal variation in soil and bedrock water contents, followed by net surface recharge.The effect of water table evaporation was not significant at all sites. DFA Model performance was better at sites with smaller depths to water table ( The developed MODFLOW based model was able to reproduce measured water table elevation, with average Nash-Sutcliffe (NSE) > 0.9 and Root Mean Square Error (RMSE) -3. Soil water content before and after the incremental raises in canal stage were predicted to be significantly different (p<0.001), sitesvthat had surface elevation less than 2.0 m NGVD29 were predicted to experience root zone saturation.
General Note:
In the series University of Florida Digital Collections.
General Note:
Includes vita.
Bibliography:
Includes bibliographical references.
Source of Description:
Description based on online resource; title from PDF title page.
Source of Description:
This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility:
by Isaya Kisekka.
Thesis:
Thesis (Ph.D.)--University of Florida, 2013.
Local:
Adviser: Migliaccio, Kati White.
Local:
Co-adviser: Munoz-Carpena, Rafael.

Record Information

Source Institution:
UFRGP
Rights Management:
Applicable rights reserved.
Classification:
lcc - LD1780 2013
System ID:
UFE0045647:00001


This item has the following downloads:


Full Text

PAGE 1

1 MODELING INFLU E NCES OF CANAL STAGE RAISES ON GROUND WATER AND SOIL WATER IN THE C 111 BASIN OF SOUTH FLORIDA By ISAYA KISEKKA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLM ENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2013

PAGE 2

2 2013 Isaya Kisekka

PAGE 3

3 To my wife Racheal Murungi and my daughters Esther Namanda Kirabo and Abigail Namanda Birungi

PAGE 4

4 ACKNOWLEDGMENTS I am very grateful for the unconditional support I have continuously received from my wife, daughters, parents and siblings. I would also like to thank everyone at the University of Florida department of agricultural and biological en gineering in Gainesville and at the Tropical Research and Education Center (TREC) Homestead for the encouragement and support all of you have given me over time. I would like to give special thanks to Dr Migliaccio for her encouragement and patience in gu iding me and nurturing me over several years Dr Migliaccio, seeing your high level of efficiency in whatever you do has made a positive impact on me. My thanks also go to all my other committee members whose continuous counsel enabled this project to bec ome a success. I would like to thank Dr Muoz Carpena for being such an inspirational role model. To Dr Schaffer thank you for your continuous encouragement and sharing with me your experience with EnviroScan. Dr Li thank you for all the encouragement a nd support especially with field data collection. To Dr Boyer thank you for your continuous encouragement to look at problems more broadly and for sharing new ideas I would also like to thank Ms Tina Dispenza for positively contributing to the success of this project. Finally, I would like to acknowledge the organizations and people that supported this study including : South Florida Water Management District for providing funding and Mr. Vito Strano and Mr. Sam Accursio for allowing us to use their lands

PAGE 5

5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 8 LIST OF FIGURES ................................ ................................ ................................ ........ 10 LIST OF ABBREVIATIONS ................................ ................................ ........................... 14 ABSTRACT ................................ ................................ ................................ ................... 16 CHAPTER 1 BACKGROUND ................................ ................................ ................................ ...... 18 Rationale ................................ ................................ ................................ ................. 18 Literature Review ................................ ................................ ................................ .... 21 Hydrological Monitoring of Water Table and Soil Water ................................ ... 21 Dynamic factor Analysis ................................ ................................ ................... 24 Dynamic Factor Model ................................ ................................ ............... 24 Interpretation of DFA results ................................ ................................ ...... 26 Groundwater Flow Modeling ................................ ................................ ............. 27 Unconfined Aquifer Groundwater Flow Equation ................................ ....... 28 Analytical Modeling of Groundwater Flow ................................ .................. 32 Numerical Modeling of Groundwater Flow ................................ ................. 33 Unsaturated Flow Modeling ................................ ................................ .............. 35 Limitations of Richards Equation ................................ ................................ 38 Soil Water Retention Characteristics and Hydraulic Conductivity in WAVE ................................ ................................ ................................ ..... 39 2 DYNA MIC FACTOR ANALYSIS OF SURFACE WATER MANAGEMENT IMPACTS ON SOIL AND BEDROCK WATER CONTENTS IN SOUTHERN FLORIDA LOWLANDS ................................ ................................ ........................... 46 Introduction ................................ ................................ ................................ ............. 46 Materials and Methods ................................ ................................ ............................ 51 Experimental Site ................................ ................................ ............................. 51 Soil and Bedrock Water Content Monitoring ................................ ..................... 53 Measurement and Estimation of Hydrologic Variables ................................ ..... 56 Dynamic Factor Analysis ................................ ................................ .................. 59 Simple Pr edictive Regression Model for Soil Water Content ............................ 62 Results and Discussion ................................ ................................ ........................... 63 Visual Exploratory Analysis of Experimental Time S eries ................................ 63 Response and Explanatory Variables ................................ .............................. 65 Common Trends ................................ ................................ ............................... 66

PAGE 6

6 Relative Contribution of Explanatory Variables ................................ ................ 68 Predicting Soil and Bedrock Water Contents Using a Simplified Dynamic Factor Analysis Based Model ................................ ................................ ........ 71 Assessing the impact of proposed operational changes in C111 canal stage management on soil and bedrock water content ................................ ........... 73 Conclusions ................................ ................................ ................................ ............ 75 3 SENSITIVITY ANALYSIS AND PARAMETER ESTIMATION FOR AN APPROXIMATE ANALYTICAL MODEL OF CANAL AQUIFER INTERACTION APPLIED IN THE C 111 BASIN ................................ ................................ ............. 92 Introduction ................................ ................................ ................................ ............. 92 Materials and Methods ................................ ................................ ............................ 98 Study Area ................................ ................................ ................................ ........ 98 Hydrological Data Monit oring ................................ ................................ ........... 99 Analytical Model ................................ ................................ ............................. 101 Global Sensitivity Analysis ................................ ................................ .............. 104 Mo rris method ................................ ................................ .......................... 105 ................................ ................................ ........................ 107 GLUE Methodology ................................ ................................ ........................ 107 Resul ts and Discussions ................................ ................................ ....................... 108 Parameter Screening: Morris Method ................................ ............................. 108 ................................ .................... 110 Parameter Estimation for STWT1 Using GLUE ................................ .............. 111 Model Fit ................................ ................................ ................................ ......... 112 Comparison of Estimat ed Parameters to Literature Values ............................ 113 Conclusion ................................ ................................ ................................ ............ 114 4 SIMULATING WATER TABLE RESPONSE TO PROPOSED CHANGES IN SURFACE WATER MANAGEMENT IN THE C 111 AGRICULTURAL BASIN OF SOUTH FLORIDA ................................ ................................ ........................... 125 Introduction ................................ ................................ ................................ ........... 125 Materials and Methods ................................ ................................ .......................... 128 Study Area ................................ ................................ ................................ ...... 128 Hydrologic Time Series ................................ ................................ .................. 129 Groundwater Flow Simulation ................................ ................................ ........ 130 Boundary conditions ................................ ................................ ................ 131 Space and time discretization ................................ ................................ .. 133 Sensitivity analy sis and parameter estimation ................................ ......... 134 Model validation ................................ ................................ ....................... 135 Model Application: Canal stage Operational Adjustment Scenarios ............... 135 Assessing Aquifer Response to Large Storms ................................ ............... 136 Results and Discussion ................................ ................................ ......................... 137 Sensitivity Analysis and Parameter Estimation Results ................................ .. 137 Model Calibration and Validation ................................ ................................ .... 139 Model Application Results ................................ ................................ .............. 141

PAGE 7

7 Results of Aquifer Response to Large Storms ................................ ................ 142 Conclusion ................................ ................................ ................................ ............ 145 5 MODELING SOIL WATER RESPONSE TO SURFACE WATER MANGEMENT 111 BASIN USING WAVE CONSIDERING AND MEASUREMENT UNCERTAINITY ................................ ................................ ....... 162 Introduction ................................ ................................ ................................ ........... 162 Material and Methods ................................ ................................ ........................... 166 Study Area and Experimental Set Up ................................ ............................. 166 Numerical modeling with WAVE ................................ ................................ ..... 168 Parameterization and Sensitivity Analysis ................................ ...................... 170 Estimating uncertainty in measured soil and bed rock water content ............. 173 Model Applications ................................ ................................ ......................... 175 Results and Discussion ................................ ................................ ......................... 175 Sensitivity Analysis and Par ameterization of Vadose Zone Model ................. 175 Soil Water Prediction ................................ ................................ ...................... 178 Evaluation of Soil Water Response to Proposed Incremental Raise s in Canal Stage ................................ ................................ ................................ 180 Conclusion ................................ ................................ ................................ ............ 182 6 CONCLUSIONS AND FUTURE RESEARCH ................................ ....................... 212 Chapter Summaries ................................ ................................ .............................. 212 Chapter 2 ................................ ................................ ................................ ........ 214 Chapter 3 ................................ ................................ ................................ ........ 215 Chap ter 4 ................................ ................................ ................................ ........ 217 Chapter 5 ................................ ................................ ................................ ........ 218 Suggested Future Research ................................ ................................ ................. 22 0 Soil Water Measurement Technologies ................................ .......................... 220 Rapid Water Table Raise Phenomenon and Groundwater Flooding in Biscayne Aquifer ................................ ................................ ......................... 221 Data on Crop P honological and Growth Characteristics ................................ 222 Coupling Vadose Zone Models to Sensitivity and Parameter Estimation Tools ................................ ................................ ................................ ........... 222 Climati c impacts ................................ ................................ ............................. 223 LIST OF REFERENCES ................................ ................................ ............................. 224 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 235

PAGE 8

8 LIST OF TABLES Table page 1 1 Summary of selected analytical solutions to the Boussinesq equation ............... 41 2 1 Dynamic Factor Analysis (DFA) models tested b ased on the following goodness of fit measures: AIC, BIC and C eff ................................ ...................... 78 2 3 Dynamic Factor Analysis results for model 14 with 5 common trends and 2 explanatory variables implemented with non stan dardized time series .............. 80 3 1 Monitoring groundwater wells with descriptors in the C111 basin. ................... 116 3 2 Summary of STWT1 model unce rtain parameters and their probability distributions ................................ ................................ ................................ ...... 116 3 3 GLUE estimated parameters for the canal aquifer interaction model STWT1 .. 116 4 1 Water table elevation monitoring sites with descriptors ................................ .... 147 4 2 Goodness of fit statistics for model calibration for water table elevation predictions using MODFLOW ................................ ................................ ........... 147 4 3 Results for comparison of water table elevation before and after a 6, 9 and 12 cm incremental raise in C 111 canal stage. ................................ ................. 148 5 1 Init ial parameter ranges used in WAVE for simulating soil water content at five sites within the C 111 basin ................................ ................................ ....... 185 5 2 Crop coefficient (Kc) and leaf area index (LAI) for sites 3 and 6 ....................... 185 5 3 Crop coefficient (Kc) and crop leaf area index (LAI) for site 1 .......................... 186 5 4 Crop coefficient (Kc) and crop leaf area index (LAI) for site s 2 and 4 ............... 186 5 5 Improved Morris screening results for WAVE model applied at site 1 for soil layer (L1: 10 and 20 cm) and limestone layer (L2: 30 and 40 cm). ................... 187 5 6 Improved Morris screening results for WAVE model applied at site 2 for soil layer (L1: 10 and 20 cm) and limestone layer (L2: 30 and 40 cm). ................... 188 5 7 Impr oved Morris screening results for WAVE model applied at site 3 for soil layer (L1: 10 and 20 cm) and limestone layer (L2: 30 and 40 cm). ................... 189 5 8 Improved Morris screening results for WAVE mode l applied at site 4 for soil layer (L1: 10 and 20 cm) and limestone layer (L2: 30 and 40 cm). ................... 190

PAGE 9

9 5 9 Improved Morris screening results for WAVE model applied at site 6 for soil layer (L1: 10 and 20 cm) and limestone layer (L2: 30 and 40 cm). ................... 191 5 10 WAVE parameters obtained from calibration at different sites (October 1, 2010 to December 31, 2011) ................................ ................................ ............ 192 5 11 Goodness of fit statistics without consideration of measurement uncertainity .. 193 5 12 Goodness of fit statistics considering measurement uncertainity ..................... 194 5 13 Comparison of volumetric soil water content before and after a 6, 9 and 12 cm incremental raise in canal C 111stage. ................................ ....................... 195

PAGE 10

10 LIST OF FIGURES Figure page 1 1 Map of C 111 basin showing farmlands east of canal C 111, Taylor Slough, Florida Bay and water table and soil water monitoring sites. .............................. 43 1 2 Picture of an EnviroScan installed in the field at C 111 project site. ................... 44 1 3 Showing conceptual model of Biscayne aquifer interacting with canal C 111 for simplici ty this model assumes that C 111 fully penetrates the aquifer. ......... 45 2 1 Map of C111 basin showing farmlands east of canal C 111, Taylor Slough, Florida Bay and three monitoring sites ................................ .............................. 81 2 2 Map of C111 basin showing farmlands east of canal C 111, Taylor Slough, Florida Bay and three monitoring along transect perpendicular to C111 canal. 82 2 3 Temporal variation in scaled frequency at three sites (i.e., T500, T1000 and T200 0 ................................ ................................ ................................ ................ 83 2 4 Temporal variation in hydrologic factors evaluated for their influence on s oil and bedrock water contents at the study site during the period August 2010 to June 2012. ................................ ................................ ................................ ...... 84 2 5 Common trends with 95% confidence interval describing unexplained temporal variation in scal ed frequency. ................................ .............................. 85 2 6 Fitted Dynamic Factor Model and observed temporal variation in scaled frequency in gravely loam soils and limestone bedrock at a site located 500 m along a transect from C111. ................................ ................................ ........... 86 2 7 Fitted Dynamic Factor Model and observed temporal variation in scaled frequency in gravely loam soils and limestone bedrock at a site located 1000 m along a transect from C111. ................................ ................................ ........... 87 2 8 Fitted Dynamic Factor Model and observed temporal variation in scaled frequency in gravely loam soils and limestone bedrock at a site located 2000 m. ................................ ................................ ................................ ....................... 88 2 9 Performance of a simple model for predicting scaled frequency as a function of canal stage and net recharge at specific elevations in parentheses NGVD29 at a site located 500 m. ................................ ................................ ....... 89 2 10 Performance of a simple model for predicting scaled frequency as a function of canal stage and net recharge at specific elevations in parentheses NGVD29 at a site located 2000 m. ................................ ................................ ..... 90

PAGE 11

11 2 11 Boxplots showing soil and rock water content as measured using scaled frequency at site T500 before and after 6, 9 and 12 cm increase canal at structure S18C along C111. ................................ ................................ ............... 91 3 1 Map of the study area showing University of Florida and South Fl orida Water Management District experimental sites and the SFWMD canal network. ....... 117 3 2 Morris screening results for ST WT1 model applied in the Biscayne aquifer .... 118 3 3 Improved Morris sampling screening results for STWT1 model applied in the Biscayne aquifer ................................ ................................ .............................. 119 3 4 Sobol indices for the canal aquifer interaction model STWT1 at 6 groundwater observation wells. ................................ ................................ ........ 120 3 5 Posterior probability density functions (PDF) and cumulative dens it y functions (CDF) at C 111AW ................................ ................................ ........... 121 3 6 Posterior probability density functions (PDF) and cumulative density functions (CDF) at C 111AE ................................ ................................ ............ 122 3 7 Dotty plots from GLUE analysis showing change in the likelihood measure NSE over the range of model parameters ................................ ....................... 123 3 8 STWT1 predicted and measured water table elevation time series along a transect from canal C 111. ................................ ................................ ............... 124 4 1 Study area showing groundwater monitoring sites, agricultural lands adjacent to Everglades National Park, and canal network. ................................ ............. 149 4 2 Composite scaled sensitivities for the parameters selected for estimation in the model. ................................ ................................ ................................ ......... 150 4 3 Validation goodness of fit indicators fr om FITEVAL for MODFLOW simulations at well 2. ................................ ................................ ........................ 151 4 4 Validation goodness of fit indicators from FITEVAL for MODFLOW simulations at well 3. ................................ ................................ ........................ 151 4 5 Validation goodness of fit indicators from FITEVAL for MODFLOW simulations at well 4. ................................ ................................ ........................ 152 4 6 Validation goodness of fit indicators from FITEVAL for MODFLOW simulations at well 5. ................................ ................................ ........................ 152 4 7 Validation goodness of fit indicators from FITEVAL for MODFLOW simulations at well 6. ................................ ................................ ........................ 153

PAGE 12

12 4 8 Temporal variation in water table elevation in reference to ground surface under current canal stage operation criteria at spillway S18C for observation well 2. ................................ ................................ ................................ ............... 154 4 9 Temporal variation in water table elevati on in reference to ground surface elevation under current canal stage operation criteria at S18C for observation wells well 4, well 5, and well 6. ................................ ................................ ......... 155 4 10 Temporal variation in water table el evation in reference to the root zone under proposed incremental raises in canal stage operation at S18C for observation well 4. ................................ ................................ ............................ 156 4 11 Temporal variation in water table elevation in referen ce to the root zone under proposed incremental raises in canal stage operation at S18C for observation well 5. ................................ ................................ ............................ 157 4 12 Temporal variation in water table elevation in reference to the root zon e under proposed incremental raises in canal stage operation at S18C for observation well 6. ................................ ................................ ............................ 158 4 13 Aquifer response to Tropical Storm Isaac at observation wells south of the spillway at S17 7. ................................ ................................ ............................... 159 4 14 Canal aquifer system response to large storms of various sizes for well s north of the spillway at S177 ................................ ................................ ............ 160 4 15 Cana l aquifer system response to large storms of various sizes for well s south of the spillway at S177 ................................ ................................ ........... 161 5 1 Study area showing soil water monitoring sites, agricultural lands adjacent to Everg lades National Park, and canal network within the C 111 basin of south Miami Dade County, Florida. ................................ ................................ ............ 196 5 2 Depiction of the discretizing of the soil profile and location of the EnviroScans 197 5 3 Sobol indices on the vertical axis and parameters for the WAVE model on the horizontal axis as applied to simulate volumetric soil water content at 10, 20, 30 and 40 cm depths at site 1. ................................ ................................ .... 198 5 4 Sobol indices on the vertical axis and parameters for the WAVE model on the horizontal axis as applied to simulate volumetric soil water content at four monitoring depth 10, 20, 30 an d 40 cm at site 2. ................................ .............. 199 5 5 Sobol indices on t he vertical axis and parameters for the WAVE model on the horizontal axis as applied to simulate volumetric soil water content at four monitoring d epth 20, 30 and 40 cm at site3. ................................ ..................... 200

PAGE 13

13 5 6 Sobol indices on the vertical axis and parameters for the WAVE model on the horizontal axis as applied to simulate volumetric soil water content at four monitoring depth 10, 20, 30 and 40 cm at site 4. ................................ .............. 201 5 7 Sobol indices on t he vertical axis and parameters for the WAVE model on the horizontal axis as applied to simulate volumetric soil wat er content at four monitoring depth 10, 20, 30 and 40 cm at site 6. ................................ .............. 202 5 8 Visual comparison of WAVE simulated and EnviroScan measured volumetric soil water content at site 1. ................................ ................................ ............... 203 5 9 Visual comparison of WAVE simulated and EnviroScan measured volumetric soil water content at site 2. ................................ ................................ ............... 204 5 1 0 Visual comparison of WAVE simulated and EnviroScan measured volumetric soil water content at site 3 w. ................................ ................................ ............ 205 5 1 1 Visual comparison of WAVE simulated and EnviroScan measured volumetric soil water content at site 4. ................................ ................................ ............... 206 5 1 2 Visual comparison of WAVE simulated and EnviroScan measured volumetric soil water content at site 6. ................................ ................................ ............... 207 5 1 3 Example FITEV AL output for WAVE model validation at site 1. ....................... 208 5 1 4 Simulated volumetric soil water content under different C 111 canal stage management scenarios at four different depth at site 3. ................................ ... 209 5 1 5 Simulated volumetric soil water content under different C 111 canal stage management scenarios at four different depth at site 4. ................................ ... 210 5 1 6 Simulated volumetric soil water content under different C 111 canal stage management scenarios at four different depth at site 6. ................................ ... 211

PAGE 14

14 LIST OF ABBREVIATIONS Factor loading Mean Standard deviation a Inverse o f Air Entry Value C eff Nash Sutcliffe coefficient of efficiency CERP Comprehensive Everglades Restoration Plan Cv Coefficient of variation DFA Dynamic Factor Analysis E Wate r Table Evaporation ENP Everglades National Park ETo Potential Evapotranspiration GLUE General Likelihood Uncertainty Estimator K Hydraulic conductivity Kc Crop Coefficient LAI Leaf Area Index lam Pore Connectivity Parameter M Common trends MAE Mean Abso lute Error n Curve Shape Parameter NGVD National Geodetic Vertical Datum NSE Nash Sutcliffe Coefficient of efficiency RMSE Root Mean Square Error Rnet Net Recharge SFWMD South Florida Water Management District

PAGE 15

15 Smax Maximum Root Water Uptake tr Residual S oil Water Content ts Saturated Soil Water Content USACE U.S. Army Corps of Engineers WAVE Water and Agrochemicals in the soil, crop and Vadose Environment WTE Water Table Elevation

PAGE 16

16 Abstract of Dissertation Presented to the Graduate School of the Universit y of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy MODELING INFLU E NCES OF CANAL STAGE RAISES ON GROUNDWATER AND SOIL WATER IN THE C 111 BASIN OF SOUTH FLORIDA By Isaya Kisekka August 2013 Chair : Kati W. Migliaccio Cochair: Rafael Muoz Carpena Major: Agricultural and Biological Engineering The study purpose was to develop modeling tools to investigate how raises in canal stage proposed under the C 111 spreader canal project could impact groun dwater and soil water content in the C 111 basin. The objectives were to: 1) evaluate the impact of surface water management on soil and limestone water content using Dynamic Factor Analysis (DFA), 2) estimate hydraulic parameters governing canal aquifer interaction using analytical model ing 3) simulate water table response to the proposed incremental raises in canal stage and 4) simulate soil and bedrock water content dynamics in response to proposed changes in canal stage. Canal stage significantly (t > 2) drives temporal variation in soil and bedrock water contents, followed by net surface recharge. The effect of water table evaporation was not significant at all sites. DFA Model performance was better at sites with smaller depths to water table (< 1 m) highlighting the effect of micro topography on water content dynamics. The following hydraulic parameters were estimated: specific yield (0.07 to 0.14) horizontal hydraulic conductivity (11,000 to 14,300 m/day) aquifer

PAGE 17

17 thickness (13.4 to 18.3 m) and ca nal leakance (99.8 to 279 m) The estimated values were within the range of values estimated using more complex methods at nearby sites. The developed MODFLOW based model was able to reproduce measured water table elevation, with average Nash Sutcliffe (NS E) > 0.9 and Root Mean Square Error (RMSE) <0.05 m. Model predicted that incremental raises in canal stage resulted in significant differences (p<0.05) in water table elevation. Increases in canal stage of 9 and 12 cm resulted in occasional root zone satur ation of low elevation sites and shortening of growing season. The model also predicted that lowering canal stage prior to large storm s reduced water table intrusion into the root zone. C onsidering measurement uncertainty in the WAVE based models resulted in NSE values less than zero to 0.89 and RMSE values between 0. 22 and 1.61 m m 3 Soil water content before and after the incremental raises in canal stage were predicted to be significantly different (p<0.001) sites that had surface elevation less than 2 .0 m NGV D 29 were predicted to experience root zone saturation.

PAGE 18

18 CHAPTER 1 BACKGROUND Rationale rest oration project named the Comprehensive Everglades Restoration Plan (CERP) was formulated CERP was approved by the United States Congress under the Water Resources Development Act (2000). One of the 68 components that comprise CERP is the C 111 spreader c anal project whose goal is to reduce the negative impacts of canal C 111 (i.e., by reducing groundwater seepage into C 111) on Everglades National Park (ENP) and Taylor Slough (Fig ure 1 1 ) which is a natural drainage feature that conveys fresh water to Flo rida Bay while maintaining existing levels of flood protection in the adjacent agricultural and urban areas (U.S. Army Corps of Engineers [USACE] and South Florida Water Management District [SFWMD], 2011). As part of the C 111 spreader canal project, struc tural modifications and operational adjustments involving incremental raises in canal stage are planned along canal C 111 separating ENP and agricultural production areas to the east of the canal. The intention of raising canal stage is to create a hydrau lic ridge along the eastern boundary of ENP which will cause water to flow west towards ENP and southwest towards Florida Bay. The increase in canal stage will occur by changing surface water management at the gated spillway located at structure named S18C (Fig ure 1 1 12 cm (USACE and SFAWMD, 2011).

PAGE 19

19 It is anticipated that the planned rise in C 111 canal stage will affect water table levels in the adjacent agricultural areas. Ear lier research indicated that there is substantial interaction between the highly permeable Biscayne aquifer and water level in canals (Bolster et al., 2001 ; Genereux and Slater, 1999; Ritter and Muoz Carpena 2006). The hydraulic connection between Biscay ne aquifer and canal C 111 causes the shallow water table system to fluctuate with respect to changes in canal stage. Using the drain to equilibrium assumption, Barquin et al. (2011) also showed that water table elevation in the Biscayne aquifer significan tly influenced soil and limestone water contents in a fruit orchard with soil and bedrock formations that are very similar to those in the farmlands east of C 111 Therefore, raising canal stage in C 111 could result in increased water table elevation and increased soil and bedrock water contents or even saturation of the root zone which could affect the production of winter vegetables predominately grown in farmlands east of C 111 Saturation of the root zone could impact yield potential by impairing root growth due to anoxia, reducing stomatal conductance, and reducing net CO 2 assimilation (Schaffer, 1998). In addition to physiological stress, having the soils saturated for long periods could also impact growing season and market dates. Vegetable productio n in Miami Dade County, a substantial proportion of which is located along the extensive eastern boundary of ENP, is a significant contributor to both the local and state economies. In 2011 vegetable produces in Florida harvested vegetables from 185,000 ac res that were worth over 1.5 billion dollars (DACS, 2012). Green beans, sweet corn, squash, tomatoes and sweet potato are the dominant vegetables grown in the area east of C 111 There is need to quantify the impacts of

PAGE 20

20 hydrological modifications and surfa ce water management changes proposed under the C 111 spreader canal project on agricultural land use at field scale because large regional hydrology models have discretization that might not be suitable for resolving small scale micro topographic differenc es within the landscape. However t o date, no comprehensive data set considering the response of water table and soil water to changes in canal stage has been collected and synthesized to simulate this phenomenon in the farmlands east of canal C 111 The study goal was to use monitoring and modeling to investigate how the proposed raises in canal stage along C 111 could impact agricultural production through raised water table elevation and saturated root zone in farmlands east of the canal The following studies were conducted: 1) the impact of surface water management on soil and bedrock water content in the lowland agricultural areas located within the C 111 basin was evaluated using Dynamic Factor Analysis (DFA) 2) hydraulic parameters governing canal aquifer interaction were estimated using an approximate analytical model of groundwater flow and sensitivity analysis and parameter estimation techniques, 3) water table response to the proposed incremental raises in canal stage as well as other large stor m events was simulated using the groundwater flow model MODFLOW and 4) soil and limestone water content dynamics were mechanistically simulated in response to proposed changes in canal stage using the vadose zone model Water and Agrochemicals in the soil, crop and Vadose Environment (WAVE) In chapters 2 to 5 each study is presented in detail as an independent study.

PAGE 21

21 Literature Review Hydrological Monitoring of Water Table and Soil Water Water table elevation is monitored using groundwater observation wel ls. M ulti parameter pressure transducer s ( e.g., Levelogger, Gold Solinst Canada Ltd., 35 Todd Rd, Georgetown, Ontario, Canada) are commonly used to measure groundwater level at various intervals e.g., 15, 30, or 60 minutes Atmospheric corrections are need ed to remove the error s introduced by atmospheric pressure Groundwater level data could be manually downloaded from the level logger or w h ere possible telemetry can be used to access measured data Soil water content in the field can be measured directly or indirectly using a surrogate. Direct measurement of soil water content involves taking a physical sample of soil and weighing it before and after oven drying. The mass of water lost on drying is a direct measure of soil water content (IAEA, 2008). When weight is normalized by dividing with oven dry mass, the technique is called thermo gravimetric (or simply gravimetric) and the soil water content has units of M M 1 alternatively the mass of water lost can be converted to volume by dividing with the dens ity of water, followed by dividing the volume of water with the volume of the soil sample in this case the technique is called thermo volumetric (or simply volumetric) and the soil water content has units of L 3 L 3 Although accurate (0.00028 L 3 ), direct measurement of soil water content is destructive, labor intensive, slow, and does not allow for repetitive measurement of soil water at the same point (Muoz Carpena, 2004). Alternatively, several methods (e.g., neutron moisture meters, thermal sensors, t ime domain reflectometer, capacitances sensors, conductivity sensors, and tensiometers) have been developed for measuring soil water content under field conditions

PAGE 22

22 Unfortunately, none of the alternative methods measures soil water content directly. They i nstead measure a surrogate (e.g., count of slow neutron around a source of fast neutrons in the case of a neutron probe or frequency of an oscillating circuit in case of capacitance probes) which changes when soil water changes. A calibration relationship between the volumetric soil water content and surrogate value must be established. The main advantage of indirect methods is they permit continuous logging of soil water data at the same point and they are nondestructive. The main disadvantage of the indir ect methods is that if they are not well calibrated they can produce erroneous data (Gabriel et al., 2010). Previous research has shown that neutron moisture meters are probably the most accurate indirect methods for measuring volumetric water content but their wide spread use is limited due to regulations and potential health risks (Evett et al., 2006; IAEA, 2008). An alternative to the neutron probe for volumetric water content measurement are methods based on electromagnetic properties. The most widely u sed of the electromagnetic based methods are those based on soil dielectric constant (e.g., capacitance probes). Capacitance probes consist of an electronic circuit called an oscillator and capacitors made up of two hollow metal cylindrical electrodes arr anged coaxially separated by a few millimeters using insulating plastic (IAEA, 2008). The electronic oscillator produces alternating current whose frequency is influenced by the circuit components. There are several commercial brands of capacitance probes ; an example is shown in Figure 1 2 The capacitor(s) form part of the oscillating circuit and electrodes are placed close to the internal walls of the access tube in order to ensure that the fringing electric field of the capacitor(s) would be influenced b y the soil bulk

PAGE 23

23 dielectric permittivity (i.e., a measure of how much an electric field affects and is affected by the bulk soil medium), which is in turn influenced by the water content of the soil (Dean et al., 1987 ). The frequency of oscillation of the a lternating current decreases as the soil water increases. The change in frequency is used as a surrogate for soil water content measurement. One of the drawbacks of capacitance probes is that the volume sensed is very small compared to that of other measur ement techniques such as neutron moisture probes. Another potential source of error is the air gaps between the capacitance probes and the soil usually created during installation. The problem of air pockets is minimized or completely eliminated by ensurin g that the holes drilled for the access are vertical and by using a soil or cement slurry between the access tube and the soil. There also are slight variations in the oscillation frequencies of each sensor. However, this can be overcome through normalizat ion of the measured frequency using f requency readings when the sensor in the access tube is surrounded by water and air Earlier investigations (Evett et al., 2006 ; Gabriel et al., 2010) have shown that factory calibrations are not accurate for all soil s therefore soil specific calibration is necessary. Soil specific calibration of capacitance probes may be conducted under field or laboratory conditions. Soil specific calibration equations are typically obtained by fitting statistical equations between s caled frequency reading s from the capacitance probe and measured volumetric soil content. The volumetric soil water content maybe measured directly using gravimetric methods or by use of a surrogate such as soil suction readings from tensiometers. Gabriel et al. (2010) in their study on comparison of laboratory and field methods for calibrating capacitance probes in Xerochrept soils in

PAGE 24

24 Madrid, Spain observed that manufacturer calibration equation overestimated volumetric soil water compared to the locally d eveloped calibration equation. However, they noted that despite the overestimation of volumetric soil water content, the manufacturer equation was able to reproduce so il water dynamics. Therefore, if the goal is to measure relative differences in soil wa default is sufficient Dynamic factor Analysis DFA is a dimension reduction multivariate time series analysis technique that is used to estimate underlying common patterns (common trends) in short time series as well as the e ffect of explanatory variables on response variables. The advantage of DFA over other traditional dimensional reduction techniques (e.g., Factor Analysis or Principal Component Analysis) is that DFA accounts for the time component. This allows the underlyi ng hidden effects driving the temporal variation in the observed time series data to be detected (Zuur et al., 2003). DFA does not require observed time series to be long and stationary. Although non stationarity could be handled through de trending, trend s in the times series could hold necessary information required to explain the temporal dynamics in the observed variable (Ritter et al., 2009). In addition, DFA can handle missing values in the time series (i.e., DFA does not require data sets to be regul arly spaced). Missing values are not uncommon especially when time series data are obtained from unattended automatic data logging field instruments (e.g., multi sensor capacitance probes for soil water monitoring). Dynamic Factor Model D FA uses a d ynamic f actor m odel (Equation 1 1 ) to describe a set of N observed time series (Ltkepohl, 1991; Zuur et al., 2003; Ritter and Muoz Carpena, 2006).

PAGE 25

25 ( 1 1 ) The goal in DFA is to keep M as small as possible while still obtaining a good mode l fit. Including relevant explanatory variables helps to reduce some of the unexplained variability in the observed time series. Equation 1 1 is expressed mathematically as : ( 1 2 ) ( 1 3 ) where is a vector containing the set of N time series being modeled (response variables), is a vector containing the common trends (same units as the response variables), are factor loadings or weighting coefficie nts that indicate the importance of each of the common trends to each response variable ( unit less ), is a constant level parameter for shifting time series up or down, is a vector containing explanatory variables and are weighting coefficients for the explanatory variables (regression parameters) which indicate the relative importance of explanatory variables to each response variable (inverse units to convert into resp onse variable units), and and are independent, Gaussian distributed noise with zero mean and unknown diagonal covariance matrix. The elements in the covariance matrix represent information that cannot be explain ed by the common trends or the explanatory variables. The unknown parameters and are estimated using the Expectation Maximization (EM) algorithm that is described in Dempster et al. (1977) and Shumway

PAGE 26

26 and Stoffer (1982). The common trends in ( Equation 1 3 ) are modeled as a random walk (Harvey, 1989) and were predicted using the Kalman filter and EM algorithms. The regression parameters in Equation 1 2 are estimated using the same procedure as used in linear regres sion (Zuur et al., 2003). Interpretation of DFA results The results from the DFA are interpreted in terms of the canonical correlations ( ), factor loading ( ), regression parameters ( ) and a greement between modeled and observed time series of the response variable The goodness of fit between modeled and observed values of the response variable are quantified using statistical indicators such as the Nash Sutcliffe coefficient of efficiency ( N SE (BIC). NSE provides an estimate of how well a model predicts an observed data set, while AIC and BIC are relative measures of the g oodness of fit of a statistical model. A model with the NSE closest to 1 and lowest AIC and BIC is the preferred DFA model. Cross correlations between the response time series and common trends were measured using Values of close to unity impl y that the common trend is highly associated with response time series. Typically canonical correlations are classified as follows: | |>0.75, 0.5 0.75, and 0.3 0.5 as high, moderate, and weak correlation s, respectively. The influence of the explanatory variables on response time series is quantified using the magnitude of the coefficients and their associated standard errors which are used with a t test to assess whether the response variable and explanatory variables are significantly related.

PAGE 27

27 Groundwater Flow Modeling Groundwater flow is an example of fluid flow in porous media and can be described using the physical laws governing fluid mechanics such as conservation of mass and c onservation of momentum (Delleur, 1999). The diffusion or groundwater equation ( Equation 1 4 ) is the fundamental mathematical expression used to model groundwater flow and is derived for a representative elementary volume (REV) i.e., the smallest volume ov er which a measurement can be made that will yield a value representative of the whole medium. By pe rforming a mass balance on fluxes flowing in and out of the REV, followed by converting the mass fluxes into volume fluxes using groundwate r flow equation is derived. Under steady state, the groundwater equation transforms into the Laplace equation ( Equation 1 5 ). Detailed information on the derivation of the diffusion equation can be found in Delleur (1999). Some assumptions are made in the derivation of E quation 1 4 for example in the REV concept (Bear, 1972) it is assumed that the average properties of the aquifer are which is a realistic assumption in groun dwater flow since the flow is normally slow with Reynolds number less than unity i.e., laminar flow. Another assumption i s that water is incompressible i.e., the density of water does not vary with change in pressure. Despite the many assumptions, the diff usion equation has been shown to accurately model groundwater flow and is widely used in hydrology models to describe groundwater flow e.g., it is used in MODFLOW, a widely used groundwater flow model developed by the USGS (Harbaugh et al., 2000). ( 1 4 )

PAGE 28

28 ( 1 5 ) w here h(x, y, z, t) is the spatially distributed hydraulic head [ L ], S s specific storage (i.e., the amount of water a portion of the aquifer would release from storage per unit mass per unit drop in hy draulic head while remaining fully saturated), K x K y and K z [ L T 1 ] are components of the hydraulic conductivity tensor and x, y and z are rectangular Cartesian coordinates. The solution to the diffusion equation could be simplified by assuming a homoge nous isotopic aquifer. Equation 1 4 is normally solved numerically considering three types of boundary conditions: a) Dirichlet boundary condition (i.e., specifie d head at a boundary and is normally used when groundwater is in contact with a surface water body), b) Neumann boundary condition (i.e., specifies water flux at a boundary typical examples include wells and infiltration basins) and c) Cauchy boundary condition (combination of the above two boundary conditions, relates hydraulic head to water flux e.g., a canal in contact with groundwater but the interaction between the two being restricted by a layer of sediment). F or transient flow, an initial condition h (x, y, z, 0 ) describing the hydraulic head within the domain of interest prior to start of si mulation is required and is either obtained from observation well data of from a steady solution of the system Unconfined Aquifer Groundwater Flow Equation For unconfined aquifers, the solution to the 3D diffusion equation is complicated by the presence o f a free surface water table boundary condition whose location is also unknown. To describe transient groundwater flow in unconfined aquifers requires an alternative formulation of the fundamental groundwater flow equation. The schematic in Fig ure 1 3 of B iscayne aquifer interacting with canal C 1 1 1 is used to aid discussion of

PAGE 29

29 derivation of the unconfined groundwater flow equation. The governing equation for groundwater flow in unconfined aquifers was first derived by Boussinesq (1877). By examining the wa ter balance of an unconfined aquifer Representative Elementary Dupuit Forchheimer approximations (i.e., no flow in the vertical direction, water table is only slightly inclined, length of aquifer is long compared to saturated thickness [ Figure 1 3 is not to scale as length of the aquifer is much longer than the saturated thickness] ), Bo ussinesq (1877) was able to formulate a general equation for unconfined groundwater flow ( Equa tion 1 6 ); derivation provided in Delleur (1999). ( 1 6 ) where h(x, y, t) is the hydraulic head at the surface of the water table [ L ], Sy is the specific yield/drainable porosity, R is recharge [L T 1 ], and b is the head [ L ] at the ba se of the aquifer, K x and K y [ L T 1 ] are components of the hydraulic conductivity tensor, and x and y are rectangular Cartesian coordinates. Assuming the hydraulic head at the base of the aquifer is zero and the aquifer is homogeneous and isotropic (i.e. hydraulic conductivity does not vary with location and direction) and interacts with the canal through a semipervious sediment layer in addition ignoring recharge E quation 1 6 reduces to Equation 1 7 : ( 1 7 ) Initial condition ( 1 8)

PAGE 30

30 Boundary conditions: ( 1 9) ( 1 10) ( 1 11) ( 1 12) w here l x and l y are the dimensions of the flow domain h o is the initial condition [L] and a [L] is a retardation coefficient proposed by Hantush (1965). The boundary conditions along y=0 and y=l y ( F ig ure 1 3) are defined using D irichlet boundary conditions while the boundary condition at x=0 is defined as a head dependent boundary flux. The ad vantage of Equation 1 7 is that it has no terms in the Z (i.e., vertical) dimension; however, the equation is still nonlinear in h and calculus cannot be easily applied to generate general exact solutions although approximate analytical solutions are possi ble. Therefore, to generate analytical solutions to Equation 1 7 assumptions such as a homogeneous and isotropic aquifer with small changes in water table elevation compared to the saturated thickness of the aquifer (i.e., transmissivity is approximated a s the product of the average saturated thickness of the aquifer and the hydraulic conductivity) are made. Notice also that the terms in the parentheses in Equation 1 7 can be modified using this relationship ( ) which produces a linear ized form of the Boussinesq equation (Wang and Anderson 1982). ( 1 13 )

PAGE 31

31 This form of the equation is commonly applied to conditions similar to those in south Florida (i.e., shallow unconfined aquifer with relatively large saturated t hickness compared to variations in hydraulic head). Equation 1 1 3 can be further simplified by making a change of variables ( ), which produces : ( 1 1 4) Ignoring flow in the y dimension, E quation 1 1 4 reduces to E quation 1 1 5 for which many analytical solutions have been proposed to study 1D aquifer response to stream stage fluctuations: ( 1 1 5) Initial condition ( 1 16 ) Boundary conditions: ( 1 17 ) ( 1 18 ) The linearized form of the Boussinesq equation is associated with under prediction of groundwater table head. However, if the nonlinear form of the Boussinesq equation is solved the problems of under prediction are minimized. Serrano and Workman (1998) using data from a shallow unconfined aquifer in Ohio demonstrated that there were no differences in prediction between analytical solution to the linearized and nonlinearized aquifers when Dupuit assumptions where varied. Groundw ater flow in

PAGE 32

32 the C 111 basin (i.e., farmlands east of C 111) is predominately horizontal and therefore Dupuit assumptions could be assumed as valid. Analytical Modeling of Groundwater Flow Analytical modeling of groundwater flow involves applying various mathematical techniques (e.g., Fourier series, Laplace transform Adomian decomposition method etc.) to generate closed form or approximate analytical solutions to the groundwater flow boundary value problem (i.e., groundwater flow equation with associated initial and boundary conditions). Analytical modeling of groundwater flow has attracted research attention in the past (e.g., Cooper and Rorabaugh, 1963; Hantush, 1965; Hall and Moench, 1972; Serrano and Workman, 1998; Zlotnik and Huang, 1999; Moench and Barlow, 1998 ; Parlange et al., 2000; Lal, 2001; Srivastava et al., 2006). Analytical models are normally the quickest and simplest way to obtain preliminary answers to questions related to characterization of groundwater flow Mathematical expressions (i.e ., analytical models) describing the behavior of the system can provide insight about system behavior without collecting detailed input data needed for complex numerical models. In some cases analytical solutions are used as a method to estimate aquifer p arameters and also to validate simulations from numerical models ( Lal 200 1 ). However, analytical models should be used within the limits of the assumptions under which they were developed. There are numerous analytical solutions in literature and the cha llenge could become selecting the one that is most suitable for studying a particular problem. Examples of selected analytical solutions are provided in Table 1 1 with their contribution and possible limitations. Review of literature indicates that the ana lytical models by Barlow and Moench (2000) are very useful The advantage of the analytical

PAGE 33

33 solutions proposed by Barlow and Moench (1998) is that they are more generalized than previous solutions, they also consider various stream aquifer hydraulic intera ction configurations (e.g., unconfined, confined, semifinite, finite with or without semipervious sediment layers). In addition, they considered anisotropy in the aquifer hydrogeolocal properties yet most of the early analytical solutions simply assumed is otropic aquifer properties. Their solution also attempts to incorporate the effects of transient recharge. Their solutions are presented in a unified manner which provides insight in the mathematical relations among the various aquifer configurations. Last ly the computer programs STLK1 and STWT1 developed by Barlow and Moench (1998) facilitate application of their analytical solution and the convolution algorithm for simulating aquifer response to arbitrary fluctuations in stream stage and recharge. Althou gh closed form and approximate analytical models cannot fully incorporate all the physical complexities of groundwater flow systems, they are simple, easy to implement, versatile and can provide useful preliminary results for many practical applications. Numerical Modeling of Groundwater Flow Numerical modeling of groundwater flow involves solving the partial differential equations (PDE) describing groundwater flow, together with the boundary and initial conditions using step wise approximation. The numer ical approach is normally preferred for highly irregular domain geometries, for highly heterogeneous and anisotropic hydrogeological systems or when fully coupled simulations of stream aquifer interactions are needed (i.e., suitable for more complex simula tions). Numerical solutions require discretization of the domain of interest in both space and time. To numerically obtain a solution to the boundary value problem, the continuous state variable (e.g., hydraulic head) is replaced by discrete values defined at the center of

PAGE 34

34 each block or nodes. T he most frequently used numerical techniques used to solve the groundwater flow boundary value problem include finite difference (FD) and finite element (FE) techniques with each having various schemes of implementa tion. In both techniques, the continuous partial differential equation describing hydraulic head everywhere in the aquifer is replaced by a finite number of algebraic equations describing head at each block or node. The system of equations is then solved u sing matrix techniques or iteration The solution is improved by reducing the size of the time step in much the same way as drawing a curve using a series of very short straight lines. In general, FD methods are simpler conceptually and easier to implement and program compared to FE methods although the later are preferred for complex domain configurations. Comprehensive discussion of numerical modeling of groundwater flow processes can be found in texts such as McDonald and Harbaugh (1988), Huyakorn and P inder (1983), Anderson and Woessner (1992), Konikow (1996) and Konikov and Reilly (1998). There are numerous generic model codes that could be used for numerical simulation of stream aquifer interactions (e.g., MODFLOW, MIKE SHE, FEFLOW, HydroGeoSphere, M ODHMS, etc.). The following websites among others maintain lists of generic model codes, some of which could be used for numerical simulation of canal aquifer interactions problems similar to those in south Florida : http://igwmc.mines.edu/ software/freeware_list.html and http://water.usgs.gov/software/lists/groundwater The following criteria could be used in selecting a gener ic model code for simulating groundwater flow problems : (1) code should be capable of simulating major processes (e.g., unconfined groundwater flow and stream seepage) identified at the conceptual

PAGE 35

35 stage of model development, (2) accuracy of numerical schem e implemented in numerical code needs to be known, (3) computational efficiency (i.e., the code should execute fast), (4) code availability (i.e., public domain or proprietary) and documentation, (5) availability of preprocessing an d post processing grap hical user interfaces and (6) more recently the ability for the generic code to automatically integrate calibration, sensitivity and uncertainty analysis algorithms. N umerical models are approximations of the complex natural processes and therefore good ju dgment is recommended in making the tradeoffs between accuracy and cost (i.e., in terms of input data and computation) during generic model code selection. Unsaturated Flow Modeling The unsaturated zone is the portion of the subsurface in which pore space s are partly filled with both gases (air and water vapor) and water. The theory of soil water flow through the unsaturated zone dates back to the findings of Buckingham (1907), who developed the relationship between flow rate and suction gradients in the u nsaturated zone. Richards (1931) extended the work on unsaturated soil water dynamics by developing the classical partial differential equation governing unsaturated Buckingham (E quation 1 19 ) with the p rinciple of conservation of mass. Equation 1 19 is the derivation of Richards equation in the vertical direction subject to the following assumptions : homogeneous isotropic isothermal rigid porous media hysteresis is neglected and flow is one dimensiona l : ( 1 19 )

PAGE 36

36 were q the is Darcy flux [LT 1 ], is unsaturated hydraulic conductivity [LT 1 ], t is time [T], h is soil water pressure head [h], and z is vertical co ordinate [L]. Combining Equation 1 19 with the continuity equation i .e., by substituting q in Equation 1 20 with Equation 1 19 we obtain Richards equation for unsaturated flow in porous media ( Equation 1 22). ( 1 20 ) ( 1 21 ) Equation 1 21 is known as the mixed form of Richar equation because it written in both the soil water content [L 3 L 3 ] and the soil water pressure head h [L] terms By introducing the differential soil water capacity [L 1 ] which represents the slope of the soil water retention curve and expressing hydraulic conductivity as a function of s oil water pressure head, Equation 1 21 can be written as Equation 1 22 Equation 1 22 is known as the h based formulation of the Richards equation and is written with only one unknown term h Equation 1 22 is implemented in the vadose zone model WAVE ( Van clooster et al., 1995) The h based formulation is applicable in both saturated and unsaturated conditions, in the earlier case the equation is parabolic while in the latter case it reduces to an elliptic equation because the capacitance term becomes zero. ( 1 22 ) Alternatively, Equation 1 21 can be transformed into a form known as the based formulation of equation ( Equation 1 23 ) B y introducing a term on the right

PAGE 37

37 hand side of Equation 1 21 known as the soil water diffusivity [L 2 T 1 ] produces Equation 1 23 ( 1 23 ) To completely describe a transient unsaturated flow problem requires specifying initial condition and boundary conditions. In the case of one dimensional vertical flow, initial condition comprising soil water head or soil water content values at each node in the soil profile is required. If measured values are not available, values of soil water head when the soil water is in equilibrium with t he water table can be used. The upper boundary condition is determined by infiltration and evaporation flux and can either be described as a Dirichlet or Neuman boundary condition. In the WAVE model which was used to simulate soil water dynamics in this in vestigation, the boundary condition at the top of the soil surface is usually a flux or Neuman boundary condition and the flux calculated as: ( 1 24 ) wh ere [LT 1 ] is the potential flux through the soil surface (positive upwards) under non limiting flow conditions (examples of limiting conditions are high rainfall intensities and prolonged soil evaporation without re wetting) [LT 1 ] is the potential soil evaporation rate as determined by weather conditions Rain is the precipitation rate [LT 1 ], Irr is the irrigation amount [LT 1 ], and Pond and Inc are the ponding depth at the soil surface and canopy storage respectively ( Vanclooster et al., 1995). Similarly the bottom boundary condition can also be represented as a Dirichlet or Neuman boundary

PAGE 38

38 condition. In WAVE the bottom boundary condition can be represented as: 1) groundwater table, 2) know n pressure head as a function of time, 3) free drainage and 4) as lysimeter with free outflow at the bottom. In E q uations 1 21 to 1 23 ; the hydraulic conductivity, differential soil water capacity and diffusivity are all non linear functions of soil water pressure head which makes equation non linear and difficult to produce exact solut ions using calculus. Analytical solutions are only possible for specific boundary conditions (Philip, 1957). Thus, equation is usually solved using numerical techniques. A quation is provided in Celia et al. (1990) which note that numerical solutions based on the standard h based form generally may yield poor results, characterized by large mass balance errors and incorrect estimates of infiltration depth. On the other hand, There are several numerical solvers for equation e.g., WAVE, HYDRUS (1D, 2D and 3D), and SWACROP (Soil WAter and CROP production model) Selection for the appropri ate solver depends on the problem to be solved The criteria outlined under numerical modeling of groundwater flow can be applied in selecting generic computer solvers for unsaturated flow. Limitations of Richards Equation equation the basis fo r most unsaturated flow models holds for only stable flow conditions. Unstable flow conditions are not uncommon e.g., under abrupt increases in hydraulic conductivity with depth such as found in south Florida w h ere the thin scarified gravely loam soil lay er has lower hydraulic conductivity compared to the underlying limestone bedrock ( Muoz Carpena et al., 2008). equation is also

PAGE 39

39 not readily applicable under conditions in which the soil profile has preferential flow through non capillary macro po res. This condition is also present in south Florida the highly gravel soil layer and limestone bedrock are both characterized by having macro pores which may cause preferential flow (Kumar, 2002). Compression of air ahead of the wetting front also intro duces errors when simulating unsaturated flow using equation. Version 3 of WAVE has algorithms to simulate unsaturated flow in soils with dual porosity and also includes parametric hysteresis models characterizing soil water retention Soil Water Retention Characteristics and Hydraulic Conductivity in WAVE One of the most commonly used empirical equation to describe the soil water characteristic curve is the van Genuchten equation (van Genuchten, 1980): ( 1 25 ) where is volumetric soil water content, is residual soil water content, is saturated soil water, is soil suction [L], is related to the inverse of the suction at the air entry value, is a parameter related to the measure of the pore size distribution, and is a shape parameter. In an attempt to develop a closed form solution for hydraulic conductivity, van Genuchten (1980) related the parameters n and m through the following expression m=1 1/n Under non hysteresis conditions, Equation 1 26 is among the many soil water retention models implemented in WAVE and is also the most accurate ( Viaene et al., 1994 ). The hydra ulic conductivity curve of soils can be expressed as a function soil water head, volumetric soil water content or effective saturation. The parameters of the hydraulic conductivity curve are determined using

PAGE 40

40 optimization fitting techniques using hydraulic conductivity data measured at various values of soil water head or volumetric soil water content. Different models are available in WAVE to describe hydraulic conductivity e.g., Garder (1958), Brooks and Coorey, (1964) van Genuchten, (1980) and Mualem (197 6).

PAGE 41

41 Table 1 1. Summary of selected analytical solutions to the Boussinesq equation Reference Assumptions Limitations Barlow et al. and Moench ( 1998 ) Aquifer is homogeneous and of uniform thickness. Aquifer can be anisotropic provided that the principal d irections of the hydraulic conductivity tensor are parallel to the coordinate axes Lower boundary of each aquifer type is horizontal and impermeable. Water is instantaneously exchanged between the vadose zone and aquifer in response to a decline elevation of the water table. Change in saturated thickness of the aquifer due to stream stage fluctuations is small compared with the initial saturated thickness The semipervious stream bank material, if present, has negligible capacity to store water. Initially, the water level in the stream is at the same elevation as the water level everywhere in the aquifer. Assumes linear system (i.e., not suitable for systems characterized by large variations in water table head). Assumes full canal penetration and may not be application for very shallow and wide streams. Ignores capillarity and hysteresis in the soil layer above the water. May not be applicable for highly heterogeneous systems with complex boundary conditions. Zlotnik and Huang (1999) Aquifer is homogeneous, isotropic, and semi infinite Water table is initially horizontal and at equilibrium with stream. Streambed partially penetrates the aquifer Drawdown in the aquifer is small compared with the saturated aquifer thickness. Leakage across the streambed is vertical and occurs at the bottom. Flow in the portion of the aquifer under the stream is confined Assumes linear system (i.e., not suitable for systems characterized by large variations in water table head) Ignores recharge May not be applicable for h ighly heterogeneous systems with complex boundary conditions.

PAGE 42

42 Table 1 1 Continued Reference Assumptions Limitations Serrano and Workman (1998) Assumes flow is predominantly horizontal. Length of aquifer is long compared to saturated thickness Aquifer is homogeneous, isotropic, and infinite in lateral extent. Stream and aquifer are in perfect hydraulic connection Stream fully penetrates the aquifer Ignores the effect of the low permeability stream bank sediment layer which could lead to over prediction of aquifer responses. Assuming full stream penetration may also led to over prediction May not be applicable for highly heterogeneous systems with complex boundary conditions. Assumes constant recharge Hall and Moench (1972) Aquifer is homogeneous, isotrop ic, and can be semi infinite or finite in lateral extent Changes in water table levels are small compared to the saturated thickness The semipervious stream bank material, if present, has negligible capacity to store water. The stream fully penetrates the aquifer Flow is predominantly horizontal Aquifer storativity could be substituted with specific yield in order to approximate unconfined aquifer flow. Initially, the water level in the stream is at the same elevation as the water level everywhere in the aq uifer Assuming full stream penetration may also led to over prediction Assumes linear system (i.e., not suitable for systems characterized by large variations in water table head). Assumes full canal penetration and may not be application for very shallow and wide streams. Ignores recharge May not be applicable for highly heterogeneous systems with complex boundary conditions.

PAGE 43

43 Figure 1 1 Map of C 111 basin showing farmlands east of canal C 111 Taylor Slough, Florida Bay and water table and soil wa ter monitoring sites

PAGE 44

44 Figure 1 2 Picture of an EnviroScan installed in the field at C 111 project site (Photo courtesy of Isaya Kisekka) Oscillating electric circuit Cylindrical electrodes comprising the capacitor Access tube

PAGE 45

45 Figure 1 3 Showing conceptual model of Biscayne aquifer interacting with canal C 111 for simplicity this model assumes that C 111 fully penetrates the aquifer.

PAGE 46

46 CHAPTER 2 DYNAMIC FACTOR ANALYSIS O F SURFACE WATER MANAGEMENT IMPACTS O N SOIL A ND BEDROCK WATER CONTENTS I N SOUTHERN FLORIDA LOWLANDS Introduction In an attempt to correct some of the undesi restoration project named the Comprehensive Everglades Restoration Plan (CERP) is currently under implementation. CERP was approved by the Unit ed States Congress under the Water Resources Development Act (2000). One of the 68 components that comprise CERP is the C111 spreader canal project whose goal is to reduce the impacts of C111 (i.e., reduce groundwater seepage into C111) on Everglades Natio nal Park (ENP) and Taylor Slough which is a natural drainage feature that conveys water to Florida Bay while maintaining existing levels of flood protection in the adjacent agricultural and urban areas (U.S. Army Corps of Engineers [USACP] and South Florid a Water Management District [SFWMD], 2011 ). As part of the C111 spreader canal project, structural modifications and operational adjustments involving incremental raises in canal stage are planned along one of the major canals (i.e., C111) separating ENP a nd agricultural production areas to the east of the canal. The increase in canal stage will occur by changing surface water management at the gated spillway located at structure named S18C ( Figure 2 1 ) in the form of incremental raises in canal stage of up to 12 cm. It is anticipated that the planned rise in C111 canal stage will affect water table levels in the adjacent agricultural areas. Earlier research indicated that there is an Copyright for this chapter of the dissertation belongs to Journal of Hydrology Kisekka et al. 2013a

PAGE 47

47 interaction between the highly permeable Biscayne aquifer and water level in canals (Genereux and Slater, 1999). The hydraulic connection between Biscayne aquifer and canal C111 causes the shallow water table system to fluctuate with respect to changes in canal stage. Using the drain to equilibrium assumption, Barquin et al. (2 011) showed that water table elevation in the Biscayne aquifer significantly influenced soil and limestone bedrock water contents in a fruit orchard with soil and bedrock formations that are very similar to our current study site. Therefore, raising water table elevation could result in increased soil and bedrock water contents or greater saturation of the root zone which could affect the production of winter vegetables predominately grown in the area. Saturation of the root zone could impact yield potentia l by impairing root growth due to anoxia, reducing stomatal conductance, and reducing net CO 2 assimilation (Schaffer, 1998). In addition to physiological stress, having the soils saturated could render movement of machinery difficult and also impact growin g season and market dates. However, it is not known to what extent the proposed structural modifications and operational adjustments along canal C111 would impact water table elevations and thus soil and bedrock water contents in agricultural areas east of the canal. Vegetable production in Miami Dade County, a substantial proportion of which is located along the extensive eastern boundary of ENP, is a significant contributor to both the local and state economies. In 2011 vegetable produces in Florida har vested vegetables from 185,000 acres that were worth over 1.5 billion dollars (DACS, 2012). Green beans, sweet corn, squash, tomatoes and sweet potato are the dominant vegetables grown in the area east of C 111 There is need to quantify the impacts of hy drological modifications and surface water management on agricultural land use at

PAGE 48

48 field scale because large regional hydrology models have discretization that might not be suitable for resolving small scale micro topographic differences within the landscap e. Long term monitoring and exploratory analysis of soil and bedrock water contents could characterize the effect of various drivers on the temporal variability of water contents. The soils in the agricultural areas east of C111 were created from scarifica tion of the underlying limestone bedrock hence they are very shallow and have high gravel content. Three main stresses that influence soil water content that could be included in exploratory analysis are 1) canal stage, which affects water table elevation; 2) rainfall ; and 3) evapotranspir ation. While these stresses may be assessed using physically based models of vadose zone flow and transport, implementation of unsaturated flow not an easy task since they contain numerous parameters and processes that have to be quantified (Ritter et al., 2009). In very gravelly and shallow soils such as those in south Miami Dade County, quantifying parameters such as hydraulic conductivity for further complicated by having porous gravely soils that are not homogeneous. Previous applications of WAVE, for example, in gravely soils of south Florida have indicated that a detailed description of soil hydraulic properties (e.g., using dual porosity) could result in improved robustness of vadose zone models (Duwig et al., 2003; Muoz Carpena et al., 2008). Therefore the success of applying physically based models to simulate soil and bedrock water dynamics depends largely on proper conceptualization of location specific processes and proper measurement or estimation of parameters. In this context, complementary exploratory tools such as Dynamic Factor Analysis (DFA) which are not processes based are desired as simpler preliminary

PAGE 49

49 exploratory tools that could also be used for preliminary predictions of the impact of surface water management decisions on land use. A comprehensive description of DFA and modeling can be found in Zuur et al. (2003). For purposes of aiding discussion, we only provide a brie f description of this technique. DFA is a dimension reduction multivariate time series analysis technique that is used to estimate underlying common patterns (common trends) in short time series as well as the effect of explanatory variables on response va riables. The advantage of DFA over other traditional dimensional reduction techniques (e.g., Factor Analysis or Principal Component Analysis) is that DFA accounts for the time component. This allows the underlying hidden effects driving the temporal variat ion in the observed time series data to be detected (Zuur et al., 2003). DFA does not require observed time series to be long and stationary. Although non stationarity could be handled through de trending, trends in the times series could hold necessary in formation required to explain the temporal dynamics in the observed variable (Ritter et al., 2009). In addition, DFA can handle missing values in the observed time series (i.e., DFA does not require data sets to be regularly spaced). Missing values in obse rved time series data sets are not uncommon especially when time series data are obtained from unattended automatic data logging field instruments (e.g., multi sensor capacitance probes for soil water monitoring). DFA applications are documented in litera ture from several disciplines (e.g., Geweke, 1977; Mrkus et al., 1999; Zou and Yu, 1999; Zuur et al., 2003; Zuur and Pierce, 2004; Muoz Carpena et al., 2005; Ritter and Muoz Carpena, 2006; Zuur et al., 2007; Ritter et al., 2009; Kaplan et al., 2010a; Ka plan and Muoz Carpena,

PAGE 50

50 2011).Thus, we only provide a brief review of the most relevant examples. Ritter and Muoz Carpena (2006) applied DFA and modeling to study interactions between surface water and groundwater levels within the Frog Pond agricultural area located west of canal C111 in south Florida ( Figure 2 1). Their results indicated that the two canals surrounding the Frog Pond area had the greatest influence on temporal changes in water table elevation. Their study did not address the issue of the impact of surface management decisions on soil water content. Soil water is a major concern for vegetable growers in south Florida due to the impact saturated or near saturated soil conditions have on planting dates and yield losses ( Figure 2 1). Others have applied DFA and modeling to study soil water dynamics. Ritter et al. (2009) applied DFA to analyze temporal changes in soil water status of a humid, subtropical, evergreen forest in Canary Islands, Spain. Kaplan and Muoz Carpena (2011) applied DFA t o study the complementary effects of surface and groundwater on soil water dynamics in a coastal flood plain. Thus, DFA was successfully used to identify unexplained variability in observed hydrologic time series and to assess the effect of selected explan atory variables on response variables (observed time series of interest). The difference between our study and prior studies is that we applied DFA to investigate the effect of surface water management in canals on soil water dynamics in an agricultural a rea with very shallow very gravely loam soils, and unlike in the previous studies we also considered not only the effects of potential evapo transpi ration (ET o ) but also the effect of water table evaporation given the shallow water table. We then attempted to develop a simple model, using information from the DFA, to predict soil water content from easily measured variables such as canal stage and recharge (i.e.,

PAGE 51

51 difference between rainfall and evapotranspiration). Canal stage was selected instead of water t able elevation since water table elevation data in our study area are less complete due to the limited period of record and the limited number of continuously monitored groundwater wells. Canal stage has been monitored for a longer period of record and has no foreseeable end of data collection, thus it is a more reliable measurement for long term use. We assumed that at any given time, water table elevation is approximately equal to canal stage. We concede that at certain times this assumption might not hol d e.g., immediately after or during storm events; however, due to the high permeability of the aquifer and the daily time step used, the assumption holds for the majority of the time. The goal of this study was to use DFA and modeling to investigate how th e proposed raises in canal stage along C111 could impact soil and bedrock water contents in low lying farmlands located between canals C111 and C111E. The specific objectives were to: (1) apply DFA to identify the most important factors affecting temporal variation in soil and bedrock water contents, (2) develop a simplified DFA based regression model for predicting soil and bedrock water contents as a function of canal stage, and (3) use the developed simple regression model to predict the impact of propos ed incremental raises in canal stage on soil and bedrock water contents at various elevations and distances from the canal. Materials and Methods Experimental Site The study was conducted in southern Miami Dade County, Homestead, Florida, United States in a small agricultural area approximately 17 km 2 (Fig ure 2 1 ). The area is located east of ENP between SFWMD canals C111 and C111E which are planned to

PAGE 52

52 experience increases in canal stage under the C111 spreader canal project. Canal stage upstream in the two canals is controlled by a remotely operated spillway at S177 and a culvert at S178, respectively ( Figure 2 1). C111 is the larger of the two canals and the two join to become a single canal at the southern end of the study area which is managed using a ga ted spillway at S18C. It is proposed that stage will be increased by modifying operation of S18C and thus affect canal stage in the reach of C111 between S177 and S18C. The hydrogeological system at the study site consists of the Biscayne aquifer which is a highly permeable shallow unconfined aquifer with hydraulic conductivities reported to exceed 10,000 m/day, which explains the high connectivity between the canals and the aquifer (Chin, 1991). The shallow nature of the water table implies that evaporatio n from the groundwater could impact soil water content. The topography at this site is essentially flat with elevation ranging approximately between 1.2 to 2.0 m above sea level NGVD 29. The climate is subtropical with a dry season (November to May), which is the growing season for vegetables, and a wet season (June to October). Approximately two thirds of all the rain (average annual rainfall ranges between 1100 to 1524 mm) is received during the wet season months. The soil at the study site is very shall ow (10 to 20 cm) with underlying limestone bedrock. According to Nobel et al. (1996), the soils east of C111 vary and could be classified as either Krome and Chekika very gravely loam (loamy skeletal, carbonatic, hyperthermic, Lithic Undorthents), or Bisca yne Marl (loamy, carbonatic, hyperthermic) based on their physical characteristics. We performed particle size analysis using a standard 2 mm sieve and determined that the soils contain on average of 45% fine fractions and 55% gravel. Color analysis using the Munsell soil color charts (Munsell soil

PAGE 53

53 charts, 2000) and the color guide in Noble et al. (1996) identified the study site soils to be broadly characterized as Chekika soil series. Three monitoring sites were used located at 500, 1000 and 2000 m along a transect perpendicular to canal C111, the three sites also had varying topographies and represented areas expected to experience the greatest impact from the proposed raises in canal stage. Sites were selected to capture differences in soil texture with in our study area; this was done with a soil survey map and site visits. Sites were also selected to ensure they were in privately owned agricultural low lying lands that were expected to be impacted by the rises in water table elevation. For each site: i) GPS coordinates and elevation data were collected, ii) groundwater wells were constructed and each was equipped with level loggers (Levelogger, Gold Solinst Canada Ltd., 35 Todd Rd, Georgetown, Ontario, Canada) to record water table elevation every 15 min utes, iii) multi sensor capacitance probes (MSCP) (EnviroScan probes, Sentek Technologies, Ltd., Stepney, Australia) were installed at each site to monitor soil and bedrock water contents. Monitoring site locations are shown in Figure 2 1; elevations are s hown in Figure 2 2 Differences in the length of times series at the three sites was due to differences in the dates of installation of the EnviroScan probes (i.e., probes could only be installed when water was at least 50 cm below the ground surface) and relocation of the probes due to initial poor installation. Site T500 was installed on August 25, 2010, while sites T1000 and T500 were installed on January 21, 2011. Soil and B edrock W ater C ontent M onitoring Two EnviroScan probes were installed at each s ite for a total of six. Each access tube with a diameter of 50.5 mm housed four sensors positioned at various elevations as shown in Figure 2 2. The elevations correspond to 10, 20, 30 and 40 cm from the

PAGE 54

54 ground surface at each site. The top 20 cm typically represent the scarified soil layer which is used for crop production and the lower 20 cm represent the underlying limestone bedrock in which plant roots cannot penetrate. To minimize the problem of air pockets, we used fast setting cement slurry between t he access tube and the soil. The purpose of installing two EnviroScan probes at the same location was to ensure that at least one probe was functioning at any given time. Due to the shallowness of the limestone bedrock at all the study sites, a motorized d rill was required to bore a hole that held the access tube in a vertical position. Water content data were logged every 15 minutes and were downloaded weekly and averaged daily. EnviroScans are an example of capacitance based sensors which measure frequen cy of an oscillating electrical circuit. The oscillator is coupled electrically to capacitive elements that are made of two metal cylindrical electrodes. The electrode system is arranged so the soil becomes part of the dielectric medium affected by the fri nging electromagnetic field. Volumetric soil water content affects the electrical permittivity of the soil which in turn affects the capacitance causing the oscillation frequency to shift (IAIA, 2008) since the soil dielectric constant is a combination of mineral particles (2 4), water (80), and air (1). According to Dean et al. (1987) the oscillatory frequency from the capacitance soil water sensor could be expressed as : ( 2 1) w here C b is the total base capacitance and C c is the tot al collector capacitance and these represent capacitances of internal circuit elements to which the electrodes are connected, L is the inductance of the coil in the circuit, and C is the capacitance of the soil access tube system. Therefore capacitance of the soil access tube system, C can

PAGE 55

55 g representing the geometry of the sensor as shown : ( 2 2) Differences in oscillatory frequency among sensors at the same soil and bedrock water contents were eliminated by normalizing the oscillatory frequency values using values of frequency when the sensor was surrounded by water and air. The normalized oscillatory frequency is known as the scaled frequency (SF) and is est imated as in eq. 2 3. The manufacture default calibration equation ( E q uation 2 4) can be used to convert scaled frequency to volumetric soil water content ( ( 2 3) ( 2 4) where F is the oscillatory frequency value measured by the EnviroScan sensor, F a is frequency value when the EnviroScan probe is surrounded by air, and F w is the frequency value when the Envir oScan probe is surrounded by water. To avoid location specific calibration for each sensor, we use d effect of various factors on soil and bedrock water contents and thus did not use E q uation 2 4. This approach wa s successfully applied by Ritter et al. (2009) when studying the effect of various factors on hydrologic fluxes in a forest top soil using refractive index from time domain reflectometry (TDR) as a surrogate for volumetric soil water content. Gabriel et al equation overestimated volumetric soil water compared to the locally developed calibration equation. However, they noted that despite the overestimation of volumetric soil water content, the manufacture

PAGE 56

56 water dynamics. Therefore, if the goal is to measure relative changes in water content Measurement and E stimation of H ydrologic V ariables Hydro logic variables including canal stage, water table elevation NGVD29 m, rainfall ( P ), potential evapotranspiration ( ET o ) and groundwater evaporation ( E ) were measured or estimated to assess their influence on soil and bedrock water content time series. Cana l stage data were measured at the S177 spillway for headwater (S177H) and tail water (S177T) every 15 minutes but daily averages were used. Canal stage data were measured by the SFWMD and are publically available from the online environmental database (DBh ydro; http://www.sfwmd.gov/dbhydroplsql/show_dbkey_ info.main_menu ). During the first phase of the C111 spreader canal project, the main operational adjustments will involve incre mentally raising canal stage at S18C ( Figure 2 1) which will result in increased stage in the reach of C111 between the spillways at S177 and S18C Water table elevation data were collected from three observation wells constructed at the three monitoring s ites. Water table elevation was measured by the University of Florida (UF) every 15 minutes and averaged daily using a multi parameter pressure transducer at T1000 (Levelogger, Gold Solinst Canada Ltd., 35 Todd Rd, Georgetown, Ontario, Canada ). Atmospheric corrections were included using a STS Barologger (Solinst Canada Ltd) in the well at T1000 ( Figure 2 1). Data were downloaded from the well weekly and as a quality control procedure, water table elevations were also measured manually with a Model 102 Lase r water level well meter (Solinst, Canada Ltd). Wells T2000 (C111AE) and T500 (C111AW) were installed and operated by the SFWMD and published on DBHydro.

PAGE 57

57 Gauge adjusted Next Generation Radar (NEXRAD) rainfall data used in this study were obtained from the SFWMD. The United States National Weather Service operates two NEXRAD sites close to the study site (i.e., KBYX in Key West, FL and KAMX in Miami, FL) that provide 2 km x 2 km NEXRAD rainfall data. There are tradeoffs between rainfall estimated by rain gau ges and NEXRAD. Rain gauges (e.g., tipping buckets) provide accurate point estimates of rainfall which are acceptable for frontal related rainfall events. However, in South Florida where most of the rainfall is received in summer and summer rainfall is dom inated by conventional or tropical rainfall forming processes, rain gauges may fail to accurately represent the orientation of the rainfall front or fail to capture the entire rainfall event (Pathak, 200 1 ). On the other hand, measurement of rainfall by NEX RAD relies on the raindrop reflectivity which could be affected by factors such as raindrop size and microwave signal reflection by other particles in the atmosphere. Skinner et al. (2008) showed that the best of the two measurement methods is realized by using rain gauge or tipping bucket data to adjust NEXRAD values. Ground surface reference evapotranspiration (ET o ) was computed from micrometeorological data (i.e., solar radiation, temperature, relative humidity and wind speed) obtained from a Florida Aut omated Weather Network (FAWN; http://fawn.ifas.ufl.edu/) station located approximately 10 km northeast of the study site at the Tropical Research and Education Center, Homestead, FL. The American Society of Civil Engineers (ASCE) standardized Penman Montei th equation was used to estimate ET o values (ASCE, 2005). We assumed a crop with the following characteristics transpiring at a potential rate: crop height (0.12 m), albedo (0.23), active

PAGE 58

58 leaf area index (1.44), and well illuminated leaf stomatal resistanc e (100.8 s/m). We applied the tool REF ET (Allen, 2011) to calculate the ASCE standardized ET o from weather data. Flux due to water table evaporation may influence soil and bedrock water contents. Previous studies have shown that when canal influences are negligible, direct evaporation from the water table significantly contributes to water table declines in the Biscayne aquifer (Merrit, 1996; Chin, 2008). Two types of models are available to estimate evaporation from a water table: physically based models and empirically based models. In this study, the latter was used because the former requires detailed data such as coefficient of diffusion of water vapor through the soil and vapor pressure above the soil surface which were not collected. Empirical models simply relate water table evaporation rate to the depth of the water table below the ground surface and are used in groundwater studies (e.g., MODFLOW uses this approach; Chin, 2008). We used a model similar to that proposed by McDonald and Harbaugh (1988 ; Equation 2 5). Chin (2008) modified Equation 2 5 and obtained Equation 2 6 for south Florida conditions. ( 2 5) ( 2 6) where E is water table evaporation [mm/day], E 0 (same as ET o ) is the potential evaporation rate at the ground surface [mm/day], d is the depth of the water table below the ground surface [m], d cr is the critical depth below which evaporation ceases [m], T is

PAGE 59

59 annual average air temperature [ o C] which is approximately 25 o C in south Florida, d 0 is water table depth above which water table evaporation proceeds at potential rate i.e., at the rate similar to the ground surface evapotranspiration [m]. Chin (2008) proposed parameters d 0 and d cr in Equation 2 6 at each observation well can be estimated fr om the least squares best fit of Equation 2 7 and the parameters described as Equations 2 8 and 2 9 ( 2 7) ( 2 8) ( 2 9) Dynamic F actor A nalysis DFA uses Equation 2 10 to describe a set of N observed time series (Ltkepohl, 1991; Zuur et al., 2003; Ritter and Muoz Carpena, 2006). The goal in DFA is to keep M as small as possible while still obtaining a good model fit. Including relevant explanatory variables helps to reduce some of the unex plained variability in the observed time series ( 2 10) ( 2 11) where s n (t) is a vector containing the set of N time series being modeled (response variables), (t) is a vector containing the common trends (same units as the response variables), are factor loadings or weighting coefficients that indicate the importance

PAGE 60

60 of each of the common trends to each response variable (unitless), is a constan t level parameter for shifting time series up or down, is a vector containing explanatory variables, and are weighting coefficients for the explanatory variables (regression parameters) which indicate the relativ e importance of explanatory variables to each response variable (inverse units to convert into response variable units), and and are independent, Gaussian distributed noise with zero mean a nd unknown diagonal covariance matrix. The elements in the covariance matrix represent information that cannot be explained by the common trends or the explanatory variables. The unknown parameters and were estim ated using the Expectation Maximization (EM) algorithm that is described in Dempster et al. (1977) and Shumway and Stoffer (1982). The common trends in Equation 2 11 were modeled as a random walk (Harvey, 1989) and were predicted using the Kalman filter an d EM algorithms. The regression parameters in Equation 2 10 are estimated using the same procedure as used in linear regression (Zuur et al., 2003). DFA was implemented using a statistical package called Brodgar Version 2.5.6 (Highland Statistics Ltd., New burgh, UK). The results from the DFA were interpreted in terms of the canonical correlations ( ), factor loading ( ), regression parameters ( ) and agreement between modeled and observed soil a nd bedrock water contents (i.e., expressed as scaled frequency). The goodness of fit between modeled and observed soil and bedrock water contents were quantified using the Nash Sutcliffe coefficient of efficiency ( C eff ; Nash and Sutcliffe, 1970), the Akaik AIC ; Akaike, 1974) and the Bayesi ( BIC) C eff provides an estimate of how well a model predicts an

PAGE 61

61 observed data set, while AIC and BIC are relative measures of the goodness of fit of a statistical model. A model with the C eff closest to 1 and lowest AIC and BIC is the preferred DFA model. Cross correlations between the soil and bedrock water content time series and common trends were measured using In our study close to unity implied that the common trend was highly a ssociated with water content time series. Typically canonical correlations are classified as follows: | |>0.75, 0.5 0.75, and 0.3 0.5 as high, moderate, and weak correlations, respectively. The influence of the explanatory variables o n water content time series were quantified using the magnitude of the coefficients and their associated standard errors which were used with a t test to assess whether the response variable and explanatory variables were significantl y related. DFA was implemented sequentially by varying the number of common trends M until a minimum AIC and BIC and C eff closest to one were achieved (Zuur et al., 2003). After identifying the minimum M different combinations of explanatory variables we re introduced into the analysis until a combination of common trends and explanatory variables that resulted in the most parsimonious model with best good of fit indicators was achieved. The procedure followed here is similar to that described by Ritter et al. (2009). Soil and bedrock water content time series are autocorrelated (Kaplan and Muoz Carpena, 2011) while evapotranspiration and rainfall time series are not. For example, soil and bedrock water contents at time t will depend on antecedent soil an d bedrock water contents at time (t 1) whereas the rainfall today does not depend on rainfall yesterday. Therefore in order to relate the soil and bedrock water content time

PAGE 62

62 series and evapotranspiration and rainfall time series, we calculated a new variab le called net cumulative recharge ( R net ) using : ( 2 12) where P t is the total rainfall for day t (mm) and E ot is the potential evapotranspiration on day t (mm/day). Cumulative water table evaporation was also used instead of daily va lues. To minimize multi collinearity of explanatory variables, we used mean water table elevation instead of water table elevation at each well. Before proceeding with the DFA, multi collinearity of explanatory variables was quantified by computing varianc e inflation factors (VIFs) for each explanatory variable (Zuur et al., 2007). Simple P redictive R egression M odel for S oil W ater C ontent The simple regression model was developed from a DFA model having the minimum number of common trends required to explai n underlying common patterns in the eleven time series and explanatory variables with significant influence on modeled soil water and bedrock water content time series. To enable practical use of the simple model, DFA was performed again for the identified model using non normalized/non standardized time series. After estimating the parameters through DFA the common trends were ignored in the model to derive a simple expression relating identified significant explanatory variables and soil and bedrock wate r contents. The period from August 25, 2010 to December 2011 was used to develop the regression model while the data from January 01, 201 2 to June 30, 2012 was used to validate the new simple model. The developed simple model was then applied to predict th e impact of a 6, 9 and 12 cm increase in canal stage on soil and bedrock water contents at the study sites.

PAGE 63

63 Results and D iscussion Visual E xploratory A nalysis of E xperimental Ti me S eries Visual inspection of soil and bedrock water content time series expre ssed as SF indicates that there were some common patterns in the temporal variation of soil and bedrock water contents at the three sites (T500, T1000 and T2000) along the transect perpendicular and east of canal C111. From February 2011 to July 2011, soil and bedrock water contents gradually decreased at all monitoring elevations and all sites ( Figure 2 3). The gradual decrease in soil and bedrock water contents corresponded to the decline in canal stage and water table elevation ( Figure 2 4). The period f rom April to August was characterized by pronounced drying and wetting cycles at all sites. The wetting or spikes in soil and bedrock water contents in this period correspond to the start of the rains while the drying cycles correspond to the increasing po tential evapotranspiration during the same period ( Figure 2 4). The period from late March to July corresponds to the end of the growing season and beginning of the wet season. From August 2011 to February 2012, soil and bedrock water increased correspondi ng to stage operation criteria within the canal network that enhances water storage in the system. However, there were observed differences in temporal soil and bedrock water variability at the three monitoring sites along the transect. Site T500 which is the shallowest and closest to the canal exhibited lack of temporal variation in bedrock water content at elevations less than 0.9 m NGVD29 while soil water content at 1.0 m NGVD29 exhibited temporal variation in the same period probably due to irrigation during the growing season. Site T1000 (i.e., approximately 1000 m from canal C111) exhibited the least increase in water content between March 2011 and June 2012.

PAGE 64

64 Unlike sites T500 and T2000, the trends in soil and bedrock water contents at T1000 were not identical to the temporal variation in canal stage or water table elevation suggesting micro topography within the field might be affecting soil and bedrock water contents since this site had the highest elevation along the transect ( Figure 2 2 ). At site T 2000 (i.e., approximately 2000 m from canal C111), soil and bedrock water contents for the periods between August 2010 to March 2011 and August 2011 to February 2012 were similar characterized by small temporal variation similar to those exhibited at site T500. Sites T500 and T2000 have very similar elevation (1.1 9 and 1.2 m NGVD29 respectively) implying that topography or ground surface elevation might exert a stronger influence on temporal variation of soil and bedrock water contents compared to distance from the canal. Differences also existed at the different monitoring elevations with bedrock water content generally higher at the lowest elevation at each site. Other reasons for observed differences in water content at the different sites could be a com bination of several factors such as differences in soil surface conditions, soil and bedrock heterogeneity (specifically differences in soil water retention and unsaturated hydraulic conductivity) and differences in the environments surrounding the EnviroS can access tubes. All the hydrologic variables monitored ( Figure 2 4) exhibited seasonal variations with rainfall increasing during wet season (May to October) resulting in increased water table elevation and canal stage. ET o also increased during the wet season. In turn, decreased depth to water table and increased ET o resulted in increased E ( Figure 2 4). The water table evaporation parameters for Equation 2 6 were computed following the procedure described by Chin ( 200 8 ) in which steady declines in wate r table elevation

PAGE 65

65 particularly in the dry season when canal stage was maintained relatively constant are assumed to be caused by water table evaporation. Using data from a total of six wells (i.e., the 3 wells along transect T and 3 additional wells approx imately 1 km north of the transect) within the vicinity of the study area, we obtained an average critical depth of 1.94 m which is within the range of 1.5 to 2.9 earlier reported by Chin ( 200 8 ). We obtained a value of 0.59 m for the depth above which wate r table evaporation proceeds at the potential rate which is approximately half the average value of 1.4 reported by Chin ( 200 8 ). The water table elevations at the three monitoring sites were very similar and also corresponded to the temporal variations in canal stage on the tail water side of the spillway at S177. Response and E xplanatory V ariables Visual inspection indicated that seasonality affects temporal variation of both response variables (i.e., soil and bedrock water contents at different elevation s) and explanatory variables (i.e., ET o rainfall P water table elevation, E and canal stage). We attempted to remove seasonality effects through seasonal standardization following procedures described by Salas (1993), but this approach was abandoned sinc e it resulted in poor model fit compared to the models in which seasonal effects were assumed to be masked in the common trends (i.e., average C eff < 0.7 and C eff > 0.9, respectively). The poor model fit could be attributed to loss of information resulting from seasonal standardization. Ritter et al. (2009) also reported improved DFA model fit after back transforming refractive index data from a TDR as a surrogate for soil water content compared to seasonally standardized refractive index. To facilitate in terpretation of factor loadings and comparison of regression parameters as suggested by Zuur et al. (2004), all the time series were normalized.

PAGE 66

66 Therefore, the DFA results presented in reference to objective 1 are based on normalized time series data. Prio r to performing the DFA, multi collinearity in explanatory variables was quantified by calculating Variance Inflation Factor (VIFs) for each explanatory variable. Threshold VIF of 5 was set as the highest, high values of VIF indicate multi collinearity in the explanatory variables which makes interpretation of regression results difficult (Ritter et al., 2009). As expected there was high multi collinearity between water table elevation time series for different wells (VIFs > 30), but this was considerably r educed when mean water table elevation at the three sites was used instead (i.e., VIFs < 2). There was also high multi collinearity between headwater and tail water canal stages at S177 (VIFs > 8) implying that these two time series could not be used as ex planatory variables in the same DFA model. Mean water table elevation was also correlated to canal stage S177 (VIFs >10) probably due to the high hydraulic connectivity between C111 and Biscayne aquifer. The correlation coefficient between canal stage and water table elevation time series was greater than 0.9. Common T rends We developed the DFA model by exploring common trends and explanatory variables in relation to the 11 observed water content time series. Results of the DFA model selection are summarize d in Table 2 1. We used the AIC the BIC (which penalizes more strongly for over parameterization than the AIC ) and the C eff statistic for deciding which of the DFA models with zero explanatory variables best described the response time series. Ten was t he maximum number of common trends used to describe common variability in the 11 response water content time series. However, the goal of DFA is to minimize the number of common trends while maintaining a good model fit. Several models consisting of fewer numbers of common trends and noise

PAGE 67

67 were tested and model 4 with five common trends was determined to be the model with the minimum number of common trends required to describe the 11 response time series. Model 4 was selected since using M>5 resulted in ne gligible improvement in model goodness of fit measures while increasing the number of parameters to be interpreted. The three common trends with high ( |>0.75) to moderate (0.5< <0.75) canonical correlations partic ularly at sites T500 and T2000 are shown in Figure 2 5. Common trends 2 and 3 exhibited minor cross correlation with water content time series as measured by < 0.5 at all the sites and in the interest of brevity are not presented. Vis ually, the unexplained variation in soil and bedrock water contents described by the common trends in Figure 2 5 is similar to the seasonal variation of soil and bedrock water contents at sites T500, T1000 and T2000 for the period August 2010 to August 201 1. There was greater uncertainty as shown by a large (95%) confidence interval from August 25, 2010 to January 21, 2011 which is due to missing data for sites T500 and T1000 during this period. The first common trend exhibited high positive (| | 0.75) correlation with soil and bedrock water content time series at sites T500 and T2000 with low surface elevation (1.1 and 1.2 m NGVD29, respectively) compared to the moderate to weak correlation at site T1000 with ground surface elevation of 2.17 m NGVD29. Indicating that in addition to other factors, such as irrigation during the growing season, micro topography within the field influences temporal variations in soil water content as it governs the effect exerted by the water table.

PAGE 68

68 Relative C ontribution of E xplanatory V ariables Introducing net surface recharge, water table evaporation, and mean water table elevation or C111 canal stage to model 4 resulted in the best models (11 and 13). Inclusion of explanatory variable s in the DFA model also produced regression parameters ( ) and since response and explanatory variables were normalized, the regression parameters were used to quantify the relative influence of each explanatory variable on the modeled soil and bedrock water content time series. It is worth noting that substituting mean water table elevation in model 11 with canal stage as in model 13 resulted in AIC and BIC that were not substantially different and similar goodness of fit indicator (Ta ble 2 1). Since part of the motivation for this research was to assess the effect of canal stage management on soil and bedrock water contents, further analysis was made on model 13 because canal stage data have a more consistent record compared to water t able elevation data. At the study site, canal stage can be used as a good approximation to water table elevation due to the high permeability of the aquifer. Model 13 fitted plots are shown in Fig u res 2 6 to 2 8; these figures indicate that DFA modeling w as successfully applied to describe temporal variations in soil and bedrock water contents at all three monitoring sites and elevations ( C eff > 0.9). Results in Table 2 2 indicate that net surface recharge ( R net ) had a significant influence ( t value >2 ) on the temporal variation of soil and bedrock water contents at sites T500, T1000, and T2000 but was not significant at lower elevations at sites T1000 and T2000 as shown ( t value <2 ). The significance of R net could be attributed to rainfall ( P ) patterns in the study area in which two thirds of the P was received in the wet season (SFWMD, 2011) and these large amounts of net water input to the vadose zone are sufficient to

PAGE 69

69 maintain soil and limestone bedrock near saturation, while absence of P in the dry seas on was responsible for the dry conditions. Lack of significance at lower elevations at sites T1000 and T2000 could be attributed to heterogeneity in soils and bedrock (e.g., differences in hydraulic conductivity), and differences in surface cover which inf luence ET o Water table evaporation was found to not significantly influence temporal variation of soil and bedrock water contents ( t value <2 ) at all the sites monitored. The non significant effect of water table evaporation on soil and bedrock water con tent could be attributed to the fact that there is sufficient water for evaporation due to the shallow water table. However, the negative effect was stronger at site T1000, the negative effect is due to the fact that water table evaporation is a net loss f rom the vadose zone system. The small positive water table evaporation regression coefficient at T1000 and T2000 (Table 2) could be attributed to computational numerical errors. These results are worth highlighting given the fact that meteorological based methods for estimating ET o like Penman Monteith equation are criticized for ignoring evaporation from the shallow water table meaning they might under estimate total ET o losses. These observations could be attributed to that fact ET o in such cases is not l imited by water availability but by available energy only. C111 canal stage on the tail water side at the S177 spillway ( Figure 2 1) had the strongest influence on soil and bedrock water content temporal variations ( t value >7 ) for most sites. This findin g is significant because it confirms the hypotheses that the shallow water table and canal stage are highly connected and that canal stage can be used to predict soil water content at a given location. From a hydrologic perspective,

PAGE 70

70 these results were expe cted because in this case canal stage is used as an approximation for the shallow water table which serves as the lower boundary condition for the vadose zone and therefore regulates available storage during the rainy season. Based on the relative magnitud es of the regression coefficients (Table 2 2 ), the overall contribution of canal stage on the respective soil and bedrock water content time series is higher than that of net recharge. The factor loadings ( ) for the five common tren ds are shown in Table 2 2 these represent the influence of each common trend on the modeled soil and bedrock water content time series at the different monitoring sites and elevations. Since the time series in the DFA were normalized, the coefficients and can be compared (Zuur and Pierce, 2004). The results indicate that trend 1 was very critical for describing unexplained variation in soil water dynamics at site T2000, while common trends 2 and to a lesser extent 3 were more critical for describing unexplained variation in soil water content at site T1000. Site T500 was sufficiently described by the explanatory variables and constant level parameters given their magnitudes were larger compared to the Trends 4 and 5 had minor effects at all the monitoring sites. Overall at all the sites, compared to regression coefficients and the constant level parameters, common trends had less influence on soil and bedrock water dynamics. However, since the values of the factor loadings are not zero (i.e., they account for some unexplained variability) especially at T2000 and site T1000, this implies that the information provided by the hydrologic variables used as the explanatory variables in the DFA models only account for part of the unexplained variability in the temporal variation of the soil and bedrock water contents. Other information such as irrigation,

PAGE 71

71 differences in soil surface conditions, differences in the environment surrounding the EnviroScan access tube, and variation in soil hydraulic properties not considered in this study might account for part of the remaining unexplained variability. Predicting Soil and Bedrock Water Contents Using a Simplified Dynamic Factor Analysis Based Model To enab le practical application of the DFA model, the common trends and two of the exploratory variables included in model 13 were used in a new DFA model with non standardized time series. This new model was referred to as model 14. To further simplify model 14, we ignored the common trends to derive a simple model that predicts soil and bedrock water contents as function of net recharge and canal stage expressed as : ( 2 13) where SF(X,Z,t) is the SF at distance X from the canal, at elevatio n Z and time t other terms are previously described and varies with elevation and distance from the canal. The coefficients for Equation 2 13 at all the sites and monitoring elevations are obtained from Table 2 3. The C eff in Table 2 3 are calculated ba sed on Equation 2 13 with common trends removed. As expected, performance of the simple model ( Equation 2 13) was lower as shown by the reduction in C eff (Table 2 3 and Fig ures 2 9 to 2 10) compared to the DFA models that include common trends particularly for site T1000. Since factor loadings are not zero for all the trends (Table 2 3), this suggests that the explanatory variables (net recharge and canal stage) used in the DFA model are not sufficient to explain all the observed variations in the soil and bedrock water content time series. This is particularly true at site T1000 which is affected by 4 out of the 5 common trends. Common trend number 2 appears to affect all the sites, it probably

PAGE 72

72 masks common variation such seasonal changes in rainfall, evap otranspiration and canal stage. Other common trends had minor effects at all the other sites particularly at site T1000. The difference in response at site T1000 could be attributed to differences in elevation as shown in Figure 2 2 site T1000 has a high er surface elevation and hence larger depth to water table. The results in Table 2 3 also underscore the point that the effect of canal stage is stronger at low elevation sites T500 and T2000 compared to T1000. Thus, proper interpretation of modeling resu lts in this area requires accurate quantification of micro topography. Model performance ranged from good at sites T500 and T2000 to poor at site T1000 with root mean square error (RMSE) ranging from 0.005 to 0.01. Fig ures 2 9 to 2 10 show model performanc e during the calibration and validation periods, after removing the common trends, it can be seen that the simple model misses the peaks but is able to generally predict the temporal variation in soil and rock water content. The simple model ( Equation 2 13 ) could be improved by using location specific water table elevation since canal stage is simply a good approximation of the mean water table elevation. Another simple sigmoidal regression model to predict soil and bedrock water contents from canal stage p roposed by Kaplan et al. (2010a) was tried but later abandoned due to lower C eff (i.e., averaging 0.2). Th eir approach is based on the physical concept of drain to equilibrium. However, for our study site this condition was hard to achieve since during the dry season irrigation was taking place and in the rainy season there was frequent rainfall hence by removing data points corresponding to rainfall or irrigation, very few data points were left to develop a useful sigmoidal model for predicting soil and be drock water content from canal stage.

PAGE 73

73 Assessing the impact of proposed operational changes in C111 canal stage management on soil and bedrock water content The low lying agricultural areas east of canal C111 are anticipated to experience the greatest impac t from the proposed changes in C111 stage operation (i.e., canal stage increases of 6, 9, and 12 cm); a simple DFA based regression model (Equation 2 13 ) was proposed to predict the soil and bedrock water contents as a function of canal stage. We considere d the period from January 01, 2012, to June 30, 2012 for the analysis. Increases in canal stage were computed by simply adding the proposed incremental rises in canal stage to the daily canal stage recorded at S177T while P and ET o from the original datase t were not changed. The results from using this simplified DFA based model ( Fig ures 2 11 and 2 12) indicate that the proposed increases in canal stage were predicted to have changes in daily mean SF for the study period (i.e., which is used as a surrogate for soil and bedrock water contents) of <1% at all sites and all elevations monitored. The range in daily SF differences was 0.065 to 0.024 and 0.075 to 0.041 at sites T500 and T2000 respectively, which indicates that the simple model over predicted an d under predicted SF on certain days during the study period. However, note that the daily differences in SF are not substantially large; this may be attributed to already high values of soil and bedrock water contents observed in the area. On an event ba sis the potential to flood or saturate the root zone would depend on the size of the storm and storm contingency planning for lowering of canal stage in anticipation of heavy storms. Since we showed using DFA that soil and bedrock water contents were sign ificantly affected by canal stage and net recharge.

PAGE 74

74 The simple model used in this evaluation was more accurate at sites T500 and T2000 and therefore results at these two sites would be considered with less uncertainty. Soil and bedrock water responses to incremental raises in canal stage were not computed for site T1000 since results at this site would be considered less accurate (greater uncertainty) because model performance was very poor at this site. S oil and bedrock water contents were more noticeable at the highest elevation ( Figures 2 11 and 2 12 ) However, at the lowest elevations monitored the difference between mean SF before and after all increments was zero at T500. These observations could be attributed to the fact that low elevation sites are normally close to saturation. For example, at site T500 (0.7) when water elevation was above the sensor (implying saturated conditions), SF was recorded as 0.786 compared to average SF of 0.775 for the study period meaning small changes in water table may not result in substantial changes in soil water content since the pores are already near saturation. Note that the simple model developed should be applied with the following limitations in mind. The model does not account for water input from irrigation a nd therefore would under predict soil and bedrock water content during the growing season, the model also uses canal stage as an approximation for water table elevation at a specific location although the two are usually close there may be deviations espec ially after large rainfall events, it ignores water content drivers that were masked in the common trends, and lastly the simple model ignores the effect of E which might vary based on micro topography within the field as well as differences in land surfac e cover conditions. Finally, although the simplified DFA based model is empirical in nature, the results suggest it can be used as a preliminary tool to relate the potential

PAGE 75

75 impacts of surface water management decisions on soil and bedrock water contents i n low lying farmlands adjacent to canal C111. This is because during the duration of the study, we were able to capture a wide range of variation in canal stage and water table elevation e.g., on June 10, 2011 we recorded canal stage and groundwater table elevation of 0.14 m NGV29 which is lower than the optimum design stage of 0.6 m for the reach of C111 between S 18C and S177 under current canal stage operational criteria. O n October 09, 2011 we recorded canal stage and groundwater levels as high as 0.9 a nd 1.02 m NGVD29 which is close to the level supposed to trigger the spillway to open at S177 under current operational criteria. Conclusions The response of soil and bedrock water contents to incremental raises in canal stage proposed under the C111 spre ader canal project whose goal is to restore the hydrology of ENP while maintaining flood protection in the adjacent agricultural areas was investigated using DFA. The study objectives were to use DFA to identify the important factors driving temporal varia tion in soil and bedrock water content above the shallow water table at the study site, develop a simple model for predicting soil water content as a function of canal stage and assess the effect of the proposed incremental raises in canal stage on soil an d bedrock water contents. Five was the minimum number of common trends required to account for the unexplained variation in the eleven observed soil and bedrock water content time series while producing an acceptable model fit. Introduction of explanatory variables i.e., net recharge, water table evaporation, and canal stage or water table elevation to the DFA model resulted in lowering AIC and BIC values while C eff values did not substantially change. Evaluation of the regression coefficients indicated tha t net recharge and canal stage had significant

PAGE 76

76 effects on temporal variation of soil and bedrock water contents while the effect of water table evaporation was non significant. Based on the magnitude of the regression coefficients, canal stage had the grea test influence on the temporal variation of soil and bedrock water contents at all elevations and distances from the canal at the locations monitored. The effect of canal stage and mean water table elevation in the DFA model was similar confirming the high hydraulic connectivity between the canal and Biscayne aquifer. Based on the high connectivity between surface water in the canal and Biscayne aquifer, a simple DFA based regression model (DFA model in which the common trends were removed), was developed to predict soil and bedrock water contents as a function of canal stage and net recharge at various elevations. The performance of the simplified regression model was described as good to acceptable at sites with low elevation (i.e., water table elevation within 1m from the ground surface) and poor at the location at with water table depth greater than 1.5 m. These findings highlight the effect of micro topography within the field on soil water content. The study also revealed that factor loadings were not zero for all the common trends suggesting that the explanatory variables (net recharge and canal stage) used in the DFA model are not sufficient to explain all the observed variations in the soil and bedrock water content time series. The effect of the pr oposed 3 incremental raises in canal stage on soil and bedrock water content was simulated using the developed simple DFA based regression model for a total of 181 days beginning January 01, 2012. The results based on the data collected indicate that the p roposed raises in canal stage would result in negligible changes in average soil and bedrock water contents at low elevations

PAGE 77

77 monitored in this study based. Changes in soil water content near the ground surface were more noticeable. The DFA based regressi on model developed is limited in its prediction ability to the range of canal elevations and net recharge by which it was developed. The uncertainty in predictions could be minimized by continuously updating the regression coefficients and constant level p arameters as more data on response and explanatory variables are collected. The results of the regression model could be further evaluated using physically based modeling approaches. The approach used in this study could be applied to any system in which d etailed physical modeling would be limited by inadequate information on parameters or processes governing the physical system.

PAGE 78

78 Table 2 1 Dynamic Factor Analysis (DFA) models tested based on the following goodness of fit measures: AIC, BIC and C eff Mode l No. of common trends Explanatory variables No. of parameters AIC 1 BIC 2 C eff 3 Step I (DFA model with K=0) 1 2 None 98 2690.5 2041. 8 0.68 2 3 None 107 4654.2 3945.9 0.84 3 4 None 115 5830.2 5068.9 0.88 4 5 None 122 6901. 5 6093.8 0.97 5 6 No ne 128 7028. 8 6181.4 0.97 6 8 None 137 7263.9 6357.0 0.97 Step II (DFA model with K>0) 7 5 R net 4 133 7018.6 6138. 2 0.97 8 5 R net E 5 144 7797.5 6844. 3 0.98 9 5 S177T 6 133 734 1 0 6460.5 0.97 10 5 S177T, R net 144 7542. 7 6589.4 0.97 11 5 R net E, MWT 7 155 8052.4 7026. 4 0.98 12 5 MWT, R net 144 7444.0 6490. 8 0.97 13 5 R net E, S177T 155 7922. 4 6896. 3 0.98 1 AIC Akaike information criterion 2 BIC Bayesian Information Criterion 3 C eff Nash Sutcliffe coefficient calculated based a ll the nine observed time series 4 Cumulative net surface recharge 5 R net Cumulative water at table evaporation 6 S177T Canal stage in C111 7 MWT Mean water table evaporation

PAGE 79

79 Table 2 2. Dynamic Factor Analysis results for model 13 with 5 common trends and 3 explanat ory variables s n C eff 1 T500 (1.0) 0.05 0.02 0. 04 0.02 0.03 0.28 (0.6) 0.34 (6.9) 0.00 (0.0) 0.24 (8.8) 0.93 T500 (0.9) 0.05 0.06 0.03 0.04 0.05 0.37 (0.5) 0.24 (3.5) 0.14 ( 0.3) 0.29 (8.3) 0.94 T500 (0.8) 0.04 0.06 0.03 0.03 0.04 0.34 (0.6) 0.20 (3.2) 0.17 ( 0.5) 0.22 (7.5) 0.90 T500 (0.7) 0.03 0.02 0.01 0.00 0.01 0.13 (0.6) 0.18 (7.1) 0.09 ( 0.7) 0.13 (9.1) 0.90 T1000 (1.97) 0.04 0.16 0.13 0.02 0.00 0.95 (0.9) 0.47 (3.1) 0.53 ( 0.7) 0.62 (8.7) 0.85 T1000 (1.87) 0.04 0.20 0.11 0.01 0.01 0.82 (0.8) 0.38 (2.1) 0.61 ( 0.8) 0.70 (8.5) 0.8 1 T1000 (1.77) 0.01 0.50 0.01 0.00 0.00 0.00 (0.0) 0.44 (1.1) 0.23 (0.1) 0.77 (4.6) 0.67 T2000 (1.11) 0.10 0.04 0.06 0.07 0.01 0.07 (0.1) 0.13 (2.0) 0.05 (0.1) 0.50 (11.6) 0.99 T2000 (1.01) 0.13 0.05 0.06 0.02 0.06 0.09 ( 0.1) 0.03 (0.3) 0.03 (0.0) 0.68 (13.2) 0.90 T2000 (0.91) 0.17 0.03 0.06 0.01 0.01 0.12 ( 0.1) 0.05 (0.4) 0.22 ( 0.3) 0.71 (12.4) 0.93 T2000 (0.81) 0.16 0.04 0.03 0.02 0.02 0.31 ( 0.3) 0.08 (0.8) 0.02 (0.0) 0.46 (8.8) 0.96 Factor loading corresponding to common trend 1to 5 and observation, n Constant level parameter in dynamic factor model with associated t value in parenthesis Regression parameter corresponding to the 3 explanatory variables (net recharge [R net ], water table evaporation [E], and canal sta ge in C111 [C111stage]) with associated t value in parenthesis C eff is Nash Sutcliffe coefficient 1 Site name nomenclature; T is refers to transect name T, number refers to distance from canal and number in parentheses refers to elevation NGVD29 m n number of observations

PAGE 80

80 Table 2 3. Dynamic Factor Analysis results for model 14 with 5 common trends and 2 explanatory variables implemented with non standardized time series s n 1T500 (1.0) 0.003 0.000 0.000 0.000 0.000 0.72 0.14 0.06 0.73 0.70 T500 (0.9) 0.001 0.004 0.000 0.000 0.0 00 0.72 0.11 0.04 0.61 0.62 T500 (0.8) 0.001 0.004 0.000 0.000 0.000 0.75 0.09 0.02 0.51 0.56 T500 (0.7) 0.001 0.002 0.002 0.000 0.000 0.76 0.07 0.01 0.81 0.74 T1000 (1.97) 0.003 0.005 0.002 0.000 0.001 0.73 0.10 0.02 0.61 0.15 T1000 (1.87) 0.002 0.003 0.001 0.000 0.001 0.74 0.05 0.01 0.51 0.13 T1000 (1.77) 0.001 0.003 0.001 0.002 0.000 0.77 0.02 0.00 0.25 0.11 T2000 (1.11) 0.000 0.003 0.002 0.000 0.000 0.76 0.08 0.06 0.70 0.61 T2000 (1.01) 0.000 0.003 0.000 0.000 0.001 0.76 0 .05 0.05 0.60 0.67 T2000 (0.91) 0.000 0.002 0.000 0.000 0.001 0.77 0.03 0.04 0.67 0.63 T2000 (0.81) 0.001 0.001 0.000 0.000 0.001 0.80 0.02 0.02 0.65 0.61 Constant level parameter in dynamic factor model Regression parameter corresponding to the 2 explanatory variables (net recharge [R net ], and canal stage in C111 [C111stage]) 1 Nash Sutcliffe coefficient are calculated after ignoring common trends 2 Nash Sutcliffe coefficient during validation

PAGE 81

81 Figure 2 1 Map of C 111 basin showing farmlands east of canal C 111 Taylor Slough, Florida Bay and three monitoring sites along transect perpendicular to C111 canal.

PAGE 82

82 Figure 2 2 Map of C 111 basin showing farmlands east of canal C 111 Taylor Slough Florida Bay and three monitoring along transect perpendicular to C111 canal.

PAGE 83

83 Figure 2 3 Temporal variation in scaled frequency (i.e., soil and bedrock water contents ) at three sites (i.e., T500, T1000 and T2000 with soil and bedrock water contents m onitored at different elevations using EnviroScan probes) along a transect perpendicular to C111 on the tail water side of the spillway at structure S177 during the period August 2010 to June 2012.

PAGE 84

84 Figure 2 4. Temporal variation in hydrologic factors e valuated for their influence on soil and bedrock water contents at the study site during the period August 2010 to June 2012.

PAGE 85

85 Figure 2 5. C ommon trend s with 95% confidence interval describing unexp lained temporal variation in scaled frequency as a sur rogate for soil and bedrock water content and the canonical correlation for quantifying the correlation

PAGE 86

86 between water time series and the common trends in the nomenclature for site names the number represents distance from the canal in m, and the numbers in the parenthesis represent elevation NGVD 29 m. Figure 2 6 Fitted Dynamic Factor Model (DFM) and observed temporal variation in scaled frequency (used as a surrogate for soil and rock water ) in gravely loam

PAGE 87

87 soils and limestone bedrock at a site locat ed 500 m along a transect from C111 and the numbers in the parentheses indicate elevations. Figure 2 7 Fitted Dynamic Factor Model (DFM) and observed temporal variation in scaled frequency (used as a surrogate for soil and rock water ) in gravely loam soils and limestone bedrock at a site located 1000 m along a transect from C111 and the numbers in the parentheses indicate elevations.

PAGE 88

88 Figure 2 8 Fitted Dynamic Factor Model (DFM) and observed temporal variation in scaled frequency (used as a surroga te for soil and rock water ) in gravely loam soils and limestone bedrock at a site located 2000 m along a transect from C111 and the numbers in the parentheses indicate elevations.

PAGE 89

89 Figure 2 9. Performance of a simple model for predicting scaled frequenc y ( used as a surrogate for soil and bedrock water content ) as a function of canal stage and net recharge at specific elevations in parentheses NGVD29 at a site located 500 m along a transect from C111.

PAGE 90

90 Figure 2 10. Performance of a simple model for pre dicting scaled frequency (used as a surrogate for soil and bedrock water content ) as a function of canal stage and net recharge at specific elevations in parentheses NGVD29 at a site located 2000 m along a transect from C111.

PAGE 91

9 1 Figure 2 11. Boxplots show ing soil and rock water content as measured using scaled frequency at site T500 before and after 6, 9 and 12 cm increase canal at structure S18C along C111.

PAGE 92

92 CHAPTER 3 SENSITIVITY ANALYSIS AND PARAMETER ESTIMATION FOR AN APPROXIMATE ANALYTICAL MODEL OF CAN AL AQUIFER INTERACTION APPLIED IN THE C 111 BASIN Introduction Surface water in streams, canals, and rivers infiltrates into hydraulically connected aquifers during rising flood stage and is released back into surface water bodies during flow recession. T his interaction between groundwater and surface water bodies occurs in virtually all types of landscapes (Winter et al., 1998). However, it is particularly critical in low elevation coastal watersheds (e.g., southeastern Atlantic and Gulf coasts of the Uni ted States) with shallow water tables due to the high risk of groundwater flooding coastal landscape has a shallow water table aquifer (Biscayne aquifer) that is hydraulically connected to an extensive canal network primarily constructed for flood control (Chin, 1990; Genereux and Guardiario, 1998). As part of flood control management in south Florida, pre storm contingency planning normally involves artificially lowering the water table to cre ate sufficient storage for the forecasted storm (Bolster et al., 2001). Canal stage lowering has to consider other fresh water uses such as ecosystem restoration in the Everglades and control of salt water intrusion. Optimizing flood control and water man agement decisions requires continuously improving our knowledge of factors such as canal bed conductance, aquifer specific yield, and aquifer saturated thickness as well as horizontal and vertical hydraulic conductivi ti es which may be used to characterize and predict aquifer responses to canal Co pyright for this chapter of the dissertation belongs to Trans. ASABE Kisekka et al. 2013b

PAGE 93

93 stage management. There are many ways of determining these physical parameters e.g., using standard pumping tests and slug tests for quantifying hydraulic conductivity and specific yield (Freeze and Cherry, 1979) and canal drawdown tests for assessing aquifer hydrogeological parameters and canal bed conductance (Bolster et al., 2001). For regional studies, pumping and slug tests are expensive, in addition they are particularly challenging for very permeable aquifers such as the Biscayne aquifer due to large sizes of pumps and water conveyance pipes required to produce a large enough drawdown to be accurately measured (Fish and Stewart, 1991). Another approach that has been used for characterizing hydrogeological param eters in surface water groundwater interaction problems is parameterization of numerical or analytical models. In this study we focus on the later, since numerical models when applied within a global sensitivity and parameter estimation framework could be computationally expensive (Vesselinov et al., 2012). Many analytical models have been developed to describe the interaction between surface water bodies such as canals and groundwater aquifers (e.g., Hall and Moench, 1972; Serrano and Workman, 1998; Zlotn ik and Huang, 1999; Moench and Barlow, 2000; Lal, 2001; Hantush, 2005). The main difference among these analytical models is the simplifying assumptions made in deriving a solution to the groundwater flow equation. In most analytical models, aquifer respon se to arbitrary stage and recharge fluctuations is simulated using the convolution integral (Olsthoorn, 2008). For this study the approximate analytical solution developed by Barlow and Moench (1998) was selected because it is more generalized for many can al aquifer configurations and boundary conditions and also accounts for the effect of both arbitrary canal stage and

PAGE 94

94 recharge variations. In addition, computer programs are available to facilitate implementation of the generalized unit step response functi on with convolution integral (available at: http://water.usgs.gov/ogw/staq/ ). Examples of prior studies involving estimation of hydrogeological parameters using analytical models of canal aquifer interaction include Bolster et al. (2001), Lal (2006), and Ha et al. (2007). Bolster et al. (2001) determined Biscayne aquifer specific yield using data from canal drawdown experiment and the analytical model developed by Zlotnik and Huang (1999) for partially penetra ting streams with low permeability bed sediment layer in the absence of recharge. Lal (2006) determined bulk aquifer and canal resistance dimensionless parameters for Biscayne aquifer using a coupled canal aquifer analytical model. Ha et al. (2007) applied an analytical model of river aquifer interaction to estimate aquifer diffusivity and river bed resistance in the Mangyeong R iver floodplain, South Korea. In prior investigations related to Biscayne aquifer, sensitivity analysis was implemented using local sensitivity analysis techniques. Local sensitivity analysis is limited as (1) parameter interactions are not considered ; (2) parameter importance is assessed only in the vicinity of mean parameter value, and (3) it is unreliable for nonlinear models (Frey and Patil, 2002). Global sensitivity analysis techniques overcome these limitations and provide more information on the response of linear and nonlinear model output to variations in model input factors. Global sensitivity analysis methods explore the e ntire parametric space of the model simultaneously for all model input factors which allows them to capture first and higher order parameter effects. Different global sensitivity analysis methods may be selected depending on the objective of the analysis, number of uncertain input factors,

PAGE 95

95 and computing time required for a single forward model simulation ( Saltelli et al., 2000 ; Muoz Carpena et al., 2007). Saltelli et al. (2004) proposed that robust statistical frameworks for model evaluation should be based on global analysis techniques meeting the following requirement: 1) are model independent i.e., can work with various models without the need for modification, 2) contain a screening method for qualitative identification of a subset of important model inp ut factors, 3) contain a method that can quantitatively decompose model output variance in terms of first and higher order input factor effects, and 4) can allow for uncertainty analysis through construction of output probability density function (PDFs). T he method of Morris (Morris, 1991) provides a robust screening method and yet it requires few model simulations while variance based higher order or interaction input fa ctor effects (Sobol, 1993). Morris (1991) proposed an effective screening sensitivity measure for identifying important parameters in models that have many parameters. The Morris (1991) method aims at identifying input factors whose effect on the model output is negligible, linear/additive and nonlinear/involving interactions with other factors. The original Morris (1991) method is based on computing for each input factor using the one factor at a time (OAT) approach a number of incremental ratios called elementary effects (EE). S tatistics are then calculated from the distributions of EE for each input factor which are used to infer model output sensitivity to different parameters. Campolongo et al. (2007) observed that the sampling method employed in the original Morris (1991) method which is based on random sampling of the input factor space could lead to limited or non optimum coverage of the input factor space particularly for models with large

PAGE 96

96 number of input factors. Campolongo et al. (2007) then pro posed an improved sampling strategy which aims at better scanning of the input factor space without increasing the number of model executions. The philosophy behind the improved sampling strategy was to select r trajectories in such a way as to maximize th eir dispersion in the input factor space. The improved sampling strategy guarantees a global dispersion of the selected trajectories but comes at a high er computational cost. Ind ices (TSI).TSI of a given parameter includes main effects or first order effects and all interaction effects involving a particular parameter (Sobol, 1993; Chan et al., 1997; y indices is decomposing the input out relationship (i.e., model output) into summands of increasing dimensionality (Chan et al., 1997; Saltelli et al., 2000). The resulting equation is complex and is solved using Monte Carlo numerical integration to obtai n sensitivity indices. One of the draw backs of this method is computational cost especially for models with many uncertain input factors. Another key component of model evaluation is parameter estimation or model calibration using measured system response s. Beven (2006) demonstrated that there are multiple models structures and optimum parameter sets that could be used for simulating a hydrologic system. This phenomenon is known as equifinality and it arises from the fact that models are imperfect represen tations of natural systems due to various sources of uncertainty. Several methods have been developed for parameter uncertainty and among the most popular is the General Likelihood Uncertainty Estimator (GLUE). The GLUE methodology rejects the idea of one single optimal

PAGE 97

97 solution and adopts the concept of equifinality of model input factors and model structure; thus, there are multiple model structures and input factors that can be used to simulate a natural system (Beven and Binley, 1992). With regard to t he uncertain model parameters, in the GLUE analysis the prior set of model parameters is divided into a set of acceptable solutions and another set of non acceptable solutions. The degree of membership to either set is determined by assessing the extent to which model simulations fit observed data, which in turn is determined by a subjective likelihood function, e.g., the Nash Sutcliffe coefficient of efficiency (NSE). With regard to assessing parameter uncertainty, the outputs from GLUE are posterior PDFs and cumulative density functions (CDFs) describing parameter statistics. The goal of this study was to use global sensitivity and global parameter estimation techniques with an analytical model of canal aquifer interaction to better characterize model para meters influencing the exchange of water between canal C 111 and the Biscayne aquifer in south Florida. The objectives were to: (1) apply the Morris screening technique using two sampling approaches to identify a subset of the parameters to which model out global sensitivity analysis technique on a subset of parameters obtained from Morris to quantify first order (only due to given parameter) and total effect (parameter and its interactions with other parameters) sensitivity indices, and (3) apply the GLUE methodology to obtain values of parameter sets which produce acceptable results (i.e., parameters that result in the closest agreement between predicted and measured water table elevation measured usi ng the NSE ).

PAGE 98

98 Materials and Methods Study Area The study area was an agricultural area approximately 17 km 2 located within the C 111 basin in southern Miami Dade County, Homestead, Florida, United States ( Figure 3 1). The hydrogeologic system at the study site consists of the Biscayne aquifer bordered by two canals C 111 and C 111E separated by a distance of approximately 3.5 km ( Figure 3 1) and managed by the South Florida Water Managements District (SFWMD). Biscayne aquifer is a highly permeable unconfin ed aquifer with hydraulic conductivities reported to exceed 10,000 m/day and serves as the principle source of drinking water for over 3 million people in Miami Dade, Broward, and southern part of West Palm counties in southeast Florida (Genereux and Guard iario, 1998). The Biscayne aquifer is wedge shaped increasing in thickness from the western boundary of Miami Dade and Broward Counties to a thickness of approximately 60 m near the coast. At our study site, the aquifer consists of two formations: the Miam i Limestone Formation and the underlying Fort Thompson Formation (Fish and Stewart, 1991). Aquifer thickness at our study site has not been measured but Genereux and Guardiario (1998) reported a thickness of 13.6 m for a nearby site west of our current stu dy site with roughly one third accounted for by the Miami Limestone formation. Canal C 111 was constructed in 1967 as the principle flood control canal for south Miami Dade County and partially penetrates the Biscayne aquifer to a depth of approximately 5 m (i.e., 4 m through the Miami Limestone formation and 1 m into the Fort Thompson Limestone for mation). Flow in C 111 is south towards Florida Bay and topography is essentially flat ranging between 1.1 to 2.2 m National Geodetic Vertical Datum (NGVD) 29. The width of the canal increases from north to south with an

PAGE 99

99 average width of approximately 29 m at the S 177 spillway. Currently very little is known about hydraulic properties of canal bed sediment in the lower C 111; however, several studies have docume nted the presence of a low permeability canal bed sediment layer a mixture of carbonate mud and natural organic matter in several canals within the C 111 basin (Chin, 1991; Genereux and Guardiario, 1998; Merkel, 2000). Hydrological Data Monitoring Data fr om six groundwater observations wells were used ( Figure 3 1; Table 3 1). The observation wells were separated into two groups: group 1 with wells VC1, VC2, and AK5 and group 2 with wells AK6 (referred to as T1000 in chapter 1) C 111AE (T2000) and C 111AW (T500) ( Figure 3 1). Each group was assigned a different canal stage based on its proximity to the headwater or tail waters of the S 177 spillway ( Figure 3 1). Group 1 was assigned to the C 111 headwater canal stage while group 2 was assigned to the C 111 tail water canal stage. Observation wells C 111AE and C 111AW were constructed and maintained by the SFWMD while the other four sites (AK5, AK6, VC1, and VC2) were constructed and maintained by University of Florida (UF) IFAS. The UF wells were construct ed using a 50.8 mm (nominally 2 in) PVC casing which was inserted into a 101.6 mm bore of depth 6 m. The screen size was 0.254 mm (nominally 0.1 in) and a length of 1.5 m. The well was back filled using 20/30 silica sand filter pack up to a depth of 60 cm above the screen. A 2.5 m layer of 30/65 fine sand was placed above the filter pack. A 6 m steel rod was used by gently dropping it into the annular gap to ensure the space was uniformly backfilled. The well was filled using Portland type I cement grout to a depth of 35 cm below the ground surface, the top of the well was completed with 40 cm cast iron manhole. UF wells were equipped with level loggers (Levelogger, Gold Solinst Canada Ltd., 35 Todd Rd, Georgetown,

PAGE 100

100 Ontario, Canada) to record water table elev ation every 15 minutes although daily averages were used in the modeling. Atmospheric corrections were accounted for using a STS Barologger (Solinst Canada Ltd) in well AK6 ( Figure 3 1). Data were downloaded from the well weekly and as a quality control p rocedure, water table elevations were also measured manually with a Model 102 Laser water level well meter (Solinst, Canada Ltd). Elevations at the top of the well manholes were measured using a laser level with reference to a bench mark with elevation 1.2 1 m NGVD29 near well C 111AE. Details about the construction of the SFWMD wells can be obtained at http://www.sfwmd.gov/dbhydroplsql/show_wilma_info.report_process Water tabl e elevation data for wells C 111AE and C 111AW were processed by SFWMD and published on their online Environmental database DBHydro ( http://www.sfwmd.gov/ dbhydroplsql/show_dbkey_i nfo.main_menu ). In south Florida most of the rainfall is received from the end of May to the beginning of November and is dominated by conventional or tropical rainfall forming processes. Under such rainfall forming processes, tipping buckets may fail t o accurately represent the orientation of the rainfall front or fail to capture the entire rainfall events (Pathak, 200 1 ). To minimize the uncertainty associated with spatial variability of rainfall in south Florida, gauge adjusted NEXRAD (Next Generation Radar) rainfall data were used. Skinner et al (2008) showed that both NEXRAD data and point tipping bucket measurements had limitations but the best of the two measurement methods was realized by using rain gauge data to adjust NEXRAD values. Gauge adjust ed NEXRAD rainfall data on 2x2 km grid were obtained from the SFWMD.

PAGE 101

101 Ground surface potential evapotranspiration (ET o ) was computed from micrometeorological data obtained from a Florida Automated Weather Network (FAWN; http://fawn.ifas.ufl.edu/) station lo cated approximately 10 km northeast of the study site at the Tropical Research and Education Center, Homestead, FL. The ASCE standardized Penman Monteith equation and the REF ET tool by Allen (2011) were used to estimate ET o values. Canal stage data were measured at the S 177 spillway. Daily headwater and tail water canal stage data were used. Canal stage data were measured by the SFWMD and are publically available at the online environmental database called DBHydro ( http: //www.sfwmd.gov/dbhydroplsql/show_dbkey_info.main_menu ). Analytical Model The governing equation for two dimensional groundwater flow (in vertical plane) in a water table aquifer is expressed as Equation 3 1 (Ba rlow and Moench, 1998). The initial condition is expressed as Equation 3 2. Eq uation 3 3 represents right boundary condition (BC) as the aquifer extends to infinity. Eq uation 3 4 represents a head dependent boundary flux at the canal aquifer interface. Eq u ation 3 5 represents boundary at the water table while Equation 3 6 represents a no flow bottom BC. (3 1) (3 2) (3 3) (3 4)

PAGE 102

102 (3 5) (3 6) w here K x and K z are horizontal and vertical hydraulic conductivity [m/day], S s specific storage [day 1 ], x is distance in the horizontal direction x o
PAGE 103

103 Eq uation 3 7 is combined with the convolution integral ( Equation 3 8) to predict water table responses to arbitrary changes in canal stage and recharge (difference between rainfall and ET o ). The discretized convolution equatio n was expressed by Barlow and Moench (1998) as: (3 8) (3 9) w here j is the upper limit of time integration [ ], k is the time step [ ], is the time step size [days], (k 1) is the time rate of change of canal stage and recharge [m/day], F(k 1) is canal stage or recharge at time step (k 1) and F(k) is canal stage or recharge at time step k Biscayne aquifer was conceptualized as a semi infinite water table aquifer with a thin vadose zone into which water instantaneously entered or exited. The aquifer was also assumed to be bordered by a fully penetrating C 111 canal with a semi pervious sediment layer. Time step size was one day. Since t he solution obtained from Equation 3 8 is in the Laplacian domain, it is inverted back to real time domain using Stehfest algorithm (Stehfest, 1970). The software STWT1 developed by Barlow and Moench (1998) was used to implement the computation of water ta ble head. To simulate water table response to canal stage and recharge variations using the STWT1 model requires a total of eleven input factors (refers to parameters and inputs). For the simulation we considered parameters that we did not measure to be uncertain ( K x K z S y b K s d x o ) while the inputs ( x, h i recharge, well screen length ) that were measured were considered certain. Canal bed sediment hydraulic conductivity and thickness are included in a single parameter called canal lea kance ( Equati on 3 4).

PAGE 104

104 The two STWT1outputs are water table head and canal seepage. However, sensitivity analysis and parameter estimation were based only on water table head since canal seepage was not measured. Global Sensitivity Analysis Two global sensitivity anal ysis (GSA) methods were implemented: (1) parameter screening using Morris m ethod and (2) variance based global sensitivity analysis uncertain parameters were constructed usi ng data from literature, 2) input parameter sets were obtained by sampling the multivariate input distributions according to the sampling for quantitative determinati on of sensitivity indices), 3) STWT1 model simulation was executed for each input parameter set, 4) using model output (i.e., NSE) for all parameter sets, Morris sensitivity analysis was performed in order to obtain qualitative ranking of parameters, and 5 ) using a subset of critical parameters identified in step 4, steps 2 to 4 were repeated and quantitative first order and higher order analysis was implemented using t he software package SimLab v2.2 (SimLab, 2004). We interfaced SimLab with STWT1 using a program written in Matlab (R2012a and statistics of the uncertainty model input param eters in Table 2 together with the sampling method selected, a sample input file was generated. A sample input file is a matrix comprising of multiple input parameter sets obtained from random sampling of probability distributions. The Matlab interface pro gram was used to automatically execute the model for each parameter set and to produce outputs in the desired SimLab

PAGE 105

105 format for post processing i.e. sensitivity analysis. For each simulation, the Matlab program was also used to calculate NSE (between simul ated and measured water table elevation) and Root Mean Square Error (RMSE) which were used as the model outputs in the sensitivity analysis. Morris method The finite distribution of EE associated with each input factor ( F i ) is obtained by randomly sampli ng the model input factor space For each input factor, Morris (1991) proposed two sensitivity measures, which assesses the overall effect of the factor on (i.e., nonlinear effects). To estimate r elementary effects from F i distribution of each input factor, using a design that constructs r trajectories of ( k+1 ) points in the input factor space providing k elementary effects one for each factor. In the ori ginal Morris (1991) method, the number of model executions N is computed as: (3 10) Generally in the original Morris method r is taken in the range of 4 to 10 (Saltelli et al ., 2009). For this study r of 8 and k of 6 were used resul ting in 56 model simulations. Morris (1991) proposed plotting these two measures on a Cartesian plane to aid interpretation. Campolongo et al (2007) suggested using absolute values of elementary effects, to avoid the cancelling effects of opposite signs in case of non monotonic models which was implemented in SimLab. To compare perf ormance of the two Morris sampling methods, we repeated Morris screening using the improved sampling strategy proposed by Campolongo et al

PAGE 106

106 (2007). The improved sampling strategy is implemented by initially generating a high number of Morris trajectories ( M ~500 1000), followed by choosing the best r trajectories (e.g., r =10) with the greatest spread within the input factor space. A quantity D representing the sum of distances between couples of trajectories belonging to the same combination is calculated fo llowing procedure in Campolongo et al (2007). D is calculated for all possible combinations of r trajectories which results in high computational cost especially for large models. A combination trajectories corresponding to the global maximum D is selecte d as the best set of r trajectories. To overcome the computational cost associated with the improved sampling strategy, Ruano et al ( 2012) proposed a sampling strategy that considerably reduces the computation cost required to select optimum r trajectorie s out M by developing a procedure that does not take into account all the possible combinations of r trajectories but gets a combination of r trajectories that are as close as possible to the highest spread ones. However, this sampling strategy does not gu arantee that the final trajectories selected represent maximum distance between them but ensures that the distances are at least locally maximized. Since the issue of computational cost was not an issue in this study given the small size of the STWT1 model we applied the Campolongo et al (2007) improved sampling strategy (with r =10) that ensures a global distance D for the selected optimum set of trajectories. The improved sampling strategy was implemented using Matlab algorithms for Morris sensitivity ana lysis developed by Saltelli et al (2008) and available online at http://sensitivity analysis. jrc.it/software/index.htm

PAGE 107

107 method The difference between first order and TSI was used as a measure of interaction effects associated with an uncertain input parameter. TSI are a more reliable measure of the overall effect of a factor on model output than first order sensitivity indices especially when the interactions between the param eters are considerable (Saltelli et al ., 2000). The number of runs required to implement Sobol for computation of first order and total sensitivity indices is expressed as: (3 11) where N is number of model executions, n is samples size, and k is the number of input factors. Saltelli et al (2005) recommend n = 500 to 1000 to get stable results. For this study, we used n=512 and k depended on the number of important parameters identified from the Morris screening. GLUE Methodology The GLUE methodology was implemented in four steps: steps 1 to 3 were similar to those conducted under section global sensitivity analysis sensitivity indices and in step 4, an evaluation procedure was performed for every single simul ation performed in step 3. Simulations and thus parameter sets were rated according to the degree to which simulated water table elevation matched measured water table elevation. The NSE was used as the likelihood measure and was calculated by the Matlab c ode described under section global sensitivity methods outside of the GLUE analysis. Simulations with NSE close to 1 were accepted while simulations with NSE close to zero were rejected. Using the likelihood measure assigned to all acceptable parameter se ts, a discrete joint likelihood function was generated :

PAGE 108

108 (3 12) w here is the error variance and is the variance of the observed water table elevation NSE takes a value of 1 for a perf ect model fit, a value of less than 0 6 5 would imply that the mean value of the observed data would be a better predictor than the simulation model (Krause et al ., 2005; Stedinger et al ., 2008). Since the discrete joint likelihood function can only be illu strated in maximum 3 dimensions, scatter plots (dotty plots) were used to illustrate estimated parameters. Finally, the likelihoods were projected onto the parameter axis and discrete posterior PDFs and corresponding CDFs were generated for each parameter. We implemented the GLUE methodology using a software package called GLUEWIN (Ratto and Saltelli, 2001). Within the GLUEWIN environment, only the sample file, model output file, and likelihood file were required to estimate parameter posterior distribution s and statistics. Results and Discussions Parameter Screening: Morris Method The ranking of the relative importance of STWT1 input parameters based on Morris method and two sampling techniques are presented in Figures 3 2 and 3 3 with the greater the separ ation from the origin of the versus plane corresponding to greater importance of the parameter. The number of STWT1 model parameters identified as important for predicting water table elevation using NSE as our model output measure were reduced from six to four in both sampling techniques. Overall, Morris screening of parameters based on the random sampling technique and the improved sampling strategy by Campolongo et al (2007) were similar ( Figures 3 2 & 3

PAGE 109

109 higher for t maximum dispersion of trajectories within the model input factor space was able to capture effects of parameters on model output better than random sampling techniques implement ed in the original Morris (1991) method. It is worth noting that the location of all the STWT1 parameters within the below the imaginary 1:1 line within the were primarily first order with minimum parameter interactions. Global sensitivity analysis of the STWT1 model using the Morris method shows a strong influence of specific yield (ASY) at all six wells ( Figures 3 3 and 3 4) where ASY is furthest along the axi axis. As specific yield characterizes the increase in water table elevation or drawdown due to recharge or pumping, respectively, its strong influence on predicting water table elevation is expected. Others have reported similar f indings. Using OAT local sensitivity analysis methods and a two dimensional MODFLOW model in Biscayne aquifer, Bolster et al (2001) also indicated the sensitivity of water table elevations to specific yield especially during extreme fluctuations. Similarl y, Kisekka and Migliaccio (2012) observed that MODFLOW simulated water table elevation in the Biscayne aquifer was most sensitive to specific yield using PEST (Parameter Estimation) local sensitivity analysis. The other model parameters positioned away fro m the origin of the are aquifer saturated thickness, horizontal hydraulic conductivity, and canal leakance. Vertical hydraulic conductivity and half canal width were not important given their location within the Figures 3 2 & 3 3). T hese results confirm that groundwater flow in the Biscayne is primarily horizontal i.e., Dupuit assumptions are

PAGE 110

110 v alid Bolster et al (2001) confirmed Dupuit assumptions for Biscayne aquifer by placing piezometers at different depth at the same location an d observed that there were no distinguishable differences in head measured by piezometers at different depths. The sensitivity analysis approaches used in the previous studies did not provide qualitative relative ranking of all model parameters in the mod els applied to calculate water table elevation as was achieved using the Morris method in this study. The practical application of the Morris results is directly related to management of flood events at the study area in which canal stage is usually lowere d to create storage before a storm event and the associated transient aquifer response is primarily influenced by specific yield. Therefore, properly characterizing Biscayne aquifer specific yield and its interactions with other parameters provides improve d prediction of transient aquifer responses. Global Sensitivity Analysis: Based on Morris ranking of parameters, a subset of STWT1 important parameters (i.e., specific yield [ASY], horizontal hydraulic conductivity [AKX], aquifer thickness [AB], and canal leakance [XAA]) were used in further variance based global parameters at the six groundwater observation wells. Figure 3 4 depicts the fraction of the tota l output variance explained by each of the four parameter using both first order and total order Sobol sensitivity indices (vertical axis). For each parameter, the first order effects are presented first followed by the total order effects and the differen ce between the two represents higher order effects or parameter interactions. Results from Sobol analysis confirm and quantify results from Morris analysis. At the six wells, specific yield explained over 60% of the total variance in predicted water table elevation

PAGE 111

111 (expressed using NSE and RMSE) followed by aquifer thickness explaining approximately 20%. The results also show that the effect of canal leakance reduces as the distances from the canal increases ( Figure 3 4; e.g., wells C 111AW and C 111AE are 500 and 2000 m from canal C 111 respectively). It is worth noting that the Sobol method is more robust than the Morris method since it is based on a large number of model simulations and a less structured sampling method (Saltelli et al ., 2004; Muoz Carpena et al of all the first order parameter effects for the STWT1 model are approximately 100% ( Figure 3 4) indicating that the STWT1 model behaves as an additive model. This implies th at the STWT1 model can be efficiently calibrated if reliable data are available (Muoz Carpena et al ., 2007). Parameter Estimation for STWT1 Using GLUE Posterior distributions for the four parameters at the six observation wells were produced for the 5120 model simulation s and for the top 10% (best performing in terms of fit between predicted and measured water table elevation) model simulations based on the likelihood measure NSE. In the proceeding analysis, only posterior distributions corresponding to th e top 10% of model simulations are presented. For brevity and to facilitate graphical representation, only histograms and CDF plots for well C 111AW and C 111AE are shown in Figures 3 5 and 3 6. In the GLUE analysis, the shape of the posterior distributio ns also indicated the degree of uncertainty of the parameter estimates, sharp and peaked distributions are associated with well identifiable parameters while flat distribution indicate greater parameter uncertainty. Posterior distributions for specific yie ld (ASY) were sharp and peaked indicating less parameter uncertainty ( Figure 3 5). AKX, AB, and XAA had less sharp and wider distributions

PAGE 112

112 compared to ASY indicating more uncertainty compared to parameter ASY. This is also confirmed by the standard deviati on values in Table 3 3. Again, for brevity only dotty plot for well C 111AW ( Figure 3 7) is shown as an example and indicates that variation in NSE was greater for specific yield and to a less extent saturated thickness (AB), AKX and XAA. Model response to variation in the four parameters was similar at all the six sites probably due to the small variation in water table elevation at the different wells. The posterior CDF show mean values as well as 5% and 95% quintiles for the model parameters. GLUE estima ted parameters are shown in Table 3 3 Model Fit From the posterior distributions, parameter estimates were inferred as model values corresponding to the top 10% or best model simulations with the highest likelihood measure which was in our case NSE closes t to 1.0. As an example, Figure 3 8 shows good fit between measured and predicted water table elevation at three of the six wells along the same transect from C 111 canal. Water table elevation was predicted using the STWT1 model and parameters values obt ained from the GLUE analysis. The NSE s at all the wells are shown in Table 3 3 Accuracy of model fit decreased with increasing distance from the canal from NSE of 0.95 at 500 m from the canal to 0.81 at well C 111AW located 2000 m from the canal. The RMSE also increased from 3.5 cm at well C 111AW to 6.7 and 6.5 cm at wells AK6 and C 111AE, respectively. This model behavior would be expected because the model assumes that as distance from the canal becomes large, change in water table elevation approaches the initial water table in the system. The model appears to be more accurate within distances of 2000 m from the canal. The GLUE results show that both the model and likelihood function used in this study were realistic.

PAGE 113

113 Comparison of Estimated Parameters to Literature Values The value of specific yield estimated in this study is within range of values estimated using other methods in literature. In this study we estimated an average specific yield of 0.10. Bolster et al (2001) used a complex canal drawdow n field experiment to estimate specific yield for Biscayne aquifer near Everglades National Park as 0.15, while Muoz Carpena and Li (2003) used high temporal resolution (15 minute interval) groundwater data estimated specific yield near our study site by dividing precipitation based recharge (after a large storm) by change in groundwater head along a transect in Biscayne aquifer to obtain specific yield of 0.11 at several wells. Schroeder et al (1958) reported specific values for Biscayne aquifer ranging between 0.1 0.35. Table 3 3 shows values of horizontal hydraulic conductivity estimated in this study were within range of 7,590 to 14,000 m/day obtained by Genereux and Guardiario (1998) using a large scale canal draw down experiment. Our values were clo sest to values estimated by Fish and Stewart (1991) and Chin (1991) ( 12,187 and 12,500 m/day respectively) that used stepped drawn down pumping tests and transmissivity approach. Fisher and Stewart (1991) explained the difficulty in estimating aquifer par ameters such as hydraulic conductivity for very high yielding aquifers such as Biscayne aquifer, in their pumping test studies aquifer water levels recovered within 1 to 2 minutes which was to quickly for practical measurements of drawdown. Therefore, the simple and cheap method employed in the study provides a globally based method for estimating aquifer parameters for high yielding aquifers. Estimated aquifer thickness at our study site is also within range of that estimated by Boslter et al (2001) west of C 111 and is also close to the approximate values estimated based on Fish and Stewart (1991) total thickness of Biscayne aquifer.

PAGE 114

114 There is little literature on the values of C 111 canal leakance; however, our average value of 189 is close to the value of 148 estimated by Bolster et al (2001) for the same canal. The similarity in the values estimated in this study using GLUE to values estimated using other methods confirms and provides confidence for our results. This work demonstrates the simple but ro bust procedure for characterizing parameters governing surface groundwater interactions. The values estimated in this study can be used for predicting transient aquifer responses using more complex numerical models for flood event management, ecosystem man agement, and even water supply well field management. Conclusion Global sensitivity analysis using Morris screening method (original and improved STWT1 approximate analytica l model of canal aquifer interaction to assess the influence of parameters on the exchange of water between canal C 111 and Biscayne aquifer and to better quantify selected physical parameters. Using the qualitative Morris method, important STWT1 parameter s were ranked. Parameter ranking based on the original random sampling technique and the improved sampling strategy was the same but the magnitudes of sensitivity measures were high for the latter probably indicating better characterization of parameter ef fects. Ranking indicated that only four parameters were important for explaining aquifer response to transient canal stage and important parameter explaining transie nt aquifer responses to stresses. STWT1 was determined to be an additive model, all parameters had primarily first order effects with negligible parameter interactions implying that this model could be accurately calibrated

PAGE 115

115 using reliable measured data. Pa rameter values for the most sensitive parameters were estimated using the GLUE method. Posterior distributions from GLUE indicated sharp and narrow probability distributions for specific yield implying that this parameter could be estimated with minimum un certainty. The values of parameters estimated using GLUE resulted in good model fit especially for wells within 2000 m of the canal with 0.8
PAGE 116

116 Table 3 1 Monitoring groundwater wells with descriptors in the C111 basin Well name Distan ce from C 111 canal m Well installer Latitude Longitude Ground elevation NGVD29 m C 111AW 500 SFWMD1 25.39317 80.553724 1.21 AK6 1000 UF2 25.39283 80.549543 2.23 C 111AE 2000 SFWMD 25.39261 80.541605 1.19 AK5 2000 UF 25.40347 80.541933 2.07 VC2 1000 UF 25.41110 80.550375 1.86 VC1 1000 UF 25.41883 80.550041 2.07 1 South Florida Water Management District 2 University of Florida Table 3 2 Summary of STWT1 model uncertain parameters and their probability distributions Para meter Base value Range value or std 1 Units PDF Reference Horizontal hydraulic conductivity ( ) 1218 7 1844 m/day Lognormal Fish and Stewart (1991) Vertical hydraulic conductivity ( ) 614 78 1587 m/day Uniform Ge nereux and Guardiario (1998) stream bank leakance ( ) 217 83 360 Uniform Genereux and Guardiario (1998) Specific yield ( ) 0.15 0.05 0.57 Triangular Bolster et al. (2001) Saturated thickness( ) 13.6 8 19 m Uniform Bolster et al. (2001) Half width of canal ( ) 15.5 10 21 m Uniform Measured 1 Standard deviation Table 3 3 GLUE estimated parameters for the canal aquifer interaction model STWT1 Well ASY 1 ( ) AKX 2 (m/day) AB 3 (m) XAA 4 (m) Model NSE 5 threshold C 111AW 0.1060.032 12740156 15.832.40 165.665.85 0.95 0.96 AK6 0.1030.029 12760161 15.902.38 187.573.47 0.82 0.86 C 111AE 0.1000.028 12790160 15.972.34 203.176.79 0.81 0.90 AK5 0. 1000.027 12780160 15.962.35 203.476.59 0.66 0.70 VC2 0.1020.029 12770161 15.922.37 187.573.51 0.80 0.82 VC1 0.1020.029 12770161 15.922.36 187.373.43 0.76 0.80 1 Specific yield 2 Horizontal hydraulic conductivity 3 Aquifer thickness 4 Canal leaka nce 5 Nash Sutcliffee coefficient of efficiency

PAGE 117

117 Figure 3 1 Map of the study area showing University of Florida (UF) and South Florida Water Management District (SFWMD) experimental sites and the SFWMD canal network in the lower C111 agricultural basin

PAGE 118

118 Figure 3 2 Morris screening results for STWT1 model applied in the Biscayne aquifer, where ASY is specific yield, AB is aquifer thickness, AKX is horizontal hydraulic conductivity, XAA is canal leakan ce, AKZ is vertical hydraulic conductivity, and XZERO is half canal width

PAGE 119

119 Figure 3 3 Improved Morris sampling screening results for STWT1 model applied in the Biscayne aquifer, where ASY is specific yield, AB is aquifer thick ness, AKX is horizontal hydraulic conductivity, XAA is canal leakance, AKZ is vertical hydraulic conductivity, and XZERO is half canal width

PAGE 120

120 Figure 3 4 Sobol indices for the canal aquifer interaction model STWT1 at 6 groundwa ter observation wells in the Biscayne aquifer, where ASY is specific yield, AB is aquifer thickness, AKX is horizontal hydraulic conductivity, and XAA is canal leakance.

PAGE 121

121 Figure 3 5 Posterior probability density functions (PDF) and cumulative density functions (CDF) for STWT1 model parameters estimated using GLUE, the dots on the cumulative density function represent 5% and 95% quartiles for specific yield (ASY), horizontal hydraulic conductivity (AKX), aquifer thickness (AB) an d canal leakance (XAA).

PAGE 122

122 Figure 3 6 Posterior probability density functions (PDF) and cumulative density functions (CDF) for STWT1 model parameters estimated using GLUE, the dots on the cumulative density function represent 5% and 95% quartiles for specific yield (ASY), horizontal hydraulic conductivity (AKX), aquifer thickness (AB) and canal leakance (XAA).

PAGE 123

123 Figure 3 7 Dotty plots from GLUE analysis showing change in the likelihood measure NSE (Nash Sutcliffe coefficient) over the range of model parameters where the grey dots represent the 5120 model runs at well C 111AW.

PAGE 124

124 Figure 3 8 STWT1 predicted and measured water table elevation time series along a transect from cana l C 111 at three wells located at distances 500, 1000 and 2000 m for wells C 111AW, AK6 and C 111AE respectively.

PAGE 125

125 CHAPTER 4 SIMULATING WATER TABLE RESPONSE TO PROPOSED CHANGES IN SURFACE WATER MANAGEMENT IN THE C 111 AGRICULTURAL BASIN OF SOUTH FLORIDA In troduction The C 111 canal constructed in 1966 is the southernmost canal of the central and south Florida canal syst em and serves a 259 square kilometer basin The primary function of the C 111 canal system is to p rovid e flood protection and drainage for a gricultural areas along the eastern boundary of Everglades National Park (ENP). Past dredging of the C 111 canal redirected water flow causing water to flow east from ENP into C 111 ( Figure 4 1). This resulted in reduced flows in Taylor Slough thereby impa cting water quality, fisheries and ecology of Florida Bay ( U.S. Army Corps of Engineers [USACP], and South Florida Water Management District [SFWMD] 20 11). Taylor Slough is a natural drainage feature of the Everglades that empties its fresh water into Flor ida Bay ( Figure 4 1 ) The re direction of water flows to the east results in approximately 6.4 million cubic meters of water a day to be removed from the Taylor Slough system ( US Army Corps of Engineers 2009). To address some of the unintended consequenc es of the canal system, h ydrological modifications are occurring in south Florida as part of the Comprehensive Everglades Restoration Plan (CERP), whose overall goal is to restore the natural ecosystem that was negatively impacted by an extensive canal net work originally constructed to allow for development and provide flood protection (USGS, 1999). One of the 68 components that make up the CERP is the C 111 spreader canal project (U.S. Army Corps of Engineers [USACP] and South Florida Water Management Dist rict [SFWMD], 20 11 ). Through operational adjustments and structural modifications the goal of the C 111 spreader canal project is to restore the quantity, timing and

PAGE 126

126 distribution of water delivered to Florida Bay via Taylor S lough to levels as near as pos sible to pre drainage conditions while maintaining flood protection for nearby agricultural lands. In addition, there is a goal to restore hydroperiods that support pre drainage vegetation patterns in ENP. To achieve the objectives, o perational adjustments are proposed that include incremental ly rais ing canal stage by 3 cm per year up to a maximum of 12 cm at structure S 18C which is a gated spillway ( Figure 4 1). It is anticipated that the raise in C 111 canal stage will affect water table levels in the a djacent agricultural lands ( Figure 4 1). Earlier research has indicated substantial interaction between the highly permeable Biscayne aquifer and surface water in south Florida canals (Graham et al., 1997; Genereux and Slater, 1999; Lal, 2001; Ritter and M uoz Carpena, 2006). The hydraulic connection between Biscayne aquifer and C 111 causes the shallow water table system to fluctuate with respect to changes in canal stage. An increase in water table elevation, due to raised canal stage could result in prol onged root zone saturation or temporary groundwater flooding ( groundwater flooding occurs in low lying areas when the water table rises above the land surface [USGS, 2000] ) which could affect agricultural production in lands adjacent to ENP. Prolonged satu ration of the root zone or short term groundwater flooding could impact yield potential through impaired root growth caused by anoxia, reduced stomatal conductance and net CO 2 assimilation (Schaffer, 1998). It is not known how the proposed operational adju stments (involving incremental raises in canal stage) along C 111 canal would impact water table elevation which would in turn impact optimum crop growth in adjacent farmlands.

PAGE 127

127 MODFLOW, a widely used numerical groundwater flow computer code from the Unite d States Geological Survey (USGS), has previously been used in investigations of canal aquifer interactions in south Florida ( Wilsnack et al., 2000; Bolster et al., 2001 ; Saier s et al. 2004; Hughes et al., 2012) In MODFLOW modeling, various approaches ex ist for representing a surface water body either as a head dependent boundary flux using the river package ( McDonald and Harbaugh, 1988 ) or by using more complex approaches that implicitly couple a numerical open channel flow model to MODFLOW such as MODBR ANCH developed by Swain and Wexler (1996). Although MODFLOW based groundwater flow models have been used to simulate Biscayne aquifer in south Florida ( Hughes et al., 2012), most are regional and they lack the spatial resolution to address water resources issues at field scale particularly groundwater flooding issues in agricultural fields that are influenced by small scale micro topography. For example Brion et al. (2001) used the South Florida Regional Simulation Model in the south Florida E verglades with a grid size of 3.2 km x 3.2 km. The study purpose was to investigate through monitoring and modeling the effect of the proposed incremental raises in C 111 canal stage on water table elevation levels in agricultural areas adjacent to ENP. The objectives w ere to: (1) develop a MODFLOW based model for simulating groundwater flow within the study area, (2) apply the developed model to determine if the proposed changes in canal stage result in significant changes in water table elevation root zone saturation or flooding (i.e., groundwater flooding) and (3) assess aquifer response to large rainfall events and explore the effect of pre storm canal stage drawdown in the mitigation of groundwater flooding of agricultural lands

PAGE 128

128 Materials and Methods S tudy Area The agricultural study site included approximately 17 km 2 located within the C 111 basin in southern Miami Dade County, near Homestead, Florida, United States ( Figure 4 1). The hydrogeologic system at the study site is described in detail in Kisekka et al. ( 2 013a ) and consists of the Biscayne aquifer bordered by two canals C 111 and C 111E approximately 3.5 km apart ( Figure 4 1). The canals are managed by the SFWMD. Aquifer thickness at our study site was not measured but Genereux and Guardiario (1998) report ed a thickness of 13.6 m for a site west of our current study site with roughly one third accounted for by the Miami Limestone formation. Kisekka et al. ( 2013b ) applied inverse modeling using a quasi canal aquifer interaction model and estimated Biscayne a quifer thickness at our study site to range between 13.5 and 18.2 m with an average of approximately 16 m. Canal C 111 was constructed in 196 6 as the principle flood control canal for south Miami Dade County and partially penetrates the Biscayne aquifer t o a depth of approximately 5 m (i.e., 4 m through the Miami Limestone formation and 1 m into the Fort Thompson Limestone formation). Flow in C 111 is south towards Florida Bay and topography is essentially flat ranging between 1.0 to 2.2 m National Geodeti c Vertical Datum (NGVD) 29. The width of the canal increases towards the south with an average of approximately 29 m at the S177 spillway ( Figure 4 1). Currently little is known about hydraulic properties of canal bed sediment in the lower C 111; however, presence of a low permeability canal bed sediment layer which is a mixture of carbonate mud and natural organic matter in several canals within the C 111 basin ha s been documented (Chin, 1991; Genereux and Guardiario, 1998; Merkel, 2000). Using inverse mod eling

PAGE 129

129 and a quasi canal aquifer interaction model as demonstrated Kisekka et al. ( 2013b ) estimated the ratio of canal bed thickness to bed sediment hydraulic conductivity as 0.015 ( ranging between 0.009 and 0.020) days which is close to the 0.029 days esti mated by Bolster et al. (2001) for nearby canal L 31W ( Figure 4 1). Biscayne aquifer hydraulic conductivity was also estimated by Kisekka et al. ( 2013b ) as 12,768 m/day which is within range of values estimated by other s at nearby sites (Fish and Stewart, 1991; Genereux and Guardiario, 1998). Specific yield at our study site was estimated as 0.102 ( ranging between 0.07 and 0.13) by Kisekka et al. ( 2013b ) which is close the value of 0.15 estimated using a large scale canal drawdown by Bolster et al. (2001). Canal aquifer interaction hydraulic parameters will be determined using inversing modeling in the present study. Hydrologic Time Series Data from six groundwater observations wells were used ( Table 4 1 ). Data were collected from August 2010 to February 201 3 Observation wells 4 (C111AE) and 6 (C111AW) were maintained by the SFWMD while the other s ( 1 [VC1] 2 [VC2] 3 [AK5] and 5 [AK6] ) were maintained by University of Florida (UF) IFAS (Kisekka et al., 2013b ) UF wells were equipped with level loggers (Le velogger, Gold Solinst Canada Ltd., 35 Todd Rd, Georgetown, Ontario, Canada) to record water table elevation every 15 minutes although daily averages were used in modeling. Atmospheric corrections were accounted for using a STS Barologger (Solinst Canada Ltd) in well 5 ( Figure 4 1). Data were downloaded weekly and as a quality control procedure, water table elevations were also measured manually with a Model 102 Laser water level well meter (Solinst, Canada Ltd). Elevations at the top of the well manholes were measured using a laser level with reference to a SFWMD bench mark with elevation 1. 19 m NGVD29 near well

PAGE 130

130 4 Water table elevation data for wells 4 ( C 111AE ) and 6 ( C 111AW ) are processed by SFWMD and published on DBHydro ( http://www.sfwmd.gov/dbhydroplsql/show_dbkey_info.main_menu ). Canal stage data were measured along the reaches of C 111 and C 1111E surrounding the study area by the SFWMD Canal stage data are publically avai lable at the online environmental database called DBHydro ( http://www.sfwmd.gov/dbhydroplsql/show_dbkey_info. main_menu ). Groundwater Flow Simulation MODFLOW a widely used modul ar three dimensional finite difference groundwater flow simulator developed by the USGS was selected to simulate groundwater flow in the agricultural lands adjacent to C 111 canal. The governing equation for saturated flow in porous media implemented in M ODFLOW is (McDonald and Harbaugh, 1988; Harbaugh et al. 2000): ( 4 1) where h [L] is the hydraulic head or water table elevation, S s [L 1 ] is the specific storage of the porous media, K xx K yy and K zz [L T 1 ] are hydraulic conductivi ty along the x, y, and z directions, t is time [T], W [T l ] is a source/sink term, with W > 0 for flow into the aquifer and W < 0 for flows out of the aquifer. Due to its computational efficiency and the improved ability to control the conversion between w et and dry cells, the Preconditioned Conjugate Gradient (PCG) package was used to solve the finite difference equations at each time step of the MODFLOW stress period. The hydrogeologic system was modeled as a one layer unconfined aquifer with 2D horizonta l flow. For cases such as this, MODFLOW modifies Eq uation 4 1 by

PAGE 131

131 substituting the specific storage with the specific yield and allows transmissivity to vary based on the changes in aquifer saturated thickness. The assumption of predominately horizontal flo w was based on earlier investigations by Genereux and Guardiario (1998) that showed generally zero difference between piezometers installed at various depths into the Biscayne aquifer. Boundary conditions The following boundary conditions were used in the simulation: canals stage, evapotranspiration, and recharge. The bottom boundary was described as a no flow boundary consistent with observed 2D horizontal flow in the study area. Canals C 111 and C 111E formed the west, east, and south boundaries of the f low domain ( Figure 4 1). C 111 is the larger of the two canals with an average width of 29 m near the gated spillway at structure S177. Both canals partially penetrate the Biscayne aquifer with C 111 having an average depth of approximately 5 m. Water leve ls in C 111E are controlled using a gated culvert at structure S178 ( Figure 4 1). C 111E joins C 111 at the southern tip of the flow domain to become one canal. Surface water groundwater interactions were simulated using the RIVER (RIV) package. The RIV package was selected as a simple and adequate representation of the interaction between the C 111 canals and Biscayne a quifer Canal stage data for reaches of C 111 and C 111E surrounding the study area were obtained from DBHydro. In the RIV package both canal stage and canal conductance (Eq uation 4 2) control the extent of water exchange between the aquifer and the canals ( 4 2)

PAGE 132

132 where C is canal conductance [L 2 T 1 ], K s is the hydraulic conductivity of the low permeability bed sedim ent [LT 1 ], W is the width of the canal [L], L is the length of the canal reach [L], and d is the thickness of the sediment layer [L]. The canal conductance multiplier in MODFLOW was set to range between 702 and 1560 m 2 /day for headwater and tail water rea ches of C 111 based on estimates of the K s /d ratio by Kisekka et al. ( 201 3b ). Given the substantially smaller size of C 111E compared to C 111, canal conductance multiplier for C 111E was set to values ranging from 200 to 500 m/d with lower values assigned to the headwater side of the S178 gated culvert. Given the relatively flat topography, the average of tail water canal stage at S177T and the headwater stage at S18C were used to represent the west and south boundary conditions for all cells downstream of S177 while head water canal stage at S177H was used to represent canal stage for all cells north of S177. Similarly canal stage at S178 (tail waters) and S18C (headwaters) represented the east boundary condition for all cells. Canal stage data were meas ured by the SFWMD and are publically available on DBHydro. The northern boundary was described as a general head boundary using groundwater levels from a monitoring well installed by University of Florida (i.e., well 1 ). Evapotranspiration was simulated us ing the EVT package in MODFLOW (McDonald and Harbaugh, 1988) in which the elevation of the evapotranspiration surface was set to 1.0 m and the evapotranspiration extinction depth to 0.9 m based on ranges reported in Chin (2008) and water table elevation re corded during the study period. We assumed that for water table depths less than 1 m from the land surface, evapotranspiration occurred at the potential rate which was computed from micro meteorological data obtained from a Florida Automated Weather Networ k (FAWN;

PAGE 133

133 http://fawn.ifas.ufl.edu/) station located 15 km northeast of the study site. The ASCE standardized Penman Monteith equation and the REF ET tool by Allen (2011) were used to estimate ET o values. Recharge to the aquifer was simulated using the RCH package. To minimize the uncertainty associated with spatial variability of rainfall in south Florida, gauge adjusted NEXRAD (Next Generation Radar) rainfall data were used (Skinner et al., 2008). Space and time discretization The finite difference grid used with the MODFLOW simulations consisted of a single layer covering approximately 17 km 2 The model layer was discretized into 69 rows (running east to west) and 46 columns. Nodal spacing for the columns ranged from 53.5 m to 105.6 m from west to east w ith the smallest spacing closest to the canal since this is where greater changes in hydraulic head would be expected. Nodal spacing for the rows was constant over the model domain and set to 100.6 m. Further reductions in discretization did not appear to improve simulations results. All the spatial discretization was implemented using a pre and post processor for MODFLOW called MODFLOW GUI PIE version 4.34.00 an Argus One Plug In Extension (PIE) (Winston, 2000). The model simulated conditions from 25 Augu st 2010 to 28 February 201 3 The time step and stress period sizes were set to one day; the multiplier was also set to one day. The period from 25 August 2010 to December 2011 was used to calibrate the model, while the data from 01 January 2012 to 28 Febru ary 201 3 was used to validate the model. It was assumed that canal stage did not change during each stress period which was reasonable because 24 hour variation s in canal stage are small unless a large rain event occurred or an operational change in canal stage management was

PAGE 134

134 implemented. Initial conditions over the model domain were obtained from observation well data at the start of the simulation and interpolated over the model domain using Argus ONE interpolation utilities. Sensitivity analysis and para meter estimation Sensitivity analysis and parameter estimation were performed using the sensitivity and parameter estimation (PES) processes in MODFLOW 2000 (Hill 1998; Hill et al., 2000). PES calculated parameter values that minimized a weighted least sq uares objective function using nonlinear regression. The objective function was minimized using the modified Gauss Newton (also known as the Levenberg Marquardt method) as well as prior information on the parameter estimates (Hill et al., 1998). To reduce problems associated with inverse modeling such as insensitivity, instability and non uniqueness, only parameters identified through sensitivity analysis to have greatest influence on model output were estimated. The sensitivity equation method was used in the sensitivity analysis package. Output from MODFLOW 2000 also includes inferential statistics such as dimensionless scaled sensitivities (DSS) and composite scaled sensitivities (CSS). These inferential statistics measure the amount of information provi ded by the observations and the uncertainty with which the parameters values are estimated (Hill, 1998). DSS are typically used to compare the importance of different observations for estimation of a single parameter. CSS are calculated for each parameter using DSS for all the observations and indicate the amount of information provided by the observations for the estimation of a single parameter.

PAGE 135

135 Model validation Model validation was implemented using a statistical model evaluation tool called FITEVAL (Rit ter and Muoz Carpena, 2012). FEITEVAL computes a non dimensional goodness of fit indicator C eff (Nash Sutcliffe coefficient of efficiency), a dimensional goodness of fit indicator RMSE (Root Mean Square Error) as well as model prediction uncertainty range s. FITEVAL computes a 95% confidence interval based on a goodness of fit probability density function estimated using bootstrap technique. FITEVAL also provides some reference values as guides for judging model performance. The model is judged to be very g ood if C eff > 0.9, good if C eff is between 0.8 and 0.9, acceptable for C eff between 0.65 and 0.8 and unacceptable for C eff < 0.65. Model A pplication: Canal stage O perational A djustment S cenarios Before application of the model, graphical exploration of the temporal variation in water table elevation in reference to the root zone was completed to determine if under present canal stage operational criteria water table elevation extended into the root zone The developed model was then applied to evaluate the effect of the proposed incremental raises in canal stage on water table elevation. Incremental raises in canal increments of 3 cm up a maximum of 12 cm ( U.S. Army Cor ps of Engineers and SFWMD, 2011). For numerical simulation purposes, i ncremental raises in canal stage were mimicked by adding the proposed increments of 6, 9 and 12 c m to measured canal stage. Only tail water canal stage at S177 and S178 were modified. Ca nal stage of the head waters at S177 and S178, rainfall, and evapotranspiration from the period of record were used. The initial condition was taken as the interpolated surface for water table elevation at the start of the simulation. Graphical analysis wa s used to determine if

PAGE 136

136 the proposed increments in canal stage would result in root zone saturation and flooding at any of the sites analyzed. The Two sample equal variance t Test was used to determine if the water table elevation before and after the incre mental rises in canal stage were significant. Assessing Aquifer Response to Large Storms When a large storm is forecasted, t he SFWMD uses data products from the National Hurricane Center (NHC) to make pre and post storm operational plans. These include mak ing forecasts of quantitative precipitation that are accurate with in 2 4 days prior to the storm and corresponding regional canal level lowering to ensure continued flood protection. During the storm even t, the SFWMD continues to monitor flood control stru cture s as well as storm position and intensity. During T ropical S torm Isaac the SFWMD requested the USACE to put C 111 in pre storm mode in order to minimize potential impacts USACE approved pre storm drawdown request and gate openings and pumping were in itiated Aug ust 23 2012 (Strowd, 2012) T he period August 21 to 30, 2012 was chosen for the analysis of Biscayne a quifer response to large storms as this period corresponded to Tropical Storm Isaac (> 60 mm total rainfall in one day). To simulate aquifer response to large storm events, MODFLOW was used with a small time step of 15 minutes. A stress period size of one day was also used to match available tail water canal stage and precipitation data at S177 and S178 ( Figure 4 1). Model simulations of aquife r response to recharge were checked using : ( 4 3)

PAGE 137

137 where S y is the aquifer specific yield, recharge refers to input from rainfall and head difference refers to the change in water table elevation resulting from the recharge. The MODF LOW model was also applied to assess aquifer response to two, five, ten and 25 year return period storms. Maximum daily rainfall depth for the return periods were obtained from isohyetal maps for central and south Florida developed by Pathak (2001). Pathak (2001) obtained maximum daily depth of 114, 168, 203, and 254 mm for two five ten and 25 year return period storms respectively for our study area. Based on analysis of over 113 years of rainfall data, large storm s (i.e., 2 to 25 year return storms ) in south Florida occur between August and October. With October being a transitional month between the wet and dry seasons and also corresponds to the time when growers begin to prepare the land and plant winter vegetables. For this reason, the period from October 25 to November 5, 2012 was selected to explore canal aquifer system responses to large storms. V arious canal drawdown scenarios that would minimize flooding in the agricultural lands were also explored. Drawdowns were implemented by incrementally reducing canal stage 48 hours prior to a forecasted large storm in the reaches of C 111 and C 111E surrounding the study area. The desired scenario was when the water table elevation did not exceed the elevation of the bottom of the root zone. Results and Discussion Sensitivity Analysis a nd Parameter Estimation Results The CSS for our study using MODFLOW summarized in ( Figure 4 2) indicated that water table elevation measurements provided more information in the estimation of specific yield and hydraulic co nductivity. The CSS also indicate that water table elevation data alone did not provide sufficient information for accurate estimation of

PAGE 138

138 canal bed conductance in the reaches of C 111 and C 111E surrounding our study site. The need to have different types of data during parameterization of groundwater flow models was noted by earlier investigators ( Saier s et al., 2004 ; Zechner and Frielingsdorf, 2004). Zechner and Frielingsdorf (2004) observed that to accurately parameterize a canal aquifer interaction mode l with many parameters, in addition to groundwater head observations other observations such as canal seepage and pore water solute concentration provided more information for parameter estimation and improved model prediction. However Saier et al. (2004) using different combinations of observed data including groundwater head, aquifer discharge to the canal and groundwater chloride concentration noted that inverse solution uniqueness was not required for accurate prediction of groundwater head but was re quired for prediction of seepage. Their results i mply that water table head observations are sufficient for calibrating models whose goal is prediction of groundwater head but models for predicting other state variables such as seepage should be calibrated with more than one type of observation. During parameter estimation, the least squares objective function was minimized after five iterations Based on data from 5 observation wells a hydraulic conductivity value of 12,115 m/day was estimated. This value was within the range of 7,590 to 14,900 m/day observed by Genereux et al. (1998) based on a large scale canal draw down experiment and close to 12,768 m/day estimated by Kisekka et al. (2013b). Specific yield was estimated as 0.184 which is close to an est imate of 0.15 determined by Bolster et al. (2001) using data from a large scale canal draw down experiment and to a mean of 0.102 estimated by Kisekka et al. ( 201 3b ). Information from observations

PAGE 139

139 was not sufficient to accurately estimate canal conductanc e along the reach of C 111 on the headwater side of S177. Canal bed conductance multiplier for the longest and largest reach i.e., the reach between S177 and the point where C 111 joins C 111E to become a single canal was estimated as 1,965H m 2 /day, where H represents the length of the reach in meters. Canal bed conductance multiplier for C 111E was less than that of C 111 (i.e., 27H m 2 /day tail water side of S178 and 10H m 2 /day head water side of S178). There are no readily available values for canal bed c onductance for the reaches of C 111 and C 111E considered in this investigation, however, for purposes of comparison, Genereux and Guardiario (1998) estimated canal bed conductance of 720H for the nearby L 31W canal which is located near C 111 along the ea stern boundary of ENP. Model Calibration and Validation Calibration (August 25, 2010 to December 31, 2011) results reproduce d observed WTEs at five observation wells (2 to 6) with an average C eff greater than 0.9 (Table 4 2 ). Temporal variations in WTE s howed seasonal increases and decreases in WTE. The dry season was characterized by decrease in WTE due to low rainfall while increase WTE in the wet season was due to increase in rainfall and changes in canal stage management. The RMSE ranged from 1.0 cm t o 7.0 cm with the lowest value observed at well 6 and the highest value at well 4 Study site topography is essentially flat implying that small variations in hydraulic head govern which direction water flows, therefore it was desired to achieve the lowest RMSE possible (e.g., < 6 cm). However, it was not possible to obtain RMSE < 6 cm at all wells due to limitations e.g., uncertainties in model parameters, model structure, and model input, all of which introduce uncertainties in model simulations. There co uld also be errors in the observed data

PAGE 140

140 used for model calibration. This type of error was minimized by checking level logger data with manual measurements during each download. Under the C 111 spreader canal project the smallest incremental raise in canal stage at S18C is 3 cm. However, the RMSE of our model predictions are larger than 3 cm at four out of the five observation wells within the study area domain, for this reason only the 6, 9 and 12 cm incremental raises in canal stage were further analyzed for their effect on water table elevation. FITEVAL summary of the goodness of fit statistics for validation of model predictions at all the observation wells are shown in Figures 4 3 to 4 7 Overall the agreement between simulated and observed water tabl e elevation s was very good (Ceff > 0.9 and 1 cm < RMSE < 5 cm) with the exception of site well 4 at which model performance was determined to be acceptable (0.65 < Ceff < 0.8). The over prediction at observation well 4 could be attributed to several facto rs e.g., heterogeneity in hydrogeological conditions and uncertainty in model input parameters and observed data. The very good performance of the model at all the other sites indicates boundary conditions were sufficient to describe groundwater flow. The results also indicate that describing canal aquifer interactions using the simple RIV package (Harbaugh et al., 2000) in MODFLOW was adequate for our study site The good performance of the RIV package could be attributed to the underlying assumptions in t he RIV package being valid for our study site e.g., there was negligible change in canal stage during each stress period which was set as one day. Our result s are also within range of model coefficient of efficiency (a measure of agreement between measured and predicted values) obtained by prior investigators Bolster et al. (2001) applied MODFLOW to

PAGE 141

141 Biscayne a quifer and obtained a model coefficient efficiency (calculated in a similar as C eff ) of 0.99. Saiers et al. (2004) using there numerical model of gr oundwater flow and solute transport in the Biscayne a quifer obtained goodness of fit model coefficient efficiency of 0.95. The results also indicated that general groundwater flow was in the south east direction, which implies that a large increase in hydr aulic head west of C 111 could increase rate of groundwater flows to the eastern side of the canal. Based on the period evaluated, model validation results indicated that with the exception of well 4 the model developed for the study area was accurate and not biased implying it could be used to further investigate the impact of proposed incremental raises in canal stage on water table elevation Model Application Results Visual exploration of temporal variation in water table elevation in reference to the root zone ( Figures 4 8 to 4 9 ) revealed that under current canal stage management criteria for the period August 25, 2010 to February 28 2012, water table elevation occasionally extended into the root zone at well 6 and well 4 study sites which also had t he lowest ground surface elevation. The root zone for all sites was approximately the first 20 cm from the ground surface. At well 5 and well 3 sites, where land surface elevation exceeded 2 m, water table elevation was not observed to enter the root zone. Thus, surface topography might influence water table fluctuations into the root zone more than distance from the canal. Results from model applications ( Figures 4 10 to 4 12 ) revealed that the water table elevation before and after 6, 9, and 12 cm increme ntal rises in canal stage were predicted to be significantly different (p<0.001) for monitoring well 4 well 6 and well 5 sites For well 2 and well 3 sites, water table elevation before and after the proposed incremental raises in canal stage were predic ted to not be

PAGE 142

142 significantly different (p>0.05) except for the 12 cm increase at well 3 (Table 4 3) The lack of significant difference in water table levels before and after the incremental raise in canal stage for wells 2 and 3 could be attributed to tha t fact canal stage was not changed north of S177 and S178 (Figure 4 1). The model predicted increase in water table elevation for wells 4, 5 and 6 corresponding to a 6 cm rise in canal stage to range between 4.5 and 6.0 cm, while the model predicted incr eases corresponding to 9 and 12 cm canal increases were 7.0 to 9.0 cm and 11.0 to 12.0 cm, respectively. The almost equal increase in water table elevation predicted from the incremental rises in canal stage can be attributed to the high hydraulic connecti on between Biscayne a quifer and the C 111 canal network. Visual analysis ( Figures 4 1 0 to 4 1 2 ) show s that low elevation lands (as found at well 4 and well 6 sites ) we re predicted to have a shorter growing season with canal stage increases beyond 6 cm resu lting in longer periods of saturated conditions in the root zone. For example, at well 4 and well 6 sites after a 12 cm raise in canal stage, saturated conditions we re predicted to persist until late October or early November. Typically land preparation s tarts in late September and planting in October. For high elevation sites such as well 5 the proposed increases in WTE were predicted not to cause root zone saturation or flooding under conditions similar to those experienced during the study period. Resu lts of Aquifer Response to Large Storms Event analysis was conducted for the period from August 21, 2012 to August 30, 2012 which corresponded to Tropical S torm Isaac. The aquifer responded to the storm with increasing water table elevation A pproximately two days past before the water table elevation recede d back to pre storm levels ( Figure 4 1 3 ) From Equation 4 3, the

PAGE 143

143 three days of heavy rainfall during Tropical Storm Isaac would have theoretically result ed in a change in water table elevation of at le ast 0.6 m or water table elevation of at least 1.4 m NGVD29. However, the simulated water table elevation only reached a maximum of 1.2 m NGVD29 ( Figure 4 1 3 ). As indicated under model validation, performance was ranked as very good at well 6 and well 5 sites and acceptable at well 4 implying the model adequately represented the physical process es in the system. The deviation between the predicted water table elevation and what would theoretically be expected could be attributed to the pre and post T ropic al S torm Isaac canal drawdown that was undertaken by the SFWMD and USACE. This included regional lowering of canal stage particularly by operating canals C 111 under pre storm mode and opening the flood control structures (Strowd, 2012). T ropical S torm Isa ac occurred when the fields at well 5 well 4 and well 6 were fallow so no crop damage occurred. If a similar event were to occur when vegetable crops were present and with no pre storm canal drawdown, sites with lower elevations (e.g., well 4 and well 6 ) would likely experience yield loss due to root zone saturation The event would also delay entry into the field by any machinery for agricultural activity. Higher elevation sites (e.g., well 5 ) were expected to be less impacted by such a storm as the wate r table would still be below the root zone. This further illustrates the need for detailed topographic data and field scale simulation of canal aquifer system to better relate locations with potential risk of groundwater flooding. This model could be used to further explore drawdown scenarios for this area prior to major storm events. Analysis of canal aquifer system response and exploration of various canal drawdowns scenarios that would minimize the impact of root zone saturation and

PAGE 144

144 groundwater flooding in agricultural lands due to large storms revealed that micro topography within the fields was a major factor. Figures 4 1 4 and 4 1 5 show that 3 out of the 4 sites analyzed for their response to two five ten and 25 year return period storms experienced g roundwater flooding if canal drawdown was not implemented before the storm. With the exception of well 5 site with high surface elevation (2.23 m NGVD29), all the other sites experienced various degrees of groundwater flooding ( Figures 4 1 4 & 4 1 5 ). A ten and 25 return period storm caused groundwater flooding at well 2 well 3 and well 6 sites Sites with ground surface elevation less than 1.2 m NGVD29 experienced groundwater flooding from all storm sizes analyzed. For agricultural purposes it is desired th at the water table elevation does not extend into the root zone since this condition could create anoxic conditions that result in root and / or plant death. Exploration of canal drawdown scenarios revealed that a 20 cm drawdown in canal stage 48 hours pri or to a forecasted storm of 114 mm in 24 hours (2 year return period storm) would eliminate the risk of groundwater flooding at all the sites as shown in Figures 4 1 4 and 4 1 5 A 25 cm drawdown was effective in mitigating the impacting of root zone saturat ion from a 5 year return period storm at all sites, while drawdowns of 30 and 40 cm were effective for 10 and 25 year return period storms respectively. It is worth noting that the influence of the 48 hour canal stage drawdown prior to a forecasted storm was dependent on the distance from the canal. As shown by the depressions in the drawdown graphs ( Figures 4 1 4 and 4 1 5 ) for well 6 well 5 and well 3 sites which are 500, 1000, and 1000 m from C 111 canal respectively. Overall these results predict tha t canal drawdown is effective as a pre storm technique for ensuring continued flood protection of agricultural lands within the C 111 basin. However, it is

PAGE 145

145 critical to remember that management decisions should be made in view of the uncertainty associated with model predictions as shown in Figures 4 3 to 4 7 Also, the size of the drawdown should match forecasted storm depth and post storm activities should ensure the canal drainage continues to provide other services such as control of salt water intrusion Conclusion The effect of the proposed incremental raises in canal stage on water table levels along a section of a major canal draining south Florida (i.e., C 111) and aquifer response to large storms was investigated using MODFLOW and graphical analysis The incremental raises in C 111 canal stage are part of a large scale ecosystem restoration project whose goal is to restore the hydrology of ENP The MODFLOW model predicted that the incremental raises in canal stage resulted in significant differences in water table elevation in the agricultural lands south of the spillway at S177. Water table elevation s were predicted to not be significantly different for the lands north of S177. For the 9 and 12 cm increases in canal stage, water table elevations were predicted to occasionally extend into the root zone for 3 out of the 5 well sites graphically analyzed. Well 3 and well 5 sites (with ground surface elevation exceeding 2 m) were predicted to not be affected by any of the incremental raises in stage excep t for 12 cm increase at well 3 The impact of operational changes in canal stage management on the root zone saturation and groundwater flooding depended on land surface topography and depth of rainfall events. Thus micro topography within the field can ha ve a greater influence on soil water content than distance from the canal. Based on graphical analysis, low elevation lands (with surface elevation<2 m NGVD29) could

PAGE 146

146 have shorter growing seasons if canal stage is increased beyond 6 cm due to potential satu ration of the root zone. The MODFLOW model was able to mimic the rise and fall of the water table similar to that measured for Tropical Storm Isaac F urther explor ation of canal aquifer system response to 2, 5, 10 and 25 year return period storms and canal drawdowns suggested that if crops are present during storms greater than a 2 year return period storm, yield losses could occur if pre storm canal drawdown is not implemented at least 48 hour prior to the forecasted storm particularly in low elevation sit es. Overall the study concludes that canal drawdown is effective as a pre storm technique for ensuring continued flood protection of agricultural lands within the C 111 basin.

PAGE 147

147 Table 4 1. Water table elevation monitoring sites with descriptors 1 Site name Distance from canal C 111 (m) Ground surface elevation (m) NGVD29 Latitude Longitude Well 1 1000 2.07 25.41883 80.550041 Well 2 1000 1.86 25.41110 80.550375 Well 3 2000 2.07 25.40347 80.541933 Well 4 2000 1.19 25.39261 80.541605 Well 5 1000 2.2 3 25.39317 80.553724 Well 6 500 1.21 25.39283 80.549543 Table 4 2 Goodness of fit statistics for model calibration for water table elevation predictions using MODFLOW Well Ceff 1 Calibration RMSE 2 Calibration (cm) Well 2 0.97 0.98 4.0 5.0 Well 3 0.9 4 0.96 4.7 5.7 Well 4 0.86 0.91 6.0 7.0 Well 5 0.93 0.95 4.6 5.3 Well 6 0.99 1.00 1.0 1.2 Nash Sutcliffe coefficient of efficiency Root mean square error

PAGE 148

148 Table 4 3 Results for comparison of water table elevation before and after a 6, 9 and 12 cm incr emental raise in C 111 canal stage obtained from a two sample assuming equal variance t February 201 3 ). 6 cm increase Well Mean WTE 1 Before Mean WTE After Variance 2 p value C111AW 0.76 0.8 1 0.01 <0.001 A K6 0.75 0.8 2 0.0 1 <0.001 C111AE AK5 VC2 0.75 0.9 1 1.00 0.8 4 0.9 1 1.00 0.02 0.0 2 0.0 2 <0.001 0. 88 0. 37 9 cm increase C111AW 0.76 0.83 0.01 <0.001 AK6 0.75 0.8 5 0.0 1 <0.001 C111AE AK5 VC2 0.75 0.9 1 1.00 0.8 6 0.93 1.0 0 0.02 0.0 2 0.0 2 <0.001 0. 14 0. 94 12 cm increase C111AW 0.76 0.8 6 0.01 <0.001 AK6 0.75 0.8 8 0.0 1 <0.001 C111AE AK5 VC2 0.75 0.9 1 1.00 0.8 9 0.9 4 1.0 1 0.02 0.0 2 0.0 2 <0.001 0.0 02 0. 30 1 Mean water table elevation in meters National Geodetic Vertical Datum of 1929 2 Pooled variance

PAGE 149

149 F igure 4 1 Study area showing groundwater monitoring sites, agricultural lands adjacent to Everglades National Park (ENP) and canal network within the C 111 basin of south Miami Dade County, Florida

PAGE 150

150 Figure 4 2. Composite scaled sensitivities for the p arameters selected for estimation in the model were Sy is specific yield, H is hydraulic conductivity, C1 is canal bed conductance multiplier for the reach of C 111 on the head water side at S177, C2 is canal bed conductance multiplier for the C 111 reach between S177 and the point where C 111 joins C 111E to become a single canal, C3 is canal bed conductance multiplier for reach of C 111E on the tile water side of S178 and C4 is canal bed conductance multiplier for the reach of C 111E on the headwater si de of S178.

PAGE 151

151 Figure 4 3 Validation goodness of fit indicators from FITEVAL for MODFLOW simulations at well 2 for the period January 01, 2012 to February 28 201 3 Figure 4 4 Validation goodness of fit indicators from FITEVAL for MODFLOW simulatio ns at well 3 for the period January 01, 2012 to February 28 201 3

PAGE 152

152 Figure 4 5 Validation goodness of fit indicators from FITEVAL for MODFLOW simulations at well 4 for the period January 01, 2012 to February 28 201 3 Figure 4 6 Validation goodness of fit indicators from FITEVAL for MODFLOW simulations at well 5 for the period January 01, 2012 to February 28 201 3

PAGE 153

153 Figure 4 7 Validation goodness of fit indicators from FITEVAL for MODFLOW simulations at well 6 for the period January 01, 2012 to F ebruary 28, 2013

PAGE 154

154 Figure 4 8 Temporal variation in water table elevation in reference to ground surface under current canal stage operation criteria at spillway S18C for observation well 2 (ground surface elevation of 1.86 m NGVD29) and well 3 (ground surface elevation of 2.07 m NGVD29) on the headwater side of the spillway at S177 with calibration and validation separated by a vertical dash line.

PAGE 155

155 Figure 4 9 Temporal variation in water table elevation in reference to ground surface elevation under current canal stage operation criteria at S18C for observation wells well 4 well 5 and well 6 on the tail water side of the spillway at S177 with calibration and validation separated by a dash vertical line.

PAGE 156

156 Figure 4 1 0 Temporal variation in water table elevation in reference to the root zone under proposed incremental raises in canal stage operation at S18C for observation well 4

PAGE 157

157 Figure 4 1 1 Temporal variation in water table elevation in reference to the root zone under proposed incremental r aises in canal stage operation at S18C for observation well 5.

PAGE 158

158 Figure 4 1 2 Temporal variation in water table elevation in reference to the root zone under proposed incremental raises in canal stage operation at S18C for observation well 6

PAGE 159

159 Fig ure 4 1 3 Aquifer response to Tropical Storm Isaac at observation wells south of the spillway at S177.

PAGE 160

160 Figure 4 1 4 Canal aquifer system response to large storms of various sizes for wells north of the spillway at S177, were YR refers to year and RP re fers to return period, graphs also shows that canal stage drawdown prior to the forecasted storm reduces the risk of root zone saturation and groundwater flooding

PAGE 161

161 Figure 4 1 5 Canal aquifer system response to large storms of various sizes for wells sou th of the spillway at S177, were YR refers to year and RP refers to return period, graphs also shows that canal stage drawdown prior to the forecasted storm reduces the risk of root zone saturation and groundwater flooding

PAGE 162

162 CHAPTER 5 MODELING SOIL WATER RESPONSE TO SURFACE WATER MAN A GEMENT IN C 111 BASIN USING WAVE CONSIDERING MEASUREMENT UNCERTAINITY Introduction The C 111 basin drained by the extensive C 111 canal is located in Miami Dade County south Florida and covers an area of approximate ly 259 km 2 A substantial proportion of the C 111 basin is characterized by vegetable production primarily during the dry season (October to April). Vegetable production in south Florida is a significant contributor to both the local and state economies. I n 2011 vegetable produc tion in Florida was reported on 185,000 acres with a farm gate value worth of over 1.5 billion dollars ( DACS, 2012 ). Major crops include green beans, sweet corn, squash, tomatoes and sweet potato. It is unknown what impacts structura l modifications and operational adjustments in surface water management that are currently occurring under the Everglades restoration would have on agricultural land use in farmlands east of C 111. The Biscayne aquifer and the extensive C 111 canal network located in south Florida are connected due to the high permeability of the Miami and Fort Thompson limestone formations that compose the aquifer (Bolster et al. 2001; Lal, 2001; Ritter and Muoz Carpena, 2006). Therefore, the proposed incremental rises in canal stage outlined in the C 111 spreader canal project are anticipated to affect water table levels which in turn could influence soil and limestone bedrock water content. The C 111 spreader canal project is one of the 68 components that comprise the Co mprehensive Everglades Restoration Plan (CERP) whose goal is to restore the natural hydrology of ENP (U.S. Army Corps of Engineers [USACP] and South Florida Water Management District [SFWMD], 2011). Operational changes in C 111 canal stage are planned to b e

PAGE 163

163 implemented at the gated spillway located at structure S18C ( Figure 5 1) in form of gated spillway for up to a maximum of 12 cm (U.S. Army Corps and SFWMD, 2009). So me of the risks to agricultural production associated with a raised water table in a humid shallow water table controlled environment include saturated root zone due to capillarity, which could result in anoxic conditions leading to rotting of roots and de ath of plants. In addition, there could be increased risk of prolonged root zone saturation and temporary groundwater flooding due the rapid water table responses to storm events. Earlier studies have observed disproportionate raises in water table elevati ons after intense rainfall (Gillham, 1984; Germann and Levy, 1986; Heliotis and DeWitt, 1987; Kayane and Kaihotsu, 1988). In prior studies, it was observed that the raise in water table was much larger than would be predicted on the basis of water balance calculations using fillable porosity and net groundwater recharge (Kayane and Kaihotsu, 1988). Germann and Levy (1986) attributed the rapid rise in water table elevation in response to precipitation to capillary fringe groundwater ridging in which a small addition of water to the capillary fringe resulted in a rapid and large rise in water table elevation that drops immediately after the storm. In a recent study Waswa et al. (2013) used field and laboratory observations to show that rapid water table respo nse is caused by rapid addition of energy into the capillary fringe from intense rainfall. Ritter and Muoz Carpena (2006) observed rapid water table rise as spikes in water table elevation observations after storms in the C 111 basin in south Florida and attributed it to high permeability of the aquifer. Under very shallow water table conditions (<1 m from ground surface), these spikes in water table elevation could affect

PAGE 164

164 crops that are intolerant to hypoxia in the root zone or short term flooding. Howeve r, it is not currently known how changes in water table elevation due to the proposed incremental rises in canal stage would impact soil water content in the agricultural areas adjacent to ENP in the C 111 basin. Studies conducted by Barquin et al. (2011) and Kisekka et al. (2013a) also conformed the strong influence of water table elevation on soil and limestone bedrock water content in the C 111 basin. Using the drain to equilibrium assumption Barquin et al. ( 2011) showed that due to capillarity, there was a relationship between soil water content and water table elevation and the later could be used to predict the earlier using the van Genuchten equation (Van Genuchten, 1980) assuming hydrostatic conditions. Kisekka et al. (2013a) using dynamic factor a nalysis showed that water table elevation had a significant influence on temporal variation of soil and limestone bedrock water content. The effect of water table elevation fluctuations on soil water has also been observed by others (Ramirez and Finnerty, 1996). While the study by Kisekka et al. (2013 a ) is empirical in nature, the studies by Barquin et al. (2011) and Ramirez and Finnerty (1996) are based on physics of unsaturated flow under hydrostatic conditions and cannot be applied under transient condit equation are a preferred choice for simulating transient soil water dynamics. The effect of the proposed raises in canal stage on soil and limestone bedrock water dynamics in the a gricultural lands was investigated through simulation of the groundwater vadose zone system. A vadose zone computer code called Water and Agrochemicals in the soil, crop and Vadose Environment (WAVE) developed by Vanclooster et al. (1995) that solves the o

PAGE 165

165 finite difference techniques was applied. WAVE is a comprehensive vadose zone hydrology computer code for simulating the transport of water, energy, non reactive solutes, nitrogen, and pesticides in the soil cro p continuum. WAVE was previously applied in south Florida by Muoz Carpena et al. (2008) to investigate the effect of summer cover crops on soil water retention characteristics and nitrogen leaching in Krome soil. The effect of the proposed incremental rai ses in canal stage on soil and limestone bedrock water content were incorporated by simply linking MODFLOW to WAVE i.e., using water table elevation simulations from MODFLOW as the bottom boundary condition input for WAVE. There are many sources of error which make soil water content measured by indirect methods (e.g. Neutron Moisture Meter Time Domain Reflectometer, Capacitance sensors ) uncertain under all soil conditions (IAIA, 2008) Measurement uncertainty in soil water content data is caused by seve ral factors including : 1) errors related to equipment installation and calibration, 2) errors associated with the different algorithm s that are use d to convert surrogate measurements to soil water content, and 3) errors associated with spatial variab ility of soil properties. However, in many soil water prediction model performance evaluation s ( Whiting et al., 2004; Merdun et al., 2006; Chen et al., 2012 ) goodness of fit indicators like Nash Sutcliffe ( NSE ) Willmot index (d), Root Mean Square Error ( RMSE ) and Mean Absolute Error ( MAE ) based o n calculating the pairwise error between observed and predicted soil water content are used without accounting for the uncertainty in measured data and model structure. A ccurate evaluation of model performance needs to consider these two sources of uncertainty whenever possible in order to provide a more realistic assessment of model

PAGE 166

166 performance Harmel and Smith (2007) provide a framework for quantifying uncertainty in measured data while Harmel et al. (2010) outlines a procedure for quantifying model uncertainty for models in which the predicted state variable can be assumed independent. The purpose of this study was to assess the effect of the proposed operational adjustments in surface water management in form of incr emental raises in canal stage on soil and limestone bedrock water content in farmlands located east of canal C 111. The objectives were to: (1) develop WAVE based models for the study area for simulating soil and limestone bedrock water content (2) evaluat e model performance considering uncertainty in measured soil and limestone bedrock water content and ( 3 ) apply the models to investigate the effect 6, 9 and 12 cm increment raises in canal stage on soil and limestone bedrock water content dynamics at 10, 2 0, 30 and 40 cm depths Material and Methods Study Area and Experimental Set Up The study was conducted in Miami Dade County close to Homestead, Florida within an agricultural area approximately 17 km 2 ( Figure 5 1) immediately to the east of canal C 111 The topography at the project site is essentially flat, implying that the assumption of 1D vertical flow for the unsaturated zone is valid. The soil depth is generally shallow ranging between 10 and 25 cm. Particle size analysis was performed using a sta ndard 2 mm sieve and produced 55% fine particles and 45% gravel. Soil color analysis was performed using the Munsell soil color charts as described in Kisekka et al. (2013 a ). The limestone bedrock layer is highly porous and is generally found at 20 cm dept h and below.

PAGE 167

167 Two m ulti sensor capacitance probes (EnviroScan probes, Sentek Technologies, Ltd., Stepney, Australia) for soil water monitoring were installed at 5 locations ( Figure 5 1) at distances of 500, 1000 and 2000 m from the canal. Soil water conte nt was recorded every 15 minutes but data were averaged daily for modeling purposes. Soil and limestone bedrock water content data were downloaded weekly. EnviroScans are an example of capacitance based sensors which measures frequency of an oscillating el ectrical circuit (IAIA, 2008). Detailed description of how EnviroScans measure soil water content can be found in Kisekka et al. (2013 a ). Two EnviroScan probes were installed at each location to ensure that at least one probe was functioning at any given time. Where applicable (based on quality of the data) data were also averaged at each depth from the two sensors to get a representative average value. Each probe had four sensors positioned at 10, 20, 30 and 40 cm from the ground surface ( Figure 5 2). Th e top 20 cm typically represent the scarified soil layer which is used for crop production and the lower 20 cm represent the underlying limestone bedrock in which plant roots cannot penetrate. To minimize the problem of air pockets between the EnviroScan a ccess tube and the soil matrix, we used fast setting cement slurry. Due to the shallowness of the limestone bedrock at all the study sites, a motorized drill was required to bore a hole that held the access tube in a vertical position. Attempts to calibra te the EnviroScans in the field using the standard gravimetric sampling approach (Sentek Pty Ltd, 2001) was tried but later abandoned due to the difficulty caused by several factors including a very thin soil layer (less than 20 cm) at all the sites; diffi culty obtaining soil samples adjacent to the EnviroScan access tube without interfering with the operation of the sensors; difficulty in obtaining a wide range

PAGE 168

168 of soil water content under field conditions to properly calibrate the sensors; and presence of a shallow limestone bedrock in which it was difficult to sample. Gabriel et al. (2010) in a field study compared default and calibrated volumetric soil water content from EnviroScans and concluded that although default estimates of volumetric soil water co ntent from EnviroScans did not reproduce the exact soil water content they were accurate in reproducing soil water dynamics. Numerical modeling with WAVE WAVE was used to simulate soil wat er dynamics at the different study sites. A unit area soil profile was assumed at each monitoring site and the depth was set between 200 and 220 cm to account for the variations in depth to the water table at the various locations. The soil profile comprised of two distinct soil layers, the first 10 to 20 cm (layer 1) we re set to a scarified soil layer while the remaining 180 to 200 cm (layer 2) were set to limestone bedrock ( Figure 5 2). The profile was discretized into 5 cm compartments and a numerical solution was obtained at the center of each of the compartment ( Figu re 5 2). The time settings in WAVE allowed specifying a variable time step for simulation of highly nonlinear unsaturated flow. The minimum and maximum time steps were set to 0.01 and 1 day, respectively. Due to the problems associated with unattended auto matic data logging of field instruments such as malfunction and flooding of sensors at our study site, time series of observed volumetric water content were not uniform at all the sites. At each site, half of the data was used in model calibration while t he other half was used for model validation. Data at site well 5 were not used in further analysis due to short circuiting of the RLC circuit of the EnviroScan and subsequent installation did not produce quality data, however water table elevation

PAGE 169

169 measured at this site was used in developing a groundwater flow model for the study area. In WAVE, unsaturated flow is simulated with a sink term to account for water uptake by plant roots ( Vanclooster et al. 1995 ) Maximum root wate r uptake is limited by a dimensionless reduction function proposed by Feddes et al. ( 19 78 ), which expresses the effect of pressure head on water uptake rate. A linear relationship was assumed between the reduction factor and root water uptake. The van Genu chten (1980) equation was used to describe the soil water retention characteristics The parameters of the van Genuchten equation were obtained in the laboratory using study site soil samples and plate and a curve fitting c ompute r tool called RETC. The unsaturated hydraulic conductivity was described using the van Genuchten Mualem model (van Genuchten, 1980) and starting values in the modeling were obtained from Muoz Carpena et al. (2008) WAVE solves the nonlinear Richards Nicolson finite difference scheme at the center of each compartment or node ( Vanclooster et al. 1995 ) Since flow is 1D, only the top and bottom boundary conditions were required. At the bottom, a groundwater table boundary cond ition was used. The time series of groundwater table depth which is required by WAVE to describe location of the lower boundary condition was calculated as the difference between the water table elevation (National Geodetic Vertical Datum of 1929 [NGV D 29] m) and elevation of the ground surface (NGVD 29 m). The water table elevation was simulated using MODFLOW as described chapter 4 In WAVE, the top maximum ponding depth was set at 0 mm while

PAGE 170

170 the maximum allowable pressure head (corresponding to air potenti al) at the surface was set as 980 MPa. The stress modeled was water uptake by plant roots. Although detailed crop information was not available to use the crop growth model SUCROS in WAVE, we were able to simulate root water uptake by simply describing a Leaf Areas Index (LAI), root growth depth, and crop coefficient (Kc) time series for the crop grown. We assumed that the crop grown was sweet corn since it represents one of the crops grown at the study site. Developmental and phonological characteristics of south Miami Dade County sweet corn were obtained from Muoz Carpena et al. (2008). Reference evapotranspiration data for the simulation period was obtained from the Florida Automated Network (FAWN; http://fawn.ifas.ufl.edu/) station at TREC, Homestead. At the study site, two thirds of the rainfall is received in the wet season (May to October) and approximately one third is received in the growing or dry season which includes November to April (Pathak, 2001) with October being a transitional month To m inimize the effects of spatial variability in rainfall, we used gauge adjusted NEXRAD rainfall data set at a resolution of 2 km x 2 km from SFWMD. For this study, irrigation was not included in the simulation due to lack of information on irrigation rates and amounts, although we believe irrigation has a substantial effect on soil water dynamics especially during the growing season in the top soil layer. Parameterization and Sensitivity Analysis Soil water retention characteristics for layer 1 ( Table 5 1) were obtained in the laboratory. Water retention characteristics of limestone bedrock were not directly measured due to the extremely porous nature of the material (Muoz Carpena et al., 2008). We estimated saturated soil water content (ts) for the limest one bedrock to range

PAGE 171

171 between 0.24 and 0.41 depending on the site based on readings from EnviroScan sensors at 30 and 40 cm depth s when the sensors were below the water table. Initial literature values for the other parameters i.e., residual limestone water content (tr) and empirical parameters n and a for the limestone bedrock were obtained from literature (Muoz Carpena et al., 2008). Initial pore connectivity parameter (lam) values in (Table 5 1) were also obtained from literature (Mualem, 1976) but were assumed to vary between to 0.5 and 1.5 in our study. Initial parameter values were manually adjusted within ranges established from measurement or literature in order to improve the fit between measured and predicted water content considering uncertainty i n both One crop grown in the study area was s weet c orn. We assumed that three cycles of corn were grown between October and April with no crops being grown May to September as a test crop LAI was set to vary from 1.0 to 2.9 based on 2012 measurement of L AI in a corn field adjacent to site 1. Based on field observations, sweet corn roots grow to a depth of 15 to 20 c m, in WAVE maximum root depth was set to vary from 5 to 20 c m from emergent date to harvest. The time series for crop coefficient Kc ranged fr om 0.6 to 1.1 based on measured values in sweet corn by Muoz Carpena et al. (2008) for Homestead, Florida (Table 5 2 to 5 4) The date at which roots become inactive due to senescence was set to harvest date. The relationship between the reduction factor of root water uptake and pressure head was assumed linear and the pressure head at which roots start to extract water and optimally extract water from the soil were set at 10 and 46 cm of water, respectively (i.e., pressure heads greater than 10 cm woul d limit oxygen levels in the soils through soil saturation and thus affect plants roots and consequently plant water uptake). The

PAGE 172

172 pressure head at which water uptake by the root ceases was set at 16,000 cm which corresponds to the theoretical permanent wi lting point. Sensitivity analysis was implemented in two stages: 1) the improved Morris method by Campologo et al. (2007) was performed to obtain qualitative ranking of parameters (i.e., screening parameters to which the simulated volumetric soil water con tent was most sensitive) and 2) using a subset of critical parameters from step 1, (parameter interaction) sensitivity indices. Parameters included in sensitivity analysi s for layers 1 and 2 at different sites are given in Table 5 1. For all the parameters with the exception of LAI and Kc, a uniform distribution was assumed and parameters ranges were obtained from measurements o r from literature. In order to test the sensi tivity of simulated volumetric soil water content to variations in LAI, a discrete uniform distribution was assumed in which three values representing LAI during initial plant development, mid season, and late season stage were used. LAI values in Tables 5 2 to 5 4 are based on measurements made in a sweet corn field at the current study site. A similar approach was used for testing sensitivity of simulated volumetric water content to variations in Kc. During sensitivity analysis, input parameter sets were obtained by sampling the multivariate input distributions according to Campologo et al. (2007) improved Morris was then executed for each input parameter set. Campologo et al. (2007) sensitivity analysis was implemented using Matlab algorithms (R2012a Mathworks Inc., Natick, Massachusetts) developed by Saltelli et al. (2008) and available online at

PAGE 173

173 http://sensitivityanalysis.jrc.it/software/index.htm. Matlab was used to au tomatically executed WAVE for each parameter set in the generated sample input file. The sample input file was a matrix comprised of multiple input parameter sets obtained from random sampling of probability distributions. For sample generation using Campo longo et al. (2007) method, the following settings were used: number of levels ( p ) was 4, size of oversampling ( N ) was 1000, number of trajectories ( r ) was 20 and number of parameters ( k ) was 19. Th i s resulted in a total of 400 parameter sample set s (i.e. ). The number of WAVE executions for Sobol analysis was estimated as w h ere the sample size n was 512 and k is the number of the identified critical parameters from Campologo et al. (2007) analysis. Nash Sutc liffe coefficient (NSE) and the Root Mean Square Error (RMSE) were calculated as the model output for each simulation. Estimating uncertainty in measured soil and bed rock water content Soil and bedrock water content measured in this study was considered uncertainty due to the following reason s : 1) sensor specific calibration was not performed due to having a very thin soil layer under layered by hard limestone bedrock which made collection of samples for gravimetric estimation of soil water difficulty, 2) spatial variability in soil physical properties, and 3 ) variations in EnviroScan installations specifically the thickness of the cement slurry surrounding the access tube. To accommodate all these sources of error, we quantified uncertainty for each measu red soil and bedrock water content value. Uncertainty in measured soil and limestone water content data was calculated using a correction factor (Harmel and Smith, 2007; Harmel et al., 2010) The correction factor modifies the error term (i.e., the pair wi se difference

PAGE 174

174 between measured and predicted values) in goodness of fit indicators by incorporating the distribution of the measurement uncertainty as shown in Equation 5 1 (5 1) where is the modified deviati on between measured and predicted soil and bedrock water content for point considering only measurement uncertainty, is the non dimensional correction factor (ranges between 0 and 1) for each measured soil and bedrock water content ( ) and predicted soil and bedrock water content ( ) considering measurement uncertainty, and 0.5 refers to one sided probability for ( ) at mean value assuming a symm etric distribution WAVE was calibrated by manually adjusting parameter values within the ranges in Table 5 1. These ranges were selected based on laboratory measurements and represented the range in values assessed for each parameter. Thus, these ranges were a best estimate of parameter uncertainty. Parameter values were adjusted until the predicted soil water content was within the maximum and minimum uncertainty bounds of the measured data as calculated by Equations 5 2 to 5 3 5 2 5 3 where the uncertainty bounds and are the lower and upper bounds of the uniform distribution, is the mean of the distribution for measurement set at the measured value and Cv is coefficient of variation. In this study we assumed a uniform distribution for all measurements and minor (coefficient of variation, Cv =0.02) to moderate

PAGE 175

175 ( Cv =0.08) dep ending how variable that data collected from two adjacent EnviroScan probes was. After calibration, model performance (validation) was implemented using FITEVAL (Ritter and Muoz Carpena, 2012). Validat ion was performed in two stages: 1) without consideri ng uncertainty in measured values and 2) accounting for uncertainty in measured values following procedure s describe d in Harmel et al. (2010). Model uncertainty was not included in this analysis Model Applications The validated WAVE based model s at each of the five sites were applied to simulate soil and limestone bedrock water content at different depth s under the proposed 6, 9 and 12 cm incremental raises in canal stage. Canal stage will be increased in the reach of canal C 111 between the gated spillwa ys at S177 and S18C (Figure 5 1) gated spillways (U.S. Army Corps of Engineers [USACP] and South Florida Water Management District [SFWMD], 2011). Effect of surface water managem ent on water table elevation was simulated using MODFLOW as described in chapter 4 The simulated water table elevation was then used as a lower boundary condition of the soil profile for WAVE Increases in canal stage were mimicked by incrementally adding 3 cm to the current canal stage up to a total of 12 cm. Results and Discussion Sensitivity Analysis and Parameterization of Vadose Zone Model P arameter screening results ( Tables 5 5 to 5 9 ) indicate that predicted soil water content was most sensitive to parameters of the van Genuchten equation which was the water retention model used in WAVE i.e., residual soil water content, saturated soil

PAGE 176

176 water content, curve shape parameter and the inverse of the air entry value This would be expected because soil water retention curve parameters characterize soil water retention in both soil and limestone bedrock layers. Other studies have also demonstrated that state variables predicted by WAVE are highly sensitive to soil hydraulic parameters ( Vanclooster et al., 1995; Muoz Carpena et al., 2008). At our study site the predicted soil water content showed only slight to no sensitivity to variations in all the other parameters including variations in Kc and LAI. Based on the results from improved Morris screening a subset of important WAVE parameters were identified and used in further variance based global sensitivity ( Figure s 5 3 to 5 7 ) Sobol analysis confirm Morris screening results indicating that saturated soil water content i s the most important parameter explaining variations in predicted soil water content as measured by NSE and RMSE at all sites. It should be noted that Sobol analysis is a more robust technique for sensitivity analysis compared to the Morris method because it is based on a large number of model simulations and less structured sampling ( Saltelli et al., 2004 ). T he fraction of the total variation in predicted soil and limestone bedrock water content explained by variation in each of the ten important parameter s is represented using first order and total order Sobol sensitivity indices along the vertical axis (Figures 5 3 to 5 7) The first bar represents first order effects while the second represents total order effects and the difference between the two bars represents parameter interactions. Results also show that unsaturated flow is influenced by the parameters differently in the soil and limestone bedrock layers. In the soil layer (top 20 cm), unsaturated flow was mainly governed by saturated and residual s oil water content, curve shape parameter, and

PAGE 177

177 inverse of the air entry value and the effects of parameter interactions were stronger than in the limestone bedrock layer. In the limestone bedrock layer (30 and 40 cm depth) unsaturated flow was mainly gover ned by saturated soil water content and the first order effects approached 100% indicating that WAVE behaved as an additive model within the limestone bedrock layer particularly at sites 4 and 6 w h ere sensors at 30 and 40 cm were surrounded by limestone wa ter content conditions close to saturation for the majority of the study. This is probably due to the fact that the conditions (Vanclooster et al. 1995). WAVE behaving as an a dditive model at 30 and 40 cm depth indicated that it could be calibrated using accurately measured soil and limestone bedrock water content data. WAVE was calibrated by setting parameters to which predicted soil water content was least sensitive to aver age values and manually adjusting important parameters (as identified in sensitivity analysis) within the ranges listed in Table 5 1 until we graphically achieved an acceptable fit between measured and predicted soil water content considering uncertainties in each (Table 5 10) Automated calibration using the GLUE methodology was tried but later abandoned in order to avoid over fitting model parameters to an uncertainty dataset of measured volumetric soil and limestone bedrock water content. Estimated satur ated soil and limestone bedrock water content w ere compared to saturated soil and limestone bedrock water content when the EnviroScan sensors were below the water table and the values were in close agreement For example at site 4 saturated soil water con tent ( when the sensor were below water table ) was identified as 35 (m 3 /m 3 ) and the manually estimated values for

PAGE 178

178 layers 1 and 2 were 35 and 34 respectively. We attributed difficulty of graphically achieving a perfect fit between measured and predicted s oil water content at various depth at the same site and across the different five sites to the following factors: 1) uncertainty in measured data, 2) intrinsic spatial variability in soil and limestone hydraulic parameters, and 3) exclusion of irrigation w ater applied from the conceptual model Soil W ater P rediction C omparison between simulated and measured volumetric water content from EnviroScan probes at 10, 20, 30, and 40 cm depths under current canal stage operation criteria along C 111 were plotted ( Figures 5 9 to 5 13) Visual inspection indicates that WAVE was able to reproduce temporal variations in soil water content as influenced by seasonal variations in rainfall, evapotranspiration and canal stage. It can also been seen that there were some su bstantial deviations between predicted and measured volumetric water content at some sites and monitoring depth particularly during the summer of 2011 months (April to October) which also corresponded to the lowest recorded soil water content. These deviat ion s could be attributed to several factors discussed under the section on estimating measurement uncertainty. Other investigators ( Evett et al., 2006; Rowland et al., 2011) have noted that technologies that measure soil water content by sensing electromag netic properties of the soil such as EnviroScans are prone to uncertain precision and sensed volumes issues. For example EnviroScan s measure an effective distance of only 3 5 cm from the access tube and may be affected by non isothermal conditions althoug h they are idea for unattended continuous monitoring of soil water content.

PAGE 179

17 9 Deviation between predicted and measured soil water content at site 3 during the first months of the study was due to poor installation which was subsequently re installed thus im proving data at this site. It is worth noting that transformation of measured data using the EnviroScan calibration equation developed in the laboratory by Al Yahyai et al. (2006) for gravely loam soils of south Florida was tried but gave inconsistent resu lts at various depths and sites and was abandoned. This could be due to the differences in soil conditions at our study site compared to the soil Al Yahyai et al. (2006) used to develop the calibration equation in the laboratory. Goodness of fit statistic s for model validation for the different sites and monitoring depth without and with con si deration of measurement uncertainty were calculated ( Table s 5 11 and 5 12 ) Fit between EnviroScan measured soil and bedrock water content and simulated soil and bedr ock water content were unsatisfactory ( based on goodness of fit classification in Ritter and Muoz Carpena 2012 ) for all sites and depth s with the exception of 30 and 40 cm depth s at site 1 when uncertainty in measured soil water content was not taken int o account (Table 5 11) However, when uncertainty in measured water content data was considered th ere was a n improvement in the goodness of fit statistics at all sites except site 2 (Table 5 12). At site 2 only measurement s at depth 30 cm were within range of soil water content generally expected within the study area, for example measurements at 40 cm were consistently lower than measurements at 10, 20, and 30 cm depth s at the same site even when water table elevation was close to the surface, implying tha t measurements at site 2 were very inaccurate. Therefore, data at this site were not used in further evaluations

PAGE 180

180 Therefore, we can conclude that the model s developed for the five sites in this study are useful for estimating soil water content and for pr edicting temporal variations in soil and bedrock water content under changing surface water management which was the purpose for under taking the study. Evaluation of Soil Water Response to Proposed Incremental Raises in Canal Stage Effect of the proposed changes in canal C 111 stage operation criteria were evaluated for only three sites south of the spillway at structure S177 since these are the sites expected to experience significant changes in water table elevation due to raised canal stage (Table 5 16) because they are located at lower elevation within the study area Canal stage is expected to be raised along the reach of C 111 between S177 and S18C (Figure 5 1). Results indicate that the soil and bedrock water content before and after the proposed cha nges in canal stage were significantly different (p<0. 00 6 ). The i mpact of raises in canal stage on soil and bedrock water content over time were mainly seen when there were spikes in soil and bedrock water content approaching saturation water content prob ably corresponding to precipitation events or raises on canal stage ( Figures 5 15 to 5 17 ) At site 3 after the proposed raises in canal C 111 stage, there were no substantial differences in soil water content both during the wet and dry seasons ( Figure 5 15 ). I n the top 20 cm soil layer, soil water content did not reach saturation even after the maximum proposed increment in canal stage of 12 cm. Implying that farmlands with ground surface elevation similar to that at site 3 i.e., greater than 2.0 m NGVD2 9 are predicted not experience root zone saturation after the proposed incremental raises in canal stage. However, water content started to approach saturation in the limestone layer after the rains received during the wet

PAGE 181

181 season of 2012 Although t his co ndition is not expected to hinder aeration in the root zone since the roots of the crops grown in this area never penetrate the limestone bedrock saturation at 30 and 40 cm depth might exacerbate the problem of temporary groundwater flooding due to spikes in water table elevation associated with the phenomenon of groundwater ridging. At site 4, changes in canal stage did result in observable changes in soil water content both during the wet season and dry season (Figure 5 16) S oil water content reached sa turation at 20 cm depth after increasing canal stage by 9 cm and this condition persisted up to late January of 2012 ( Figure 5 16 ), implying that crop production growing periods would be greatly reduced. Soil water content in the limestone approached satu ration even under current canal stage operation criteria but raising canal stage resulted in almost continuous limestone saturation at the 30 and 40 cm monitoring depth s These results mean that farmlands with land surface elevation similar to that of site 4 (1.19 m NGVD29) might be impacted by increases in canal stage greater than 9 cm. The response at site 6 was similar to that observed at site 4 probably due to similar elevation (1.2 m NGVD29) in which the changes in canal stage resulted in increases in soil and bedrock water content at all depth. However, for the top 1 0 cm the increase in soil water content did not lead to saturated conditions during the growing season but soil water content approached saturation during the wet season. The changes in can al stage resulted in saturated conditions at 30 and 40 cm during both the wet and dry season (Figure 5 17) Based on the period ( January 2012 to February 2013 ) investigated for potential impacts of raising canal stage on root zone soil water content the p roposed changes in C 111 canal stage resulted in significant changes in soil water content. Sites with land surface elevation

PAGE 182

182 greater than 2.0 m NGVD29 did not experience saturated conditions in the top 20 cm soil layer. Raising canal stage by more than 9 cm resulted in saturated root zone and shortening of the growing season at sites with land surface elevation less than 2.0 m NGVD29, which is critical for continued use of the land for agricultural production. These results are similar to those obtained using a different technique called dynamic factor analysis using data from the same study area (Kisekka et al., 2013 a ). The benefit of the WAVE based mechanistic approach for the study area compared to the DFA based approach used in Kisekka et al (2013 a ) is the former can be applied for exploring soil water dynamics outside the ranges of data measured during the study period and also be used to examine the physical processes govern ing the response of the state variables e.g., soil water content or transpi ration within the ranges of uncertainty associated with model predictions Conclusion One dimensional WAVE based model s w ere developed for the study area for mechanistically investigating the effect of the proposed incremental raises in canal C 111 stage o n soil and bedrock water content in the top 40 cm of the soil and limestone layer. Parameter screening using the improved Morris method indicated that predicated soil water content was most sensitive to parameters of the van Genuchten equation. Quantitativ from Morris screen by showing that soil hydraulic properties contributed the greatest variance in predicted soil water content. Specifically, s aturated soil water content was ident ified as the most important input factor affecting predicted soil and bedrock water content at our study area The model behaved nonlinearly in the top 20 cm with various

PAGE 183

183 parameter interactions, and approximated an additive (linear) model in the usually sa turated limestone layer. Goodness of fit indicators were greatly improved when uncertainty in measured data was considered during model performance evaluation. There were some variations in model performance at different sites and monitoring depth these c ould be attributed to varying accuracy in measured data These differences were attributed to variations in thickness of the cement slurry surrounding the access tube, intrinsic spati al variability in soil and bedrock hydraulic properties, effect of not in cluding irrigation water inputs and uncertainties in associated with variations in installations Overall model performance ranged from good to acceptable at all the sites analyzed. We concluded that WAVE based models developed for the study area with the exception of the model developed for site 2 were useful as exploratory tools for estimating changes in soil water content conditions in the top 40 cm depth and for predicting temporal variations in soil water content in response to surface water managemen t in the C 111 canal. The RMSE and MAE produced by the model were regarded low for practical soil water management for crop production. The model was applied to evaluate the effect of the 6, 9 and 12 cm incremental raises in canal stage on soil and bedro ck water content. The results indicated that average soil water content before and after the changes in C 111 canal stage were significantly different. Visual analysis of the soil water content time series under proposed changes in canal stage management indicated that sites with land surface elevation of less than 2.0 m NGVD29 might experience saturated conditions in the root zone and shortening of the growing season for agricultural production At depths

PAGE 184

184 greater than 20 cm raises in canal stage resulted in saturated conditions particularly for sites with low land surface elevation and during the wet season. Th e saturat ed conditions at the 30 and 40 cm depth at all the sites could exacerbate the problem of temporary groundwater flooding due to groundwater ridging suggesting that water management practices (i.e., pre storm releases) would need to be modified The model s developed in this study could be could be combined with high resolution digital elevation model s (DEM) in future studies to identify areas that should not be planted to minimize potential losses associated with saturated root zones. This study shows that in humid water table controlled landscapes, surface water management could impact soil water content which in turn influences land use and other near land surface water and energy fluxes that are influenced by near ground soil water content

PAGE 185

185 Table 5 1 Initial p arameter ranges used in WAVE for simulating soil water content at five sites within the C 111 basin assuming a uniform distribution for all parameters Description Parameter Value Layer 1 (Soil layer) Saturated soil water content (m 3 m 3 ) 1 ts1 0. 2 7 0. 53 Residual soil water content (m 3 m 3 ) 2 tr1 0.0 0 0.092 Inverse of the air entry value (cm 3 ) 2 a1 0. 00 1 0.093 Curve shape parameter (slope) ( ) 2 n1 0.5 1.5 1 Pore connectivity parameter ( ) 3 lam1 0. 5 1. 50 Unsaturated hydraulic conductivity (cm/day) 4 K1 500 1551 Maximum water uptake rate (day 1 ) 4 Smax1 0.01 0.014 Layer 2 (Limestone) Saturated soil water content (m 3 m 3 ) ts2 0.20 0. 46 Residual soil water content (m 3 m 3 ) 4 tr1 0.0 0.01 Inverse of the air entry value (cm 3 ) 4 a2 0.009 0.15 Curve shape parameter (slope) ( ) 4 n2 0.9 1.5 Pore connectivity parameter ( ) lam2 0. 5 0 4.5 Unsaturated hydraulic conductivity (cm/day) 4 K2 5000 14000 1 Measured in laboratory 2 Values obtained from fitted curve of soil water suction data measured in laboratory 3 Value s obtained from Mualem (1976) 4 Value s obtained from Muoz Carpena et al. (2008) Table 5 2. Crop coefficient (Kc) and leaf area index (LAI) for sites 3 and 6 Station Kc 1 Kc value Kc symbol LAI 2 LAI symbol 01/28/2011 Initial 0.6 Kc1 0.5 LAI1 03/15/2011 Mid season 1.1 Kc2 2.9 LAI2 05/15/2011 Late season 0.85 Kc3 1.45 LAI3 1 Crop coefficient value obtained from Muoz Carpena et al. (20 08) 2 Leaf area index value measured

PAGE 186

186 Table 5 3. Crop coefficient (Kc) and c rop leaf area index (LAI) for site 1 Station Kc 1 Kc value Kc symbol LAI 2 LAI symbol 04/21/2011 Late season 0.85 Kc1 1.45 LAI1 11/01/2011 Initial 0.6 Kc3 0.5 LAI3 12/31/2011 Mi d season 1.1 Kc4 2.9 LAI4 1 Crop coefficient value obtained from Muoz Carpena et al. (2008) 2 Leaf area index value measured Table 5 4. Crop coefficient (Kc) and c rop leaf area index (LAI) for site s 2 and 4 Station Kc 1 Kc value Kc symbol LAI 2 LAI symbol 10/31/2010 Initial 0.6 Kc1 0.5 LAI1 12/15/2010 Mid season 1.1 Kc2 2.9 LAI2 04/30/2011 Late season 0.85 Kc3 1.45 LAI3 1 Crop coefficient value obtained from Muoz Carpena et al. (2008) 2 Leaf area index value measured

PAGE 187

187 Table 5 5. Improved Morris screenin g results for WAVE model applied at site 1 for soil layer (L1 : 10 and 20 cm ) and limestone layer (L2 : 30 and 40 cm ) NSE 1 10 NSE20 NSE30 NSE40 Parameter 2 3 Residual water content L 1 (m3/m3) 1 9 2 3 4 3 5 0 0. 1 0. 1 0 0 Residual water content L 2 (m3/m3) 0 0 0. 6 0. 6 2 9 2 7 0.1 0.1 Saturated water content L 1 (m3/m3) 7 4 5 5 22 6 17 5 1 3 0 0 Saturated water content L 2 (m3/m3) 0 0 15 7 10 5 93 0 54 2 38 7 24 2 Inverse of air entry value L 1 (cm 1) 2 9 3 4 9 9 11 9 6 9 1 1 Inverse of air entry value L 2 (cm 1) 0. 6 0. 9 16 4 12 2 67 0 56 4 25 7 21 8 Curve shape parameter L 1 ( ) 5 7 4 4 16 1 13 0 0. 6 0. 9 0. 1 0.1 Curve shape parameter L 2 ( ) 0. 4 0. 7 16 9 19 3 69 9 31 6 25 4 12 3 Sat hydraulic cond uctivity L 1 (m/day) 0. 1 0.1 0.1 3 0. 7 0.3 0. 6 0 0. 1 Sa t. hydraulic conductivity L 2 (m/day) 0. 3 0. 7 1 1 2 0 0. 7 1 2 0. 2 4 .0 Pore connectivity parameter L 1 ( ) 0 0. 1 0. 3 0. 7 0. 1 0. 4 0. 0 0. 1 Pore con nectivity parameter L 2 ( ) 0. 2 0. 5 1 1 1 9 0. 8 1 3 0. 2 0. 5 Crop coefficient initial stage 0 0. 1 0. 1 0. 1 0 0. 1 0 0 Crop coefficient mid season stage 0 0. 1 0. 1 0. 1 0. 1 0. 1 0 0 Crop coefficient late season stage 0 0 0 0 0 0 0 0 Leaf area index initial stage 0 0 0 0 0 0 0 0 Leaf area index mid season stage 0 0 0 0 0 0 0 0 Leaf area index late season stage 0 0 0 0 0 0 0 0 Maximum root water uptake 0 0 0 0 0 0 0 0 1 Nash Sutcliffe Coefficient of efficient was used as model output in sensitivity analysis corresponding to each model run. 2 Absolute values of the Morris mean which assesses the overall effect of the factor on model output 3 factors

PAGE 188

188 Table 5 6. Improved Mor ris screening results for WAVE model applied at site 2 for soil layer (L1: 10 and 20 cm) and limestone layer (L2: 30 and 40 cm) NSE 1 20 NSE30 NSE40 Parameter 2 3 Residual water content L 1 (m3/m3) 2 .0 2 .0 1 .0 1 .0 0 0 Residual water content L 2 (m3/m3) 0 0 0 0 0 0 Saturated water content L 1 (m3/m3) 12 .0 7 .0 7 .0 5 .0 0 0 Saturated water content L 2 (m3/m3) 0 0 7 6 18 .0 15 .0 Inverse of air entry value L 1 (cm 1) 5 .0 6 .0 2 .0 3 .0 0 0 Inverse of air entry value L 2 (cm 1) 1 .0 1 .0 4 .0 4 .0 9 .0 9 .0 Curve shape parameter L 1 ( ) 12 .0 9 .0 6 .0 6 .0 0 0 Curve shape parameter L 2 ( ) 1 .0 0 7 .0 6 .0 16 .0 17 .0 Sat hydraulic conductivity L 1 (m/day) 0 0 0 0 0 0 Sat hydraulic conductivity L 2 (m/day) 0 1 .0 0 1 .0 0 0 Pore connectivity parameter L 1 ( ) 0 0 0 0 0 0 Pore connectivity parameter L 2 ( ) 1 .0 1 .0 0 1 .0 0 0 Crop coefficient initial stage 0 0 0 0 0 0 Crop coefficient mid season stage 0 0 0 0 0 0 Crop coef ficient late season stage 0 0 0 0 0 0 Leaf area index initial stage 0 0 0 0 0 0 Leaf area index mid season stage 0 0 0 0 0 0 Leaf area index late season stage 0 0 0 0 0 0 Maximum root water uptake 0 0 0 0 0 0 1 Nash Sutcliffe Coefficient of efficient was used as model output in sensitivity analysis corresponding to each model run. 2 Absolute values of the Morris mean which assesses the overall effect of the factor on model output 3 ons with other factors

PAGE 189

189 Table 5 7 Improved Morris screening results for WAVE model applied at site 3 for soil layer (L1: 10 and 20 cm) and limestone layer (L2: 30 and 40 cm) NSE120 NSE30 NSE40 Parameter 2 3 Residual water content L 1 (m3/m 3) 1.1 1.2 0.0 0.1 0.0 0.1 Residual water content L 2 (m3/m3) 0.1 0.1 0.9 0.7 0.7 0.6 Saturated water content L 1 (m3/m3) 2.7 2.2 0.1 0.1 0.0 0.1 Saturated water content L 2 (m3/m3) 2.1 1.5 19.1 10.5 17.2 8.9 Inverse of air entry value L 1 (cm 1) 1.0 1.3 0.6 1.0 0.3 0.5 Inverse of air entry value L 2 (cm 1) 3.3 2.9 15.2 12.8 12.9 11.2 Curve shape parameter L 1 ( ) 3.9 3.2 0.4 0.7 0.2 0.4 Curve shape parameter L 2 ( ) 3.7 3.3 18.2 10.5 15.0 7.8 Sat hydraulic conductivity L 1 (m/day) 0.1 0.1 0.2 0.2 0.1 0.1 Sat hydraulic conductivity L 2 (m/day) 0.3 0.3 0.5 0.7 0.4 0.6 Pore connectivity parameter L 1 ( ) 0.1 0.2 0.1 0.1 0.1 0.1 Pore connectivity parameter L 2 ( ) 1.0 1.8 0.9 1.3 0.7 1.1 Crop coefficient initial stage 0.0 0.0 0.0 0.0 0.0 0.0 Crop coe fficient mid season stage 0.0 0.1 0.0 0.1 0.0 0.0 Crop coefficient late season stage 0.1 0.1 0.1 0.1 0.0 0.0 Leaf area index initial stage 0.0 0.0 0.0 0.0 0.0 0.0 Leaf area index mid season stage 0.0 0.0 0.0 0.0 0.0 0.0 Leaf area index late season st age 0.0 0.0 0.0 0.0 0.0 0.0 Maximum root water uptake 0.0 0.0 0.0 0.0 0.0 0.0 1 Nash Sutcliffe Coefficient of efficient was used as model output in sensitivity analysis corresponding to each model run. 2 Absolute values of the Morris mean which assesses th e overall effect of the factor on model output 3 factors

PAGE 190

190 Table 5 8. Improved Morris screening results for WAVE model applied at site 4 for soil layer (L1: 10 and 20 cm ) and limestone layer (L2: 30 and 40 cm) NSE 1 10 NSE20 NSE30 NSE40 Parameter 2 3 Residual water content L 1 (m3/m3) 4.3 5.7 2.4 2.9 0.0 0.0 0.0 0.0 Residual water content L 2 (m3/m3) 0.0 0.0 0.3 0.4 1.4 1.4 2.0 2.0 Saturated water conten t L 1 (m3/m3) 24.2 16.7 16.0 14.7 0.1 0.1 0.0 0.1 Saturated water content L 2 (m3/m3) 0.0 0.1 11.8 11.4 69.4 49.9 120.1 105.3 Inverse of air entry value L 1 (cm 1) 7.7 6.9 4.5 5.1 0.3 0.3 0.2 0.3 Inverse of air entry value L2 (cm 1) 2.1 2.8 8.9 7.6 33.6 27.7 50.4 43.6 Curve shape parameter L 1 ( ) 14.3 10.8 9.1 6.7 0.2 0.4 0.2 0.3 Curve shape parameter L 2 ( ) 0.9 1.3 9.7 11.9 31.8 16.0 45.7 26.0 Sat hydraulic conductivity L 1 (m/d) 0.2 0.5 0.1 0.2 0.2 0.3 0.1 0.2 Sat hydraulic conductivity L 2 (m/d) 1. 1 2.5 1.3 2.5 0.7 0.8 0.5 0.7 Pore connectivity parameter L 1 ( ) 0.1 0.3 0.1 0.3 0.1 0.1 0.1 0.1 Pore connectivity parameter L 2 ( ) 0.7 1.1 1.0 1.5 1.0 1.4 0.9 1.1 Crop coefficient initial stage 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Crop coefficient mid sea son stage 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Crop coefficient late season stage 0.4 0.7 0.2 0.4 0.1 0.2 0.1 0.1 Leaf area index initial stage 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Leaf area index mid season stage 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Leaf area inde x late season stage 0.1 0.2 0.0 0.1 0.0 0.0 0.0 0.0 Maximum root water uptake 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1 Nash Sutcliffe Coefficient of efficient was used as model output in sensitivity analysis corresponding to each model run. 2 Absolute values of t he Morris mean which assesses the overall effect of the factor on model output 3 factors

PAGE 191

191 Table 5 9. Improved Morris screening results for WAVE model applied at site 6 f or soil layer (L1: 10 and 20 cm) and limestone layer (L2: 30 and 40 cm) NSE 1 10 NSE20 NSE30 NSE40 Parameter 2 3 Residual water content L1 (m 3 /m 3 ) 1.3 1.6 0.7 0.7 0.0 0.0 0.0 0.0 Residual water content L 2 (m 3 /m 3 ) 0.0 0.0 0.1 0.1 0.5 0.4 0.7 0.6 Saturated water content L 1 (m 3 /m 3 ) 10.7 7.0 9.3 6.4 0.0 0.0 0.0 0.0 Saturated water content L 2 (m 3 /m 3 ) 0.0 0.0 11.3 6.0 39.2 25.5 84.7 52.0 Inverse of air entry value L 1 (cm 1 ) 4.7 4.7 2.5 2.5 0.1 0.1 0.1 0.1 Inverse of air entry value L 2 (cm 1 ) 0.7 1.0 3.8 3.6 10.8 9.8 16.7 17.1 Curve shape parameter L 1 ( ) 7.8 5.9 4.0 4.4 0.1 0.1 0.0 0.1 Curve shape parameter L 2 ( ) 0.4 0.6 5.0 4.9 13.8 13.0 21.5 21.2 Sat hydraulic conductivity L 1 (m/day) 0.2 0.2 0.0 0.1 0.0 0.1 0.0 0.0 Sat hydrau lic conductivity L 2 (m/day) 0.6 1.0 0.4 0.8 0.3 0.3 0.2 0.2 Pore connectivity parameter L 1 ( ) 0.0 0.1 0.0 0.1 0.0 0.0 0.0 0.0 Pore connectivity parameter L 2 ( ) 0.3 0.6 0.2 0.4 0.2 0.3 0.1 0.2 Crop coefficient initial stage 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Crop coefficient mid season stage 0.0 0.0 0.0 0.1 0.0 0.0 0.0 0.0 Crop coefficient late season stage 0.1 0.1 0.1 0.1 0.0 0.0 0.0 0.0 Leaf area index initial stage 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Leaf area index mid season stage 0.0 0.0 0.0 0.0 0 .0 0.0 0.0 0.0 Leaf area index late season stage 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 Maximum root water uptake 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1 Nash Sutcliffe Coefficient of efficient was used as model output in sensitivity analysis corresponding to each mo del run. 2 Absolute values of the Morris mean which assesses the overall effect of the factor on model output 3 factors

PAGE 192

192 Table 5 1 0 WAVE parameters obtained from calib ration at different sites ( October 1, 2010 to December 31 2011 ) Parameter Site 1 Site 2 Site 3 Site 4 Site 6 Avg. Layer 1 Residual water content (tr1) 0. 09 0.0 9 0. 08 0.10 0.10 0.09 Saturated water content (ts1) 0.3 0 0.2 9 0. 32 0.34 0.30 0.31 Cur ve shape parameter (n1) 1. 09 1. 09 1.2 2 1.17 1.15 1.14 Inverse of air entry value (a1) 0.0 4 0. 1 2 0. 06 0.09 0.09 0.08 Pore connectivity parameter (lam) 0. 50 0. 61 0.62 0.54 0.62 0.58 Layer 2 Residual water content (tr2) 0.0 9 0. 09 0. 0 6 0.09 0.09 0.08 S aturated water content (ts2) 0.3 1 0. 31 0.3 1 0.35 0.29 0.31 Curve shape parameter (n2) 1.1 2 1. 11 1.1 1 1.11 1.10 1.11 Inverse of air entry value (a2) 0. 08 0. 12 0.1 0 0.10 0.09 0.10 Pore connectivity parameter (lam) 0. 50 0. 61 0.62 0.54 0.62 0.58 Sat hydra ulic conductivity (K2) 8000 8335 9307 8511 8419 8514

PAGE 193

193 Table 5 1 1 Goodness of fit statistics without cons id eration of measurement uncertainity for WAVE water content simulations by soil depth during the validation period ( site 1 [ 01/01/2012 to 0 2/ 28 /201 3]; site 2 [ 10/01/2011 to 0 2/ 28 /201 3]; site 3 [ 01/01/2012 to 0 2/ 28 /201 3], site 4 [ 10/01/2011 to 0 2/ 28 /201 3], site 6[ 01/01/2012 to 0 2/ 28 /201 3] Site 1 Depth 10 cm 20 cm 30 cm 40 cm NSE 1 0.26( 0.32 0.66) 0.35(0.03 0.59) 0.80(0.66 0.88) 0.77(0.74 0.88) RMS E 2 0.95(0.72 1.19) 0.62(0.52 0.73) 0.35(0.31 0.41) 0.48(0.41 0.54) A 3 (%) 0.0 0.0 0.5 0.5 B 4 (%) 0.0 0.0 51.4 29.1 C 5 (%) 3.1 0.5 46.4 54.3 D 6 (%) 93.9 99.5 1.7 16.1 Site 3 Depth 10 cm 20 cm 30 cm 40 cm NSE 2.79( 6.21 1.3) 0.57(0.28 0.76) 0.44(0. 27 0.57) 0.10( 0.90 066) RMSE 1.16(0.94 1.37) 0.81(0.71 0.96) 0.84(0.65 1.06) 1.35(0.94 1.72) A (%) 0.0 0.0 0.0 0.0 B (%) 0.0 0.6 0.0 0.0 C (%) 0.0 67.1 0.1 0.0 D (%) 100.0 32.3 99.1 100.0 Site 4 Depth 10 cm 20 cm 30 cm 40 cm NSE 0.25( 0.24 0.50) 0.30( 0.66 0.60) 1.95( 5.5 -0.17) 3.80( 6.20 -0.87) RMSE 1.42(1.16 1.66) 0.63(0.52 0.81) 0.99(0.73 1.22) 0.89(0.55 1.28) A (%) 0.0 0.0 0.0 3.2 B (%) 0.0 0.0 0.0 5.3 C (%) 0.0 2.6 0.0 10.4 D (%) 100.0 97.4 100.0 81.1 Site 6 Depth 10 cm 20 cm 30 cm 40 cm NSE 0.35( 1.80 0.45) 0.31( 0.99 0.66) 0.63( 2.11 0.01) 12.2( 29.6 -4.18) RMSE 0.95(0.63 1.29) 0.52(0.41 0.66) 0.63(0.43 0.80) 0.97(0.79 1.11) A (%) 0.0 0.0 0.0 0.0 B (%) 0.0 0.1 0.0 0.0 C (%) 0.2 5.2 0.0 0.5 D (%) 99.8 94.7 100.0 99.5 1 Nas h Sutcliffe coefficient 2 Root mean square error 3 A probability of fit being very good 0.9
PAGE 194

194 Table 5 1 2 Goodness of fit statistics considering measurement uncertainity for WAVE water content simulations by soil depth during the validation period (site 1 [ 01/01/2012 to 0 2/ 28 /201 3]; site 2 [ 10/01/2011 to 0 2/ 28 /201 3]; site 3 [ 01/01/2012 to 0 2/ 28 /201 3], si te 4 [ 10/01/2011 to 0 2/ 28 /201 3], site 6[ 01/01/2012 to 0 2/ 28 /201 3] Site 1 Depth 10 cm 20 cm 30 cm 40 cm NSE 1 0.78(0.50 0.92) 0.87(0.75 0.93) 0.89(0.75 0.94) 0.85(0.68 0.93) RMSE 2 0.53(0.33 0.74) 0.28(0.20 0.39) 0.26(0.20 0.40) 0.38(0.31 0.51) A 3 (%) 7 .2 23.4 47.3 15.0 B 4 (%) 31.3 70.2 51.0 58.4 C 5 (%) 42.6 5.9 1.7 26.1 D 6 (%) 18.9 0.0 0.0 0.5 Site 3 Depth 10 cm 20 cm 30 cm 40 cm NSE 1.66( 5.39 -0.3) 0.89(0.78 0.94) 0.88(0.79 0.93) 0.65(0.24 0.91) RMSE 0.70(0.70 1.30) 0.41(0.31 0.55) 0.38(0.25 0.54) 0.81(0.43 1.13) A (%) 0.0 41.2 37.2 4.6 B (%) 0.0 56.2 61.1 13.6 C (%) 0.0 2.6 1.7 33.4 D (%) 100.0 0.0 0.0 48.4 Site 4 Depth 10 cm 20 cm 30 cm 40 cm NSE 0.81(0.59 0.89) 0. 76 (0. 02 0. 9 4 ) 0.70(0.25 0.91) 0.30( 0.31 0.88) RMSE 0.71(0.53 0.93) 0 37 (0. 20 0. 60 ) 0.32(0.21 0.42) 0.34(0.13 0.55) A (%) 3.3 20.5 3.4 3.2 B (%) 54.7 23.6 19.1 5.3 C (%) 37.7 28.6 37.9 10.4 D (%) 4.3 27.3 39.6 81.1 Site 6 Depth 10 cm 20 cm 30 cm 40 cm NSE 0.78(0.49 0.93) 0. 87 (0. 39 0.97) 0. 75 (0. 52 0. 87 ) 0.20( 1.95 0 .54) RMSE 0.39(0.22 0.54) 0. 22 (0. 13 0. 37 ) 0. 24 (0. 14 0. 33 ) 0.30(0.23 0.35) A (%) 8.5 40 9 1.9 0.0 B (%) 34.6 35.0 25.6 0.0 C (%) 43.8 17.1 59.7 0.5 D (%) 13.1 0 7 12.8 99.5 1 Nash Sutcliffe coefficient 2 Root mean square error 3 A probability of fit being very good 0.9
PAGE 195

195 Table 5 1 3 Comparison of volumetric soil water content before and after a 6, 9 an d 12 cm incremental raise in canal C 111stage ob tained from a Two sample assuming equal variance t 2013). Depth VWC 1 Before VWC after 6cm VWC after 9cm VWC after 12 cm Variance 2 Variance 3 Variance 4 Site 3 10 cm 23.3 23. 6 23.7 23. 9 2.0 2 0 1. 9 20 cm 24.5 24.7 24.8 24.9 10 1 0 0. 9 30 cm 25.5 25.7 25. 8 25. 9 0. 5 0. 5 0. 5 40 cm 25.8 26.0 26.1 26.2 0. 5 0. 5 0. 5 Site 4 10 cm 29. 3 29. 9 30. 2 30.5 1.2 1. 3 1. 4 20 cm 31. 2 31.8 32.2 32.6 1. 3 1. 4 1. 5 30 cm 33.1 33. 8 34.0 34. 2 1. 4 1. 3 1. 3 40 cm 34.1 34. 5 34.6 34. 7 1.1 1. 0 0. 9 Site 6 10 cm 26. 2 26. 6 26. 8 2 7 0 0.8 0. 8 0. 8 20 cm 27. 4 27. 8 2 8 0 28.2 0.5 0. 6 0. 6 30 cm 28. 4 28.8 29.0 29.1 0. 7 0. 6 0. 6 40 cm 29. 1 29. 3 29. 3 29.4 0.5 0. 4 0. 5 1 Volumetric soil water content (m 3 /m 3 ) 2 Pooled variance afte r a 6 cm increase in canal stage 3 Pooled variance after a 9 cm increase in canal stage 4 Pooled variance after a12 cm increase in canal stage

PAGE 196

196 Figure 5 1. Study area showing soil water monitoring sites, agricultural lands adjacent to Everglades Nationa l Park, and canal network within the C 111 basin of south Miami Dade County, Florida.

PAGE 197

197 Figure 5 2 Depiction of the discretizing of the soil profile and location of the EnviroScans

PAGE 198

198 Figure 5 3. Sobol indices on the vertical axis and parameters (tr1 and tr2 are residual soil water content, ts1 and ts2 are saturated soil content, a1 and a2 are inverse of air entry value, n1 and n2 are curve shape parameter, lam is pore connectivity parameter, K2 is saturated hydraulic conductivity and 1 and 2 refer to the soil and limestone bedrock layers) for the WAVE model on the horizontal axis as applied to simulate volumetric soil water content at 10, 20, 30 and 40 cm depths at site 1.

PAGE 199

199 Figure 5 4. Sobol indices on the vertical axis and parameters (tr1 and tr2 are residual soil water content, ts1 and ts2 are saturated soil content, a1 and a2 are inverse of air entry value, n1 and n2 are curve shape parameter, lam is pore connectivity parameter, K2 is saturated hydraulic conductivity and 1 and 2 refer to the soil and limestone layers) for the WAVE model on the horizontal axis as applied to simulate volumetric soil water content at four monitoring depth 10, 20, 30 and 40 cm at site 2.

PAGE 200

200 Figure 5 5. Sobol indices on the vertical axis and parameters (tr1 and tr2 ar e residual soil water content, ts1 and ts2 are saturated soil content, a1 and a2 are inverse of air entry value, n1 and n2 are curve shape parameter, lam is pore connectivity parameter, K2 is saturated hydraulic conductivity and 1 and 2 refer to the soil a nd limestone layers) for the WAVE model on the horizontal axis as applied to simulate volumetric soil water content at four monitoring depth 20, 30 and 40 cm at site3.

PAGE 201

201 Figure 5 6 Sobol indices on the vertical axis and parameters (tr1 and tr2 are resi dual soil water content, ts1 and ts2 are saturated soil content, a1 and a2 are inverse of air entry value, n1 and n2 are curve shape parameter, lam is pore connectivity parameter, K2 is saturated hydraulic conductivity and 1 and 2 refer to the soil and lim estone layers) for the WAVE model on the horizontal axis as applied to simulate volumetric soil water content at four monitoring depth 10, 20, 30 and 40 cm at site 4.

PAGE 202

202 Figure 5 7 Sobol indices on the vertical axis and parameters (tr1 and tr2 are resid ual soil water content, ts1 and ts2 are saturated soil content, a1 and a2 are inverse of air entry value, n1 and n2 are curve shape parameter, lam is pore connectivity parameter, K2 is saturated hydraulic conductivity and 1 and 2 refer to the soil and lime stone layers) for the WAVE model on the horizontal axis as applied to simulate volumetric soil water content at four monitoring depth 10, 20, 30 and 40 cm at site 6.

PAGE 203

203 Figure 5 8 Visu al comparison of WAVE simulated and EnviroScan measured volumetric s oil water content at site 1 w h ere the vertical line separate s calibration and validation data sets.

PAGE 204

204 Figure 5 9 Visu al comparison of WAVE simulated and EnviroScan measured volumetric soil water content at site 2 w h ere the vertical line separate s calibra tion and validation data sets.

PAGE 205

205 Figure 5 10 Visu al comparison of WAVE simulated and EnviroScan measured volumetric soil water content at site 3 w h ere the vertical line separate s calibration and validation data sets.

PAGE 206

206 Figure 5 11 Visu al comparison o f WAVE simulated and EnviroScan measured volumetric soil water content at site 4 w h ere the vertical line separate s calibration and validation data sets.

PAGE 207

207 Figure 5 12 Visu al comparison of WAVE simulated and EnviroScan measured volumetric soil water conte nt at site 6 w h ere the vertical line separate s calibration and validation data sets.

PAGE 208

208 A B Figure 5 13 Example FITEVAL output for WAVE model validation at site 1 at 30 cm depth from 01 21 2011 to 02 28 2013: A) Goodness of fit not considering measuremen t uncertainty and B) Goodness of fit considering measurement uncertainty.

PAGE 209

209 Figure 5 14 Simulated volumetric soil water content under different C 111 canal stage management scenarios at four different depth at site 3 adjacent to the reach of C 111 betw een the gated spillways at S177 and S18C expected to experience raises in canal stage.

PAGE 210

210 Figure 5 15 Simulated volumetric soil water content under different C 111 canal stage management scenarios at four different depth at site 4 adjacent to the reach of C 111 between the gated spillways at S177 and S18C expected to experience raises in canal stage.

PAGE 211

211 Figure 5 16 Simulated volumetric soil water content under different C 111 canal stage management scenarios at four different depth at site 6 adjacent to the reach of C 111 between the gated spillways at S177 and S18C expected to experience raises in canal stage.

PAGE 212

212 CHAPTER 6 CONCLUSIONS AND FUTURE RESEARCH Chapter Summaries Hydrological modifications are taking place in the C 111 basin of south Florida as part of the Comprehensive Everglades Restoration Plan (CERP). These modifications are being implemented under the C 111 spreader canal project a component of CERP. C 111 spreader canal project aims to restore the hydrology of Everglades National Park ( ENP) in part by reducing groundwater seepage from Taylor Slough and ENP into C 111. To achieve this, structural and operational modifications involving incremental raises in canal stage are planned along C 111. This study used modeling to assess the impact s of the proposed operational adjustments in surface water management on ground and soil water in agricultural lands east of C 111. The agricultural lands east of C 111 are used for production of winter vegetables which is a significant contributor to bot h Miami Dade County and state of Florida economies. Therefore, C 111 spreader canal project is also tasked with ensuring continued flood protection for farmlands east of C 111. The overall contribution of this study was the generation of information throug h monitoring and development of C 111 basin specific modeling tools that can be applied to investigate how the proposed raises in canal stage along C 111 could impact agricultural production through raised water table elevation and saturated root zone. The objectives of the study were achieved through completion of four independent investigations: 1) Dynamic Factor Analysis (DFA) was used to evaluate the impact of surface water management on soil and limestone water content, 2) hydraulic parameters governi ng canal aquifer interaction were estimated using an approximate analytical model of groundwater flow and

PAGE 213

213 sensitivity analysis and parameter estimation techniques, 3) water table response to the proposed incremental raises in canal stage as well as other l arge storm events were simulated using a MODFLOW based groundwater flow model developed for the study area and 4) soil and limestone bedrock water content dynamics in response to proposed changes in canal stage were mechanistically simulated using WAVE ba sed vadose zone models. DFA modeling revealed that canal stage and recharge were the strongest drivers of change in soil water content in the root zone. These findings confirm the perception that raising canal stage under the C 111 spreader canal project c ould affect soil water levels in the root zone of farmlands east of canal C 111 Since net recharge cannot be manipulated, maintaining a favorable balance of air and water in the unsaturated zone would require future operational management of canal stage l evels to minimize entry of the water table into the root zone. Groundwater flow simulations in this study predicted that the impact of raising canal stage will depend on micro topography within the field and not on the distance from the canal. Low elevati on sites (<2.0 m NGVD29) are predicted to experience pro longed saturation of the root zone resulting in shortening of the growing season for raises in canal stage greater 9 cm To minimize impact of raising canal stage on agricultural production during th e growing season (September to April), while ensuring continued environmental benefits of reduced groundwater seepage from Taylor Slough canal stage could be raised by not more than 9 cm during the growing season and l ater in the wet season could be raised to the proposed maximum of 12 cm. The study

PAGE 214

214 findings also indicated that to minimize flooding from forecasted large storms, canal stage should be drawn down at least 48 hours before large storms. In addition to canal drawdown, to minimize losses to the gr owers associated with plant death from saturated root zones, vulnerability maps showing which parts of the landscape w ould be impacted by the different incremental raises in canal stage would provide a useful guide. These maps could be generated by integr ating high resolution DEM ( digital elevation model ) with the groundwater flow model developed for the study area I also suggest that these vulnerability maps in future be packaged in extension format to facilitate transfer of information. C hapter 2 DFA w as used as an alternative tool to physically based models to explore the relationship between different hydrologic variables and the effect of proposed changes in surface water management on soil and bedrock water contents in south Florida. T he specific ob jectives were to: (1) use DFA to identify the most important factors affecting temporal variation in soil and bedrock water contents, (2) develop a simplified DFA based regression model for predicting soil and bedrock water contents as a function of canal stage and (3) assess the effect of the proposed incremental raises in canal stage on soil and bedrock water contents. DFA revealed that 5 common trends were the minimum required to describe unexplained variation in the 11 time series studied. Introducing canal stage, water table evaporation and net recharge resulted in lower Akaike information criterion (AIC) and higher Nash Sutcliffe (Ceff) values. Results indicated that canal stage significantly (t > 2) drives temporal variation in soil and bedrock wate r contents, which was represented as scaled frequency while net surface recharge was significant in 7 out of the 11 time

PAGE 215

215 series analyzed. The effect of water table evaporation was not significant at all sites. Results also indicated that the most important factor influencing temporal variation in soil and bedrock water contents in terms of regression coefficient magnitude was canal stage. Based on DFA results, a simple regression model was developed to predict soil and bedrock water contents at various elev ations as a function of canal stage and net recharge. The performance of the simple model ranged from good (Ceff ranging from 0.56 to 0.74) to poor (Ceff ranging from 0.10 to 0.15), performance was better at sites with smaller depths to water table (< 1 m) highlighting the effect of micro topography on soil and bedrock water content dynamics. Assessment of the effect of 6, 9 and 12 cm increases in canal stage using the simple regression model indicated that changes in temporal variation in soil and bedroc k water contents were negligible (average<1.0% average change) at 500 to 2000 m from C111 (or low elevations) which may be attributed to the near saturation conditions already occurring at these sites. This study used DFA to explore the relationship betwee n soil and bedrock water dynamics and surface water stage in shallow water table environments. This approach can be applied to any system in which detailed physical modeling would be limited by inadequate information on parameters or processes governing th e physical system. However, model application is limited by the dataset used for its development due to its empirical nature and thus there is uncertainty in its ability to be used as a prediction tool. C hapter 3 The goal of this study was to better charac terize parameters influencing the 111 and Biscayne aquifer using the analytical model STWT1. A three step model evaluation framework was

PAGE 216

216 implemented as follows: 1) qualitative parameter ranking by compar ing two Morris method sampling strategies, 2) quantitative variance based sensitivity analysis using statistics using the Generalized Likelihood Uncertainty Estimator (G LUE) methodology. Results indicated that the original Morris random sampling method under estimated total parameter effects compared to the improved global Morris sampling strategy. However, parameter ranking from the two sampling methods was similar. For the STWT1 model, only four out of the six parameters analyzed were important for predicting water table response to canal stage and recharge fluctuations. Morris ranking in order of decreasing importance resulted in specific yield (ASY), aquifer saturated thickness (AB), horizontal hydraulic conductivity (AKX), canal leakance (XAA), vertical hydraulic conductivity (AKZ), and half width of canal (XZERO). summation of first order parameter effects was 1.0 indicating that STWT1 behaved as an additive model or negligible parameter interactions. We estimated parameter values of 0.07 to 0.14 for ASY, 11,000 to 14,300 m/day for AKX, 13.4 to 18.3 m for AB, and 99.8 to 279 m for XAA. The estimated values were within the range of values estimated using more complex methods at nearby sites. Nash Sutcliffe and root mean square error using estimated parameters to predict water table elevation ranged from 0.66 to 0.95 and from 4 to 7 cm, respe ctively. This study demonstrates a simple and low cost way to characterize hydrogeological parameters controlling groundwater surface interactions in any region with aquifers that are highly permeable to be efficiently characterized using standard pumping tests or canal drawdown experiments.

PAGE 217

217 Hydrogeological parameters estimated using this approach could be used as starting values in large scale numerical simulations. C hapter 4 The specific objectives were to: (1) develop a MODFLOW based groundwater flow mod el, (2) apply the developed model to determine if the proposed changes in canal stage result in significant changes in water table elevation, root zone saturation or groundwater flooding and (3) assess aquifer response to large rainfall events. Results ind icate the developed model was able to reproduce measured water table elevation, with average Nash Sutcliffe > 0.9 and Root Mean Square Error <0.05 m. The model predicted that incremental raises in canal stage resulted in significant differences (p<0.05) in water table elevation. Increases in canal stage of 9 and 12 cm resulted in occasional root zone saturation of low elevation sites. The model was able to mimic the rise and fall of the water table pre and post Tropical Storm Isaac of August 2012. The mode l also predicted that lowering canal stage at least 48 hours prior to large storm (>2 return period storm) reduced water table intrusion into the root zone. We conclude d that the impact of operational changes in canal stage management on root zone saturati on and groundwater flooding depended on micro topography within the field and depth of storm events using the MODFLOW based model. The findings of this study can be used in fine tuning canal stage operations to minimize root zone saturation and groundwater flooding of agricultural fields while maximizing environmental benefits in ENP. This study also highlights the benefit of detailed field scale simulations.

PAGE 218

218 C hapter 5 The specific o bjectives of the study were to: 1) develop one dimensional unsaturated flow models for simulating soil water dynamics at specific sites within the farmlands east of canal C 111 model s for simulating soil and limestone bedrock water dynamics at the study area, and (2) apply the model s to investigate the impact of the proposed incr emental rises in canal C 111 stage of up to 12 cm on root zone soil and bedrock water content in adjacent farm lands. WAVE a generic vadose zone computer model was applied in simulating soil and bedrock water dynamics. Model performances were evaluated with and without consideration of uncertainty in measured water content data. The validated model s for our study area w ere then applied to explore potential impacts of the proposed increases in canal stage on root zone water content. Global sensitivity analys is was performed to identify model input factors with greatest influence on predicted soil and bedrock water content. Qualitative parameter ranking using the improved Morris method revealed that predicted soil and bedrock water content were most sensitive to variations in the parameters of the van Genuchten equation which was used as the water retention model in WAVE. Quantitative variance content had the greatest influence on predicted soil and bedrock water content. Sensitivity analysis also revealed that WAVE behaved as an additive model in the bottom 20 cm monitored. The model was manually calibrated by adjusting parameters within ranges measured in the laboratory or obta ined in literature until an acceptable fit between predicted and measured soil and bedrock water content was graphically achieved. To avoid over fitting of model parameters to uncertain measured water content data automated calibration using the GLUE metho dology was abandoned.

PAGE 219

219 Model performance was greatly improved when uncertainty in measured data was taken into consideration during model performance evaluation. There were some variations in model performance at different sites and monitoring depth. These differences were attributed to: variations in thickness of the cement slurry surrounding the access tube, intrinsic spatial variability in soil and limestone hydraulic properties, effect of not including irrigation water inputs in the conceptual models, a nd uncertainties in crop coefficient time series. Overall model performance ranged from good to acceptable at four out of the five sites analyzed. We concluded that WAVE based models developed for the study area with the exception of the model developed fo r site 2 were useful as exploratory tools for estimating average soil water content conditions in the top 40 cm depth and for predicting temporal variations in soil water content in response to surface wate r management in the C 111 canal Proposed incremen tal raises in canal stage were mimicked by incrementally adding 6, 9, and 12 cm to the current canal C 111 stage up to a total of 12 cm. Soil and bedrock water content before and after the incremental raises in canal stage were significantly different (p m NGV29 experienced root zone saturation in addition to shortening of the growing period. Site with surface elevation greater than 2.0 NGVD29 did not experience root zone saturation even after the maximum incremental raise in stage of 12 but experienced prolonged saturation at 30 and 40 cm which could exacerbate the problem of temporary flooding associated with spikes in water table elevation after storms. The results indicate that micro topogr aphy within the field will be the main factor influencing which areas would be impacted by the propo sed incremental raises in canal

PAGE 220

220 stage. The model s developed in this study could be used with high resolution digital elevation model (DEM) to identify areas that should not be planted to minimize potential losses associated with saturated root zones. However, application of WAVE models should be done considering the measurement uncertainty on which they were developed. Although the goodness of fit values were acceptable, the measurement uncertainty was still high. Model improvement would depend on decreasing measurement uncertainty. This study shows that in humid water table controlled landscapes, surface water management could impact soil water content which in turn influences land use and other near land surface water and energy fluxes. Suggested Future Research Several limitations were encountered in conducting these studies including: 1) difficult in proper installation and calibrating of EnviroScans for soil and bedrock water measurement in thin gravely loam soils and shallow limestone bedrock, 2) inadequate and uncertainty on phonological and growth characteristics of the selected crop due in part to lack of a controlled crop experiment and 3) limitatio n of available vadose zone models to integrate sensitivity and parameter estimation framework modules within the model structure as well as the limitations of Richards under our study conditions. Therefore, future studies should attempt to resolv e some of these problems. Soil Water Measurement T echnologies Other investigators (Evett et al., 2006; Rowland et al., 2011) have noted that technologies that measure soil water content by sensing electromagnetic properties of the soil such as EnviroScans are prone to uncertain precision and sensed volumes issues (i.e., measure a very small volume of soil) C onditions at our study site exacerbated these problems because installation of the access tube required drilling a

PAGE 221

221 hole in the limestone bedrock and f illing the annular space with cement slurry. This introduced another material i.e., cement slurry to become part of the sensor soil matrix which introduces errors. The challenge is to know what viscosity of the slurry to use and how to maintain a very thin cement lining around the tube that minimizes the problem of air pockets but does interfere with soil water measurement. Future investigation s that include measurement of soil water content under such conditions similar to our study site should : 1) perform field evaluat ion of alternative soil water measurement technologies e.g., the TS1 smart tensiometers that over come problems of traditional tensiometers such as cavitation at pressures below 80 kPa, the TS1 can be used for continuous logging of water pote ntial and temperature data and are designed to be self refilling 2) generate field based soil water retention curves by combining TS1 with gravimetric measurements of soil water content and 3) future studies could also investigate the best material (e.g., cement or local soil slurry) for use in making the slurry used to fill the annular space during installation of EnviroScans. Vulnerability maps and the Impact of the Rapid Water Table Raise Phenomenon on Vegetable and Fruit Crop production Rapid water tab le raise phenomenon in the Biscayne aquifer should be further investigated in conjunction with high resolution DEM to identify locations prone to root zone saturation and groundwater flooding. This information could be in form of e.g., spatial maps showing areas prone to groundwater flooding This type of information could help growers minimize losses by not planting in such areas that would be affected

PAGE 222

222 Data on Crop Phonological and Growth Characteristics To properly characterize the impacts of surface w ater management on crop growth future studies need to focus on obtaining data that characterizes the phonological and growth characteristics of the crop simulated in vadose zone models such as WAVE. I suggest having a field experiment in which the followi ng data are measured or estimated : 1) crop coefficient time series, 2) leaf area index and 3) root growth. Having these data would reduce uncertainty in model predictions of unsaturated flow. Measuring or estimating these data in the field would enable eas y interpretation compared to conducting a soil column based laboratory experiment that evaluates the impact of depth to water table on plant growth. Such a controlled study would also allow measuring other water balance inputs such as irrigation which be d ifficult to measure on a private field irrigation water is not metered Coupling Vadose Zone Models to Sensitivity and Parameter Estimation Tools To quicken evaluation of vadose zone models future research efforts could also be target ed at developing inte rfaces for coupling processes based vadose zone models such as WAVE to public domain sensitivity analysis tools (e.g., SIMLAB, UCODE [ universal code for inverse modeling ] ) and parameter estimation tools (e.g., PEST [ simulation laboratory for UA/SA ] UCODE, GLUE [ generalized likelihood uncertainty engine ] SCE [ shuffled complex evolution ] ) and multi modeling analysis tools (e.g., MMA [ multimodel analysis ] ) This would substantially reduce the time spent on developing model interface computer programs each tim e one requires use a model evaluation tool The time freed by integrating model eva lua tion tools into process based models would be spent in developing more representative models of physical systems

PAGE 223

223 This approach has been attempted in proprietary vadose z one models such as HYDRUS 2D/3D but most frameworks do not yet incorporate global approaches In addition, to accounting for measurement uncertainty during model evaluation, future studies could also focus on developing methodologies for evaluating model uncertainty. Some methodologies exist for quantifying model uncertainty however; they make the assumption of predicted values being independent which might not be valid for unsaturated soil water content predictions due to autocorrelation. Climatic impacts Climate variability and change is evident therefore future studies could also focus on evaluating climate change and variability impact s on surface water in canals, groundwater table levels and soil water content For example future studies could explor e linking fully coupled models of surface water groundwater vadose zone interaction to climate models.

PAGE 224

224 LIST OF REFERENCES Akaike, H. 1974. A new look at the statistical model identification. IEEE Trans. Automat. Control 19 : 716 723. Allen, R. 2011. REF E T: Reference Evapotranspiration Calculation Software. User Manual. Version: 3.1.08 Al Yahyai R B. Schaffer, F.S. Davies, R. Muoz Carpena 2006. Characterization of soil water retention of a very gravelly loam soil varied with determination method. Soi l Science 171(2): 85 93. DOI:10.1097/01.ss.0000187372.53896.9d. Anderson, M.P., W.W. Woessner.1992. Applied Groundwater Modeling Academic Press, San Diego: 381 pp. 28. ASCE 2005. The ASCE Standardized Reference evapotranspiration Equation. Task Committ ee on Standardization of Calculation of Reference ET. Environment and Water Resources Institute of ASCE. 200 p. Barlow, P.M., and A.F. Moench 1998. Analytical solutions and computer programs for hydraulic interaction of stream aquifer systems, US Geologic al Survey Open File Report 98 415A, pp. 85. Barquin, L.P., K.W. Migliaccio, R. Muoz Carpena B. Schaffer J.H. Crane. Y.C. Li. 2011. Shallow Water Table Contribution to Soil Water Retention in Capillary Fringe of a Very Gravelly Loam Soil of South Florida Vadose Zone J 10:1 8. Bear, J. 1972. Dynamics of Fluids in Porous Media. American Elsevier Pub. Co. Beven, K. 2006. A manifesto for the equifinality thesis. Journal of Hydrology 320:18 36. Beven, K. J., and A.M. Binley. 1992. The future of distribut ed models: Model calibration and uncertainty prediction. Hydrol. Processes 6:279 298. Bolster, C., D. Genereux, J. Saiers 2001. Determination of specific yield for a limestone aquifer from a canal drawdown test. Ground Water. 39:768 777. Boussinesq, M.J. 1877. Essai sur la theorie des eaux courantes. Memoires de 23(1): 252 260. Brion, L., S. U. S. Senarat, A. M. W. Lal, and M. Belnap. 2001. Application of the South Florida Regional Simulation Model in t he Everglades. Proc., ASCE Specialty Symp.: Integrated Surface and Groundwater Management, Orlando, Fla. Brooks, R. H., and A. T. Corey. 1964. Hydraulic properties of porous medium. Hydrology paper No.3, Civ. Engrg. Dept.. Colorado State Univ., Port Collin s, Colo.

PAGE 225

225 Buckingham, E. 1907. Studies on the movement of soil moisture. Department of Agriculture Bureau of Soils Washington, DC. Bulletin 38, pp 28 61. Campolongo. F., J. Cariboni, and A. Saltelli. 2007. An effective screening design for sensitivity analy sis of large models. Environmental Modeling and Software 22:1509 1518. Celia, M.A., E.T. Bouloutas, and R.L. Zarba.1990. A general mass concervative numerical solution for the unsaturated flow equation Water Resour. Res. 26(7):1483 1496. Chan, K., A. Sal telli, and S. Tarantola. 1997. Sensitivity Analysis of Model Output: Variance Based Methods Make the Difference. In Proceedings of the 1997 Winter Simulation Conference. ed. S. Andradttir, K. J. Healy, D. H. Withers, and B. L. Nelson. Chen, M., G. R. Will goose, and P. M. Saco. 2012. Spatial prediction of temporal soil moisture dynamics using HYDRUS 1D. Hydrol. Process DOI: 10.1002/hyp.9518. Chin, D. 1990. A method to estimate canal leakage into the Biscayne Aquifer, Dade County, Florida. Water Resour. Inv estigations Rep. No. 90 4135, U.S. GeoI. Survey, Tallahassee, Fla. Chin, D. 1991. Leakage of clogged channels that partially penetrate surficial aquifers. ASCE J. Hydraulic Engineering 117 : 467 488. Chin, D. 2008. Phenomenological models of hydrologic pro cesses in south Florida. J. Hydrol 349 : 230 243. Cooper Jr, H.H., and M.I Rorabaugh. 1963. Ground water movements and bank storage due to flood stages in surface streams, USGS Water Supply Paper 1536:343 366. DACS, 2012. DACS P 01304 Florida Agriculture b y the Numbers 2012 Statistical Directory Available online at: http://florida agriculture.com/brochures/P 01304 .pdf Accessed April 28, 2013. Dean, T.J., J.P. Bell, and A.J.B. Baty 1 987. Soil moisture measurement by an improved capacitance technique. Part 1: sensor design and performance. Journal of Hydrology 93:67. Delleur, J.W. 1999. Handbook of Groundwater Engineering 2nd addition., (ed). CRC Press Taylor & Francis Group Boca Rat on, FL. Dempster, A.P., N.M. Laird and D.B. Rubin 1977. Maximum likelihood from incomplete data via the EM algorithm. J.R. Stat. Soc. Ser. 39 : 1 38.

PAGE 226

226 Duwig, C., B. Normand, M. Vauclin, G. Vachaud, S.R. Green, T. Becquer 2003. Evaluation of the WAVE model for predicting nitrate leaching for two contrasted soil and climate conditions. Vadose Zone J 2 : 76 89. Evett., S. R., R. C. Schwartz, J. A. Tolk, T. A., Howell. 2006. Soil Profile Water Content Determination: Spatiotemporal Variability of Electromagnetic and Neutron Probe Sensors in Access Tubes. Vadose Zone J 8:926 94. Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulation of field water use and crop yield. Simulation Monographs. PUDOC, Wageningen, The Netherlands. Fish, J. E., and M. Stewart. 199 1. Hydrogeology of the surficial aquifer system, Dade County, Florida. Water Resource. Investigations Rep. No. 190 4108, U.S. GeoI. Survey, Tallahassee, Fla. Freeze, R.A., and J.A. Cherry. 1979. Groundwater Prentice Hall, Englewood Cliffs, N. J. Frey, H.C ., and S.R. Patil. 2002. Identification and Review of Sensitivity Analysis Methods. Risk Analysis 22(3): 553 578. Gabriel, J. L., J.I. Lizaso Q. Miguel 2010. Laboratory versus Field Calibration of Capacitance Probes. Soil Sci Soc Am J 74 : 593 601. Gardn er, W. R. 1958. Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Sci. 8.5(4): 228 232. Genereux, D., and E. Slater, 1999. Water exchange between canals and surrounding aquifer an d wetlands in the Southern Everglades, USA. Journal of Hydrology 219:153 168. Genereux, D.P., and J.D. Guardiario. 1998. A canal drawdown experiment for determination of aquifer parameters. J. Hydr. Eng 3 (4): 294 302. Germann, P. and B. Levy 1986. Gr oundwater response to precipitation. EOS, 67: 91. Gillham, R.W. 1984. The capillary fringe and its effect on water table response. J. Hydrol ., 67: 307 324. Geweke, J.F., 1977. The dynamic factor analysis of economic time series models. In: Aigner, D.J., Go ldberger, A.S. (Eds.), Latent Variables in Socio economic Models. North Holland, Amsterdam, pp. 365 382. Graham, W.D., K.L. Campbell, J. Mossa, L.H. Motz, P.S.C. Rao, W.R. Wise,. and D.P. Genereux .1997. Water management issues affecting the C 111 Basin, Dade County, Florida. Report number CNR 1997 1002, Center for Natural Resources, University of Florida, Gainesville, 156 pp.

PAGE 227

227 Ha, K., D. Koh, B. Yum, and K. Lee. 2007. Estimation of layered aquifer diffusivity and river resistance using flood wave response model. Journal of Hydrology 337(3 4): 284 293. Hall, F.R., and A.F. Moench.1972. Application of the convolution equation to stream aquifer relationships. Water Resources Research 8 (2): 487 493. Hantush, M.M. 2005. Modeling stream aquifer interactions wi th linear response functions. Journal of Hydrology 31:59 79. Hantush, M.S. 1965. Wells near streams with semi pervious beds. J. Geophys. Res 70(12): 2829 2838. Harbaugh, A.W., E.R. Banta, M.C. Hill, and M.G. McDonald 2000, MODFLOW 2000, the U.S. Geolog ical Survey modular ground water model -User guide to modularization concepts and the Ground Water Flow Process: U.S. Geological Survey Open File Report 00 92, 121 p. Harvey, A.C. 1989. Forecasting, structural time series models and the Kalman filter Ca mbridge Univ. Press, New York. Harmel, R. D., P. K. Smith, and K. L. Migliaccio. 2010. Modifying goodness of fit indicators to incorporate both measurement and model uncertainty in model calibration and validation. Trans. ASABE 53(1): 55 63. Harmel, R. D., and P. K. Smith. 2007. Consideration of measurement uncertainty in the evaluation of goodness of fit in hydrologic and water quality modeling. J. Hydrol. 337(3 4): 326 336. Heliotis, F. D ., and C. B. DeWitt. 1987. Rapid Water Table Responses To Rainfall In A Northern Peatland Ecosystem. JAWRA Journal of the American Water Resources Association 23 : 1011 1016. Hill, M.C. 1998. Methods and guidelines for effective model calibration: U.S.Geological Survey, Water Resources Investigations Report 98 4005, 90 p. Hill, M.C., E.R. Banta, A.W. Harbaugh, and E.R. Anderman 2000. MODFLOW 2000, The U.S. Geologi cal Survey Modular Ground Water Model -Observation, Sensitivity, and Parameter Estimation Processes and Three Post Processing Programs: U.S. Geological Survey Open File Report 00 184, 220 p. Hill, M.C. 1998. Methods and guidelines fo r effective model calibration: U.S.Geological Survey, Water Resources Investigations Report 98 4005, 90 p. Hughes, J.D., C.D. Langevin, K.L. Chartier, and J.T. White 2012. Documentation of the Surface Water Routing (SWR1) Process for Modeling Surface Wat er Flow with the U.S. Geological Survey Modular Groundwater Model (MODFLOW 2005), U.S. Geological Survey Techniques and Methods 6 A40.

PAGE 228

228 Huyakorn, P.S., Pinder, G.F., 1983. Computational Methods in Subsurface flow New York: Academic Press. IAEA, 2008. Field Estimation of soil water content; A practical guide to methods, instrumentation and sensor technology. Training Course Series No. 30, International Atomic Energy Agency, Vienna, Austria. http://www pub.iaea.org /MTC D/publications/PDF/TCS 30_web.pdf Accessed 12 July 2011. Kaplan, D., and R. Muoz Carpena 2011. Complementary effects of surface water and groundwater on soil moisture dynamics in a degraded coastal floodplain forest. J. of Hydrol 398, 221 234. Kaplan, D., R. Muoz Carpena, and A. Ritter 2010a Untangling complex groundwater dynamics in the floodplain wetlands of a southeastern U.S. coastal river. Water Resour. Res. 46, W08528 10 doi:10.1029/2009WR009038. Kaplan, D., R. Muoz Carpena, Y. Wan, M. Hedgepe th, F. Zheng, R. Roberts, and R. Rossmanith 2010b Linking river, floodplain, and vadose zone hydrology to improve restoration of a coastal river impacted by saltwater intrusion. J. Environ. Quality 39 (5), 1570 1584. doi:10.2134/jeq2009.0375. Kayane I., and I. Kaihotsu. 1988. Some Experimental Results Concerning Rapid Water Table Response to Surface Phenomena. Journal of Hydrology 102: 215 234. Kisekka, I., K.W. Migliaccio, R. Muoz Carpena, Y. Khare, and T. H. Boyer. 2013b. Sensitivity analysis and pa rameter estimation for an approximate analytical model of canal aquifer interaction applied in the C 111 Basin. In revision Trans. ASABE SW 10037 2012. Kisekka, I., K.W. Migliaccio, R. Muoz Carpena,and B. Schaffer. 2013 a Dynamic Factor Analysis of Surfa ce Water Management Impacts on Soil and Bedrock WaterContents in Southern Florida Lowlands. Journal of Hydrology 55 72. http://dx.doi.org/10.1016/j.jhydrol.2013.02.035 Konikow L.F. 1996. Us e of Numerical Models to Simulate Groundwater Flow and Transport, US Geological Survey, Reston, Virginia, USA. Konikow L.F., and T.E. Reilly. 1998. Groundwater Modelling In: The Handbook of Groundwater Engineering [J.W. Delleur, ed.], CRC Press, Boca Rato n 20:1 20.40. Krause, P., D.P. Boyle, and F. Base. 2005. Comparison of different efficiency criteria for hydrological model assessment. Advances in Geosciences 5:89 97. Kumar, C.P. 2002. Modelling of unsaturated flow, National Conference on "Modern Trend s in Water Resources Development and Environmental Management", Vellore nstitute of Technology, Vellore (Tamil Nadu), str. 1 9.

PAGE 229

229 Lal, A.M.W. 2001. Modification of canal flow due to stream aquifer interaction. Journal of Hydraulic Engineering 127(7): 567 57 6. Lal, A.M.W. 2006. Determination of multiple aquifer parameters using generated water level disturbances, Water Resour. Res doi:10.1029/2005WR004218. Ltkepohl, H. 1991. Introduction to multiple time series analysis Springer Verlag, Berlin. Mrkus, L. O. Berke, J. Kovcs, and W. Urfer 1999. Spatial prediction of the intensity of latent effects governing hydrogeological phenomena. Environmetrics 10 : 633 654. McDonald, M.G., and A.W Harbaugh 1988. A Modular Three Dimensional Finite Difference Ground Water Flow Model. Techniques of Water Res. Invests. of the U.S. Geol. Survey, Book 6, Ch. A1: 586 pp. Merdun, Hasan network and regression pedotransfer functions for prediction of so il water retention and saturated hydraulic conductivity. Soil and Tillage Research 90 (1 2): 108 116. Merkel, R. 2000. Element and Sediment Accumulation Rates in the Florida Everglades. Water, Air, and Soil Pollution 122 (3 4): 327 349. Merritt, M.L. 19 96. Simulation of the water table altitude in the Biscayne aquifer, Southern Dade County, Florida, water years1945 89. Water Supply Paper 2458, US Geological Survey. Moench A.F., and P.M. Barlow 2000. Aquifer response to stream stage and recharge varia tions: I. Analytical step response functions. Journal of Hydrology 230:192 210. Morris, M.D. 1991. Factorial Sampling Plans for Preliminary Computational Experiments Technometrics 33:161 174. Mualem, Y. 1976. A new model for predicting the hydraulic cond uctivity of unsaturated porous media. Water Res. Res 12(3):513 522. Muoz Carpena R, and Y. Li. 2003. Study of the Frog Pond area hydrology and water quality modifications introduced by the C 111 Project detention pond implementation. Final Project Repor t. TREC RMC 2003 01. Available online at http://abe.ufl.edu/faculty/carpena/files/pdf/research/reports/StudyOfTheFrogPond AreaHyd rologyAndWaterQuality.pdf Accessed on 23 January 2012. Muoz Carpena, R., Z. Zajac, and Y.M. Kuo. 2007. Global Sensitivity and Uncertainty Analyses of the Water Quality Model VFSMOD. Transactions of the ASABE 50(5): 1719 1732.

PAGE 230

230 Muoz Carpena, R., 2004 Field devices for monitoring soil water content. University of Florida IFAS edis BUL 343. Muoz Carpena, R., A. Ritter, and Y.C. Li. 2005. Dynamic factor analysis of groundwater quality trends in an agricultural area adjacent to Everglades National Park J. Contam. Hydrol 80: 49 70. Muoz Carpena, R., Ritter, A., Bosch, D.D., Schaffer, B., Potter, T.L., 2008. Summer cover crop impacts on soil percolation and nitrogen leaching from a winter sweet corn field. J. Agricultural Water Management 95 : 633 64 4. Munsell Soil Color Charts. 2000. Revised Edition. Greta G. Macbeth. New Windsor, NY. Nash, J.E., and J.V. Sutcliffe 1970. River flow forecasting through conceptual models: Part 1 A. Discussion of principles J. Hydrol 10 : 282 290. Nobel, C.V., R.R.W Drew, and J.D. Slabaugh 1996. Soil survey of Dade County Area, Florida, U.S. Department of Agriculture, NRCS Report, Washington, DC. Available online at http://soildatamart.n rcs.usda.gov/manuscripts/FL686/0/ Dade.pdf Accessed on 31 May 2011. Olsthoorn, N.T. 2008. Do a Bit More with Convolution. Ground Water 46(1): 3 22. Parlange, J.Y., W.L. Hogarth, R.S. Govindaraju, M.B. Parlange, and D. Lockington. 2000. On an exact analy tical solution of the Boussinesq equation. Transport in Porous Media 39: 339 345. Pathak, S. C. 2001 Frequency Analysis of Short Duration Rainfall for Central and South Florida. ASCE Environmental and Water Resources 2001, Bridging the Gap. doi: 10.1061/4 0569(2001)227 Philip, J. R 1957. The Theory of Infiltration, Soil Sci 83 ( 345 357 ). Ramirez, J.A. and B.D., Finnerty 1996. Precipitation and Water Table Effects on Agricultural Production and Economics. Journal of Irrigation and Drainage Engineering 1996: 164 171. Ratto, M., and A. Saltelli. 2001. Deliverable 18: Model assessment in integrated procedures for environmental impact evaluation: software prototypes. Estimation of human impact in the presence of natural fluctuations EC Fifth Framework Progr am SCA. Project IST 1999 11313. Richards, L.A. 1931. Capillary conduction through Porous Mediums. Physics 1:313 318.

PAGE 231

231 Ritter, A., and R. Muoz Carpena 2006. Dynamic factor modeling of ground and surface water levels in an agricultural area adjacent to Ev erglades National Park. Journal of Hydrology 317: 340 354. Ritter, A., and R. Muoz Carpena 2012. Predictive ability of hydrological models: objective assessment of goodness of fit with statistical significance. J. of Hydrology doi:10.1016/j.jhydrol.201 2.12.004 Ritter, A., C.M. Regalado, and R. Muoz Carpena 2009. Temporal common trends of topsoil water dynamics in a humid subtropical forest watershed. Vadose Zone J 8(2) : 437 449. Rowland, R., Y.A. Pachepsky, and K. A. Guber. 2011. Sensitivity of a Capacitance Sensor to Artificial Macropores. Technical Article. Soil Science 176(1):9 14. doi: 10.1097/SS.0b013e318203722d. Ruano, M.V., J. Ribes, A. Seco, and J. Ferrer. 2012. An improved sampling strategy based on trajectory design for application of th e Morris method to systems with many input factors. Environmental Modeling & Software 37: 103 109. Saiers, E. J., D. P. Genereux, and C. H. Bolster 2004. Influence of Calibration Methodology on Ground Water Flow Predictions. Ground Water 42(1): 32 44. 32 44 Salas, J.D., 1993. Analysis and modeling of hydrologic time series p. 19.1 19.72. In D.R. Maidment (ed.) Handbook of hydrology. McGraw Hill, New York. Saltelli, A., T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Ratto, M. Saisana, and S. Tara ntola. 2008. Global Sensitivity Analysis. The Primer John Wiley & Sons. Saltelli, A., F. Campolongo, and J. Cariboni. 2009. Screening important inputs in models with strong interaction properties. Reliability Engineering and Systems Safety (94):1149 1155 Saltelli, A., K. Chan, and E.M. Scott. 2000. Sensitivity Analysis John Wiley and Sons, Ltd.: West Sussex, England. Saltelli, A., M. Ratto, S. Tarantola, and F. Campolongo. 2005. Sensitivity analysis for chemical models. Chem. Rev 105(7):2811 2827. Salt elli, A., S. Tarantola, F. Campolongo, and M. Ratto. 2004. Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models Chichester, U.K.: John Wiley and Sons. Schaffer, B.,1998. Flooding Responses and Water use Efficiency of Subtropical and Tr opical Fruit Trees in an Environmentally sensitive Wetland. Annals of Botany 81, 475 481.

PAGE 232

232 Schroeder, M.C., H. Klein, and N.D. Hoy.1958. Biscayne aquifer of Dade and Broward Counties, Florida (FGS: Report of investigations 17 p. 55 56). Tallahassee, FL. Ava ilable at: http://ufdc.ufl.edu/UF00001201/00001 Accessed on 16 January 2012. Sentek Pty Ltd. 2001. Calibration of the Sentek Pty Ltd Soil Moisture Sensors. Stepney, South Australia. Soil Survey Divisio n Staff. 1993. Soil survey manual. Serrano, S.E., and S.R. Workman 1998. Modeling transient stream/aquifer interaction with the non linear Boussinesq equation and its analytical solution, Journal of Hydrology 206(3 4) : 245 255. SFWMD, 2011. Past and Proje cted Trends in Climate and Sea Level for South Florida. Hydrologic and Environmental Systems Modeling. 3301 Gun Club Road West Palm Beach, Florida. Available at: http://my.sfwmd.gov/portal/page/portal/ xrepository/sfwmd_repository_pdf/ccireport_publicationversion_14jul11.pdf Accessed on 10 June 2012. Shumway, R.H., and D.S. Stoffer 1982. An approach to time series smoothing and forecasting using the EM algorithm. J. Time Ser. Anal 3 : 253 264. SIMLAB. 2004. Simulation Environment for Uncertainty and Sensitivity Analysis Version 2.2. Joint Research Centre of the European Commission. Brussels, Belgium. Available at: http://ipsc.jrc.ec.europa.eu/index.php?id=756#c2907 Accessed on 08 June 2013. M. Th. van Genuchten, M. 2008. Development and applications of the HYDRUS and STANMOD software packages, and related codes, Vadose Zone J .7 (2), 587 600. Skinner, C., F. Bloetscher, and C.S. Pathak 2008. Comparison of NEXRAD and Rain Gauge Precipitation Measurements in South Florida. J. of Hydrologic Engineering. 14(3), 248 260. Sobol, I.M. 1993. Sensitivity estimates for non linear mathematical models. Mathematical Modeling and Computational Experiment 4(1):407 414. Srivastava, K., S. E. Serrano, and S.R. Workman. 2006. Stochastic modeling of transient stream aquifer interaction with the nonlinear Boussinesq equation. Journal of Hydrology 328:538 547. Stedinger, J. R., R.M.Vogel, S. Lee, and R. Batchelder. 2008. Appraisal of the g eneralized likelihood uncertainty estimation (GLUE) method. Water resources research. W00B06 doi:10.1029/2008WR006822. Stehfest, H. 1970. Numerical inversion of Laplace transforms: Communications of the Association for Computing Machinery (ACM), 13(1):47 49.

PAGE 233

233 Strowd, T. B. 2012. Tropical Storm Isaac. Readiness, Response, Recovery. Water Conditions Summary and After Action Assessment September 13, 2012. SFWMD. West Palm, Beach Florida. Available at: http://www.sfwmd.gov/portal /page/portal/xrepository/sfwmd_repository_pdf/gb_pres_ts_isaac_2012_0913.pdf Accessed on 06 April 2013. Swain, E.D., and E.J. Wexler 1996. A coupled surface water a nd ground water flow model (MODBRANCH) for simulation of stream aquifer interaction: U.S. Geological Survey Techniques of Water Resources Investigations Report 6 A6, 125 p. USACE. 2009. C 111 Spreader Canal Western Project. Civil Works Review Board Briefin g. December 2009. Jacksonville, FL. Available at: http://www.usace.army. mil/Portals/2/docs/civilworks/CWRB/c_111/c_111_Briefing_slides.pdf Accessed on 28 January 2013. USACE and SFWMD. 2011. Comprehensive Everglades Restoration Plan: C 111 spreader canal western project: Final Integrated Project Implementation Report and Environmental Impact Statement. U.S. Army Corps of Engineers Jacksonville, Distr ict. Available at: http://www.usace.army.mil/Portals/2/docs/ civilworks/CWRB/c_111/c_111_Briefing_slides.pdf Accessed on 20 February 2013. USGS. 2000. Ground Water Flooding in Glacial Terrain of Southern Puget Sound, Washington. USGS Fact Sheet 111 00. Se ptember 2000. U.S. Geological Survey, WRD, 1201 Pacific Avenue, Suite 600, Tacoma, WA 98402. Available at: http://wa.water.usgs.gov/projects/pugethazards/PDF/fs111_00.pdf Acce ssed on 15 March 2013 USGS.1999. Florida Everglades. Circular 1182. U.S. Geological Survey. Last accessed 02/20/2013. Available at: http://sofia.usgs.gov/publications/circular/1182/. .Accessed on 16 August 2013 Van Genuchten, M.T. 1980. A closed form equat ion for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44: 892 898. Vanclooster, M., P. Viaene, J. Diels, and K. Christiaens 1995. WAVE: A mathematical model for simulating water and agrochemicals inthe soil and vadose environment. Reference and user's manual (release2.0). Institute for Land and Water Management, Katholieke UniversiteitLeuven, Leuven, Belgium. Vesselinov, V.V., G. Pau, and S. Finsterle. 2012. AGNI: Coupling Model Analysis Tools and High Performance Subsurface Flow and Transport Simulators for Risk and Performance Assessments. In Proc. XIX International Conference on Water Resources, Computational Methods in Water Resources (CMWR 2012), Urbana Champaign, Illinois June 17 22, 2012.

PAGE 234

234 Viaene, P., H. Vereecken, J. Diels and J. Feyen 1994. A statistical analysis of six hysteresis models for moisture retention characteristic. Soil Sci. 157(6): 345 353. Wang, F. H. and M. P. Anderson, 19 82 Introduction to Groundwater Modeling: Finite Difference and Finite Element Methods. Academic Press Waswa, G. W., A. D. Clulow, C. Freese, P. A.L. Le Roux and S. A. Lorentz. 2013. Transient Pressure Waves in the Vadose Zone and the Rapid Water Table Response Vadose Zone J doi:10.2136/vzj2012.0054. Whiting L. M., L. Li, and S. L. Ustin. 2004. Predicting water content using Gaussian model on soil spectra. Remote Sensing of Environment 89 (4): 535 552. Wilsnack, M.M., Welter, D.E., Nair, S.K., Montoya,A.M., Zamorano, L.M., Obeysekera, J., Restrepo J.I., 2000. North Miami Dade County Groundwater Flow Model, Hydrologic Systems Modeling Division, South Florida Water Management District, West Palm Beach, FL. Available at: http://www.sfwmd.gov/portal/page /portal/xrepository/sfwmd_repository_pdf/nmd_tech.pdf Accessed on 23 July 2010. Winston, R.B. 2000. Graphical User Interface for MODFLOW, Version 4: U.S. Geological Survey Open File Report 00 315, 27 p. Winter, T. C., J.W. Harvey, O.L. Franke, and W.M. Alley. 1998. Ground water and surface water: a single resource. U.S Geological Survey Circular 1139, Denver. Zechner, E., and W. Frielingsdorf 2004. Evaluating the use of canal seepage and solute concentration obser vations for aquifer parameter estimation. Journal of Hydrology 289(1 4): 62 77. Zlotnik, V.A., and H. Huang.1999. Effect of partial penetration and streambed sediments on aquifer response to stream stage fluctuations, Ground Water 37(4): 599 605. Zou, S ., and Y. Yu, 1999. A dynamic factor model for multivariate water quality time series with trends. J. of Hydrol 178 (1 4) : 381 400. Zuur, A.F., R.J. Fryer, I.T. Jolliffe, R. Dekker, and J.J. Beukema 2003. Estimating common trends in multivariate time seri es using dynamic factor analysis. Environmetrics 14 (7) : 665 685. Zuur, A.F., E.N. Ieno, and G.M. Smith, 2007. Analysing ecological data Springer Verlag, Berlin. Zuur, A.F., and G.J. Pierce, 2004. Common trends in Northeast Atlantic squid time series. J Sea Res 52 : 57 72.

PAGE 235

235 BIOGRAPHICAL SKETCH Kisekka Isaya began his education at Rev. John Foundation primary school. Towards the end of his primary level education he moved to St. James primary school located in Kampala city, Uganda. In 1992, he started hi s secondary school education at for mathematics and physics. Great encouragement from his mother, friends and teachers, together with hard work helped Isaya to excel at the U ganda Certificate of Education (UCE) national exam. This was a turning point in his life, because of his good UCE score he was admitted in one of the oldest and top ranking secondary school in Uganda called Namilyango College. At Namilyango College, he maj ored in physics, chemistry and mathematics. In 1998, he took the Uganda Advanced Certificate of Education (UACE) national examination, through great support from his family and teachers he excelled again in the UACE exam. This enabled him to win a governme nt of University Kampala, Uganda. After his undergraduate degree in 2002, he worked as an irrigation engineer with Balton (U) Ltd, a place where he learned a great deal a bout the engineering profession. Three years later, he joined the National Agricultural Research Organization (NARO) of Uganda. He has also received several professional trainings in irrigation and drainage engineering in several countries including Israel Egypt and Kenya. In 2007, he decided to join graduate school and he was lucky to be accepted at the University of Florida department were his majoring in Agricultural and Biological engineering under the supervision of Dr. Migliaccio K. White. In 2009, h e completed a Master of Engineering degree majoring Agricultural and Biological Engineering at University of Florida under the same major supervisor.