HYDROGEOLOGIC AND BIOGEOCHEMICAL HETEROGENEITY ACROSS SPATIOTEMPORAL SCALES FLUXES, FLOWPATHS AND NITROGEN TRANSFORMATIONS IN AN EOGENETIC KARST AQUIFER By WESLEY R. HENSON A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2016
Â© 2016 Wesley R. Henson
To Christian, Betty, and Elizabeth your su pp ort and love made this possible
4 ACKNOWLEDGMENTS I thank my mother for encouraging me to always put my best foot forward and for instilling in me the desire to get an advanced education. I thank my grandmother for teaching me to read and strive for g reatness. I thank my husband Christian Fey for keeping me encouraged and supporting me through this process. I thank Wendy Graham for her helpful guidance every step of the way developing this research and for her critical evaluations of my scientific meth ods and arguments. She always was willing to give me a piece of her valuable time to keep me on track and moving forward. I thank Matt Cohen for giving me perspective on how to effectively write and analyze data to make clear arguments. I thank Jon Martin for his support of all of my field work and use of his Environmental Geochemistry Lab. I thank Mike Annable for help with tracer study design and materials. I thank Rafael Munoz Carpena for giving me an understanding of the importance of quantifying uncert ainty and developing novel hypotheses. I am thankful to Rob DeRooij, who assisted me in numerical methods, model development and was a willing mentor throughout the process. I would like to thank my Water Institute Graduate Fellow Faculty and Student Cohor t for giving me a new interdisciplinary lens through which the relevance of my work was made clear. I would like to thank the Suwannee River and St. Johns Water Management Districts for providing access to sites and data and supporting this work. Finally, I want to thank the staff of the Water Institute for providing research support, especially Mary Garvin who was always willing to address administrative tasks so I could focus on my research.
5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LI ST OF TABLES ................................ ................................ ................................ ............ 8 LIST OF FIGURES ................................ ................................ ................................ ........ 10 ABSTRACT ................................ ................................ ................................ ................... 12 CHAPTER 1 INTRODU CTION ................................ ................................ ................................ .... 14 Upper Floridan Aquifer System ................................ ................................ ............... 14 Florida Springs ................................ ................................ ................................ ........ 16 Nit rogen Research Needs ................................ ................................ ....................... 17 Dissertation Context and Objectives ................................ ................................ ....... 17 2 DRIVERS OF NITRATE TRANSFORMATION HOTSPOTS ACROSS AN AQUIFER C ONFINEMENT GRADIENT IN COMPLEX KARST TERRAIN. ............ 20 Background ................................ ................................ ................................ ............. 20 Methods ................................ ................................ ................................ .................. 23 Study Site ................................ ................................ ................................ ......... 23 Data Collection ................................ ................................ ................................ . 25 Well selection ................................ ................................ ............................. 25 Geo chemical sampling ................................ ................................ ............... 26 Nitrate isotopes ................................ ................................ .......................... 27 Dissolved gases ................................ ................................ ......................... 28 Data Analysis ................................ ................................ ................................ ... 28 Nitrate isotopes ................................ ................................ .......................... 28 Dissolved gas ................................ ................................ ............................. 30 Initial nitrat e concentrations ................................ ................................ ....... 31 Isotopic enrichment ratios ................................ ................................ .......... 32 Geospatial analysis ................................ ................................ .................... 33 Spatially distributed and spatiotemporally averaged nitrate transformation data ................................ ................................ ................. 35 Results ................................ ................................ ................................ .................... 35 Nitrate Isotopes ................................ ................................ ................................ 35 Dissolved Gases ................................ ................................ .............................. 36 Selected Geochemistry of Sites with Excess Nitrogen Gas .............................. 36 Geospatial Analysis ................................ ................................ .......................... 37 Spatially Averaged and Flow Averaged Denitrification Evidence ..................... 38 Discussion ................................ ................................ ................................ .............. 39
6 Utility of Nitrate Isotopes and Dissolved Gas Measurements in Karst Aquifers: ................................ ................................ ................................ ........ 40 Conditions That Facilitate Denitrification ................................ .......................... 41 Redox state and reactant availability ................................ ......................... 41 Aquifer confinement ................................ ................................ ................... 43 Land Surface aquif er connections ................................ ............................. 45 Spatially aggregated and flow averaged measurements ........................... 46 Study Summary ................................ ................................ ................................ ...... 47 3 NITRATE TRANSFORMATION MECHANISMS AND RATES IN AN UNCONFINED EOGENETIC KARST AQUIFER ACROSS A REDOX GRADIENT ................................ ................................ ................................ ............. 59 Background ................................ ................................ ................................ ............. 59 Study Objectives ................................ ................................ ................................ ..... 61 Methods ................................ ................................ ................................ .................. 62 Study Site ................................ ................................ ................................ ......... 62 Push Pull Tracer Tests ................................ ................................ ..................... 63 Microbial Analyses ................................ ................................ ........................... 67 Results and Discussion ................................ ................................ ........................... 69 Push Pull Tracer Tests ................................ ................................ ..................... 69 Microbial Analyses ................................ ................................ ........................... 71 Study Summary ................................ ................................ ................................ ...... 75 4 EXAMINING DRIVERS OF KARST AQUIFER HYDROLOGIC RESPONSE USING STOCHASTICALLY GENERATED PREFERENTIAL FLOW PATH TEMPLATES AND GLOBAL SENSITIVITY ANALYSIS ................................ ......... 84 Background ................................ ................................ ................................ ............. 84 Karst Aquifers ................................ ................................ ................................ ... 84 Karst Aquifer Development ................................ ................................ ............... 85 Karst Aquifer Model ing ................................ ................................ ..................... 86 Karst Conduit Evolution Algorithms ................................ ................................ .. 87 Sensitivity Analyses ................................ ................................ .......................... 88 Study Objectives ................................ ................................ ................................ ..... 90 Methods ................................ ................................ ................................ .................. 92 Overview of Flow Transport and Conduit Dissolution Model ............................ 92 Hydrologic Model ................................ ................................ ....................... 93 Boundary conditions ................................ ................................ .................. 94 Preferential flow path template generation ................................ ................. 95 Dissolution model ................................ ................................ ....................... 98 Morris Method Global Sensitivity Analysis ................................ ........................ 98 Disso lution Response Metrics ................................ ................................ ........ 100 Hydrograph Response Metrics ................................ ................................ ....... 101 Transport Response Metrics ................................ ................................ .......... 102 Behavioral Monte Carlo Filtering ................................ ................................ .... 102 Results ................................ ................................ ................................ .................. 103
7 Dissolution Behaviors and Sensitivity ................................ ............................. 103 Steady State and Hydraulic Pulse Response Sensitivity ................................ 106 Transport Pulse Response Sensitivity ................................ ............................ 106 Hydraulic Pulse Response Metrics ................................ ................................ . 107 Transport Pulse Response Metrics ................................ ................................ 108 Behavioral and Non Behavioral Networks ................................ ...................... 109 Discussion ................................ ................................ ................................ ............ 111 Study Summary ................................ ................................ ................................ .... 115 5 EVALUATING THE IMPORTANCE OF CONDUIT AND SINK LOCATION ON FLOW AND TRANSPORT IN AN IDEALIZED SILVER SPRINGS MODEL ......... 133 Background ................................ ................................ ................................ ........... 133 Study Obj ectives ................................ ................................ ................................ ... 133 Methods ................................ ................................ ................................ ................ 135 Results ................................ ................................ ................................ .................. 137 Hydrologic and Transport Res ponse ................................ .............................. 138 Backward Transport Pulse Monte Carlo ................................ ......................... 140 Discussion ................................ ................................ ................................ ............ 141 Study Summary ................................ ................................ ................................ .... 145 6 SUMMARY AND CONCLUSIONS ................................ ................................ ........ 157 Spatially Distributed Nitrate Transformations ................................ ........................ 157 In Situ Denitrification Measurements in Geochemical End Members ................... 159 Effects of Karst Aquifer Complexity on Fluxes and Flow Paths ............................ 161 Effects of R andomly Distributed Preferential Flow P aths. ................................ ..... 163 Insights for Further Research ................................ ................................ ............... 164 APPENDIX A : SELECTED GOECHEMICAL DATA FOR THE ICHETUCKNEE SPR INGSHED .. 165 B : KARST DISSOLUTION MODEL OVERVIEW ................................ ....................... 170 Di ssolution Model ................................ ................................ ................................ . 170 Theory ................................ ................................ ................................ ............ 170 Flow and Reactive Solute Transport ................................ .............................. 170 Calcite Dissolution ................................ ................................ .......................... 172 Model Design ................................ ................................ ................................ . 175 The quasi steady state approximation ................................ ..................... 175 Numerical Solution of Flow ................................ ................................ ............. 177 Numerical Solution of Reactive Transport ................................ ...................... 178 LIST OF REFERENCES ................................ ................................ ............................. 180 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 201
8 LIST OF TABLES Table page 2 1 Ichetucknee Spring Basin well information. ................................ ........................ 53 2 2 Estimates of excess N 2 (N 2den ), initial NO 3 N concentrations [ , residual fraction of original NO3 (f) observed at the sampling sit e. ................................ .. 53 2 3 Phase I and Phase II geochemical parameters for the four sites with the highest excess N 2 gas ([N 2 ] den ) values. ................................ ............................... 54 2 4 Confined and unconfined aquifer Mann Whitney U test results for selected geochemical and geospatial data. ................................ ................................ ...... 54 2 5 Excess nitrogen gas Mann Whitney U test results for selected geochemical and g eospatial data. ................................ ................................ ........................... 55 2 6 Comparison of coefficients of variation (CV) in measured geochemistry. ........... 57 3 1 Selected geochemical data f or characterization of push pull tracer test sites. .... 80 3 2 Primers for Nitrogen Functional Genes showing functional genes hazA , nrfA , nirK , and nirS . ................................ ................................ ................................ ..... 80 3 3 Push pull tracer test summary . ................................ ................................ ........... 81 3 4 Estimated pore water biomass nitrogen from assimilatory reduction of nitrate to ammonium in recovered injectate for all push pull tracer tests. ...................... 83 4 1 Horizontal and vertical preferential flowpath (HPF and VPF) and matrix properties for Morris Method Global Sensitivity Analysis. ................................ . 119 4 2 Comparison of flows and radii statistics at steady flux dissolution time and one and a half times the steady flux dissolution time. ................................ ...... 120 4 3 Morris Method Glo bal Sensitivity Analysis results showing behavioral and non behavioral group statistics for dissolution. ................................ ................. 121 4 4 Morris Method Global Sensitivity Analysis results showing behavioral and non be havioral group hydrograph statistics. ................................ ..................... 125 4 5 Behavioral Morris results showing behavioral and non behavioral group statistics. ................................ ................................ ................................ ........... 126 4 6 Behavioral Morris results showing transport metric summary for behavioral and non behavioral cases. ................................ ................................ ............... 127 5 1 Hydrologic and transport pulse metrics for behavioral cases from Chapter 4 w ith values similar to published values for Silver Springs. ................................ 146
9 5 2 Behavioral network properties from cases similar to Silver Springs from chapter 4 and selected parameters for chapter 5 Monte Carlo experiments. ... 146 5 3 Output statistic convergence results for Monte Carlo simulations. ................... 147 5 4 Head transect output stati stic convergence results for Monte Carlo simulations ................................ ................................ ................................ ........ 147 5 5 Complete ensemble mean and standard deviation for hydrologic and transport metrics. ................................ ................................ .............................. 147 5 6 Behavioral ensemble mean and standard deviation for Hydrologic and transport metrics. ................................ ................................ .............................. 149 A 1 Measured geochemical parameters for Phase I and Phase II sampling in the Ichetucknee Springs Basin. ................................ ................................ .............. 165 A 2 Nutrients and dissolved carbon measurements for Phase I and Phase II sampling in the Ichetucknee Springs Basin ................................ ...................... 166 A 3 Dissolved metal concentrations for Phase I and Phase II sampling in the Ichetucknee Springs Basin. ................................ ................................ .............. 167 A 4 Geospatial statistics for Ichetucknee Springs Basin spat ial analysis. ............... 168 A 5 Site names and nitrate isotope data for Phase I and Phase II sampling . .......... 169
10 LIST OF FIGURES Figure page 1 1 Hydrogeologic units in the Sil ver Springs basin . ................................ ................. 19 2 1 Ichetucknee Springshed Sampling Site map. ................................ ..................... 50 2 2 Spatial data for Ichetucknee Springs and Cody Escarpment. ............................. 51 2 3 Measured isotopes of nitrate ................................ ................................ .............. 52 2 4 Li near regression for all sites ................................ ................................ .............. 56 2 5 15 N, dissolved oxygen and enrichment ratios. ................................ ................................ ........................ 57 3 1 Push pull tracer test site map.. ................................ ................................ ........... 79 3 2 Push pull tracer test cumulative mass recovery results. ................................ ..... 82 3 3 Microbial genetic data showing baseline conditions a nd response to tracer additions . ................................ ................................ ................................ ............ 83 4 1 Conduit Evolution Algorithm. ................................ ................................ ............ 117 4 2 Silver springshed site location. ................................ ................................ ......... 118 4 3 Silver S pring model domain showing model grid and model surface topography elevations. ................................ ................................ ..................... 119 4 4 Example lnQ vs t plot showing hydrograph ................................ ...................... 120 4 5 Simulated springflow during dissolution. ................................ ........................... 120 4 6 Conduit flow during evolution. ................................ ................................ ........... 121 4 7 Morris Method Global Sensitivity Analysis dissolution statistics ....................... 122 4 8 Morris Method Global Sensitivity Analysis hydrologic pulse response metrics . 123 4 9 Mor ris Method Global Sensitivity Analysis t ransport pulse response metrics ... 124 4 10 Ensemble hydrograph ................................ ................................ ...................... 124 4 11 Behavioral ensembl e hydrograph ................................ ................................ ..... 125 4 12 Ensemble chemograph ................................ ................................ ..................... 126
11 4 13 Behavioral ensemble chemograph ................................ ................................ ... 127 4 14 Morris Method Global Sensitivity Analysis results. ................................ ........... 128 4 15 Behavioral Monte Carlo Filtering results. ................................ .......................... 129 4 16 Two non behavioral conduit networks ................................ .............................. 130 4 17 Three behavioral conduit networks with the same transport peak time ............ 130 4 18 Behavioral example with spatially heterogeneous potentiometric surface ........ 131 4 19 Behavioral example with relatively smooth potentiometric surface ................... 132 5 1 Potentiometric surface maps for Silver Springs ................................ ................ 146 5 2 Complete ensemble hydrograph. ................................ ................................ ...... 148 5 3 Comple te ensemble head transects ................................ ................................ . 148 5 4 Complete ensemble chemograph . ................................ ................................ .... 149 5 5 Behavioral ensemble hydrograph ................................ ................................ ..... 150 5 6 Behavioral ensemble head transects ................................ ................................ 150 5 7 Behavioral ensemble chemograph ................................ ................................ ... 151 5 8 Equivalent porous media (EPM) mode l hydraulic heads and moments. . ......... 152 5 9 Behavioral Monte Carlo ensemble model hydraulic heads and moments ........ 153 5 10 Behavioral case 79. ................................ ................................ .......................... 154 5 11 Behavioral case 83 ................................ ................................ ........................... 154 5 12 Backwards transport pulse for case 79 ................................ ............................. 155 5 13 Backwards transport pulse for case 83 ................................ ............................. 155 5 14 Forward transport pulse for case 79 and 83 showing tracer mass flux. ............ 156
12 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy HYDROGEOLOGIC AND BIOGEOCHEMICAL HETEROGEN EITY ACROSS SPATIOTEMPORAL SCALES FLUXES, FLOWPATHS AND NITROGEN TRANSFORMATIONS IN AN EOGENETIC KARST AQUIFER By Wesley R. Henson August 2016 Chair: Wendy D. Graham Major: Agricultural and Biological Engineering The Upper Floridan Aquifer is hydrogeolo gically complex; limestone dissolution has led t o vertical and horizontal preferential flow paths. Aquifer nitrate contamination is of concern and is subject to regulatory limits. Research objectives were to examine spatially distributed nitrate transforma tion evidence in groundwater samples to inform where and under what condi tions transformations occur; quantify in situ aquifer nitrate transformation rates and mechanisms using tracer injection experiments combined with microbial genetic sampling of key ge nes i n nitrogen transformation pathways; examine how aquifer geologic heterogeneity influences conduit and spring development in karst systems ; and evaluate sensitivity of water and solute fluxes and flow paths to porous matrix and conduit properties in a springshed representative of Silver Springs. Excess N 2 in dissolved gas groundwater s amples confirmed denitrification in hotspots where the aquifer is unconfined. Observed denitrification under unfavorable local pore water conditions, suggests denitrificat ion may occur in localized microsites or biofilms. Suggestive relationships occurred between denitrification evidence and wetlands, however no relatio nship was found for sinks, swallets, and closed topographic
13 depressions. Local scale heterogeneity in nitr ogen cycling led to inconclusive isotopic results despite reasonable agreement between spatially averaged and spring vent measurements. Nitrate, carbon, and chloride tracer injection experiments showed nitrate losses under porewater conditions expected to be unfavorable for denitrification. Microbial genetic data confirmed multiple nitrate transformation pathways, including assimilatory reduction of nitrate to ammonium, dissimilatory reduction of nitrate to ammonium (DNRA) and denitrification . This contends with prior research, which has asserted DNRA is a negligible nitrate transformation pathway in freshwater aquifers. However, negligible spring vent ammonium concentrations support that ammonium produced by DNRA may be subsequently reduced by biologic ammo nium retention or conversion to nitrate in oxic zones. Karst geologic heterogeneity was examined by combining stochastically generated preferential flow path temp lates and simulation of karst aquifer genesis via dissolution. F irst magnitude sprin g developm ent increased with lower hydraulic conductivity and higher preferential flow path connectivity. Monte Carlo simulations of conduit generation, groundwater flow and conservative solute transport showed that exact knowledge of conduit location was less impo r tant than representation of prefe rential flow pathways for predicting solute travel paths and travel times to the spring and vulnerable springshed areas .
14 CHAPTER 1 INTRODUCTION The prosperity of our modern society depends on development of land for urba n and agricultural uses, from which loss of nitrogen to hydrologic systems is often an undesirable consequence. We must find ways to manage water resources to create prosperous communities while preserving the environment for future generations. Water mana gers are tasked with trying to meet ever reducing nitrogen loading targets in order to protect the aquifers that supply their communities and preserve downgradient aquifer connected ecosystems such as springs, rivers, lakes, and estuaries. In the United St ates , regulatory mechanisms such as Total Maximum Daily Loads (TMDL) and Numeric Nutrient Criteria (NNC) have been established to protect water resources from impairment due to excess nitrogen loading. This research aims to increase understanding about flu xes and flowpaths of nitrogen, and mechanisms and magnitudes of nitr ogen attenuation, in eogenetic karst aquifers, where significant aquifer flow occurs in both the porous matrix and in karst conduit networks. Improved understanding of these processes coul d help to more efficiently achieve regulatory nitrogen targets by determining dominant nitrogen attenuating processes over redox and aquifer confinement gradients, quantifying nitrogen reaction rates, and examining vulnerability imposed by uncertain ground water flow paths inherent in karst terrain. Improved understanding of these processes would contribute to improved vulnerable area delineation and nitrate load management. Upper Floridan Aquifer System The Upper Floridan Aquifer (UFA), an eogenetic karst aquifer in the southeast United States that covers approximately 260,000 km 2 and provides municipal, industrial,
15 and agricultural water supply for a population of approximately 1.8 million people in Florida alone [ Marella , 2008] , provides the geographic context for this study. The UFA in north Florida is mainly comprised of Ocala Limestone, ranges in thickness from 0 to 180 feet, and is underlain by a lower permeability limestone cal led the Avon Park Formation that is 800 1100 feet thick (Figure 1 1). Where the Ocala L imestone is not exposed at the surface, it is covered by the Hawthorn Format ion and a surficial aquifer of Plio P leistocene sands [ Scott et al. , 2004] . Spatially variable erosion of the Hawthorn Formation has led to variations in UFA confinement throughout North Florida. The erosional boundary of the Hawthorn Formation is known as the Cody Escarpment, and defines a critical boundary for defining UFA vulnerability to contamination. Downgradient from the Cody Escarpment, where limestone is exposed, vulnerability is increased by enhanced surface to a quifer connectivity. Nitrate contamination can be mitigated by surface water retention and biological processes in wetlands, stream channels, and lakes . Where the aquifer is unconfined, these surface water features are reduc ed significantly and karst windo ws and sinks provide direct connections between the surface and aquifer . In the UFA, aquifer high permeability, high porous matrix and fracture porosity, and conduit networks lead to highly heterogeneous hydrologic flow paths [ Hornsby and Ceryak , 1998; Martin and Dean , 2001; Florea and Vacher , 2007a] . In karst aquifers, g roundwater flow ranges from laminar (porous media flow) to turbulent (pipe flow), depending on the degree of void and conduit development [ Atkinson , 1977] . Dissolution processes create sinks, caves, and networks of conduits within the soluble porous matrix, with properties and geometries that are poorly characterized [ Ford and Williams ,
16 2013] . Karst landscape features such as sinks, sinking streams (swal lets), and breaches in aquifer confinement lead to significant surface to aquifer connectivity, and can route significant surface flow into the aquifer [ Ewers , 2006] complicating contaminant source delineation and further contributing to hydrologic flow path complexity. Furthermore, contam inant transport is highly influenced by porous matrix and conduit interactions. Florida Springs Northern Florida contains the highest density of large springs in the world; with 33 first magnitude (discharge greater than 2.8 m 3 /s) and 191 second magnitude (discharge between 0.28 2.8 m 3 /s) springs [ Scott et al. , 2004] . Springs are a unique class of ecosystem, with highly adapted pla nt and animal communities [ Munch et al. , 2006] . In addition, springs are often used for recreation and contribute significantly to the North Florida econo my [ Bonn and Group , 2004] . Thus their preservation is a high priority. Springs are often connected to extensive conduit systems and have large groun dwater contributing areas referred to as springsheds. In this way, springs are sentinels for aquifer response to human activities over large spatial and temporal scales. Accordingly, spring vent flow and water have been widely used to characterize karst aq uifers [ Phelps and Survey , 2004] , groundwater source mixing [ Toth and Katz , 2006] , and nitrogen attenuation [ Katz et al. , 2001; Cohen et al. , 2012; Heffernan et al. , 2012] . In North Florida, spring fed rivers are connected to major river systems such as the Suwa nee River and St. Johns River. Thus, human activities resulting in nitrogen loading to springs have far reaching impact on downgradient water bodies. Springs persist in an approximate steady state and can be considered chemostatic [ Odum ,
17 1956] compared to most other aquatic ecosystems. However anthropogenic factors, including extraction of groundwater for consumptive uses and nitrogen loading from agricultural and urban land uses, is contributing to significant changes in these systems. Significant increases in nitrate concentrations (i.e. 10 to 350 fold above background concentrations) have been observed in some springs [ Brown et al. , 2008] , and potentially contribute to shift s in ecosystem ecology in favor of algae over rooted plants, and increased algae in streams [ Rabalais , 2002] . These changes impair the recreational and aesthetic use of springs and threaten springs driven tourism. Nitrogen Res earch Needs It is well understood that karst aquifers are particularly vulnerable to loading of excess nitrate nitrogen from land management activities [ White et al. , 1995; Coxon , 1999] . Current management efforts assume that aquifer denitrification in the UFA is negligible over large spatial scales and that nitrate nitrogen behaves like a conservative solute [ Phel ps and Survey , 2004; Katz et al. , 2005] . However, recent research has called this assumption in question [ Heffernan et al. , 2012] , indicating the need for more information about the fate and transport of nitrogen in the UFA. Understanding nitr ogen contamination source areas and fluxes and flowpaths with the aquifer, is critical to springshed management [ Brown et al. , 2008] . In addition, understanding the assimilative capacity of the aquifer for nitrate is a research priority. Dissertation Context and Objectives The goal of this study is to examine nitrate transport, retention, and attenuation in karst a quifers using a multi disciplinary approach that combines measurements of ambient groundwater biogeochemistry, nitrate isotopes, dissolved gases; in situ measurements of nitrate transformation rates and mechanisms; physical modeling of
18 karst evolution; hydr ologic modeling flow and transport in karst systems; and global sensitivity analysis methods. Specifically , the objectives of this research were to (1) examine spatially distributed nitrate transformation evidence in Upper Floridan Aquifer (UFA) springshed groundwater samples to inform where and under what conditions transformations occur; (2) quantify in situ rates and mechanisms of nitrate transformation in the UFA using nutrient injection tracer experiments combined with microbial genetic sampling; (3) e xamine how aquifer geologic heterogeneity influences conduit and spring development in karst systems; and (4) evaluate the sensitivity of water and solute fluxes and flowpaths to porous matrix and conduit properties in a representative springshed. Chapte r 2 examines where in the springshed nitrogen transformations are occurring, what landscape features are associated with areas of higher nitrate attenuation in springsheds, and compares spatially distributed geochemistry and nitrate transformation evidence to samples previously collected at spring vents emanating from the UFA. Chapter 3 investigates in situ rates and mechanisms of nitrate transformations across a karst aquifer gradient from oxidizing to reducing conditions using combined tracer tests and mi crobial genetic data. Chapter 4 examines the evolution of conduit development in karst aquifers, investigating what geologic features lead to development of first magnitude springs in eogenetic karst systems, and explores the influence of these geologic fe atures on network development, springflow magnitude, and hydraulic and transport pulse response metrics. Chapter 5 addresses the uncertainty of hydraulic and transport pulse response contributed by uncertain preferential flow path location and orientation with fixed porous matrix and preferential flow path statistical parameters.
19 Figure 1 1. Hydrogeologic units in the Silver Springs basin [Phelps and Survey, 2004; Scott et al., 2004] .
20 CHAPTER 2 DRIVERS OF NITRATE TRANSFORMATION HOTSPOTS ACROSS AN AQUIFER CONFINEMENT GRADIENT IN COMPLEX KARST TERRAIN. Backgro und Nitrogen in nitrate form is considered an abundant, mobile contaminant in shallow groundwater [ BÃ¶hlke , 2002] . Nitrate loading is spatiotemporally variable [ Rivett et al. , 2008] . Karst aquifers are particularly vulnerable to nitrate contamination from agricultural activities [ White et al. , 1995; Coxon , 1999] , especially where the y are unconfined and shallow. In the Upper Floridan Aquifer (UFA), background nitrate concentrations have increased significantly [ Katz , 1992, 1999; Hornsby and Ceryak , 1998] increasing discharge from springs to surface waters [ Pittman et al. , 1997] and potentially contributing to significant development of algae in many spring fed rivers and downgradient aquatic ecosystems [ Jones et al. , 1996; Hornsby and Mattson , 1998] . The UFA is the uppermost aquifer in the Floridan Aquifer System [ Scott et al. , 2004] . The UFA is an eogenetic, shallow limestone karst system with high matrix permeability that is the source for many North Florida springs [ Berndt et al. , 1998] . The UFA is geologically complex; including well developed karst conduits, highly permeable l imestone and spatially varying aquifer confinement. This complexity leads to heterogeneous aquifer properties and surface to aquifer connectivity [ Hornsby and Ceryak , 1998; Katz et al. , 2009] . Understanding how landscape and aquifer chemical heterogeneity influence patterns in nutrient transport and transformation is cruc ial for developing biogeochemical models that honor terrestrial and aquatic characteristics [ Grimm et al. , 2003] and interpreting nitrate transformation evidence. Karst spring measurements have been widely used to examine processes occurring throughout the heterogeneous landscape. However, more information about
21 relationships between spring measurements and distributed denitrification measures and the role of karst landscape features on observed denitrification is needed. Prior research supports that significant nitrate attenuation is occurring in the UFA [ Katz et al. , 2004; Heffernan et al. , 2012] ; however defining where signifi cant nitrate attenuation is occurring remains elusive [ McClain et al. , 2003; Brown et al. , 2008] . In the UFA; the focus of research has been to examine temporal variation in nitrate attenuation processes (e.g. denitrification ) using spring vent geochemistry and dissolved gas analyses and isotopes [ Katz et al. , 2001, 2004; Heffernan et al. , 2012 ], nitrate loading estimates and budget methods [ Katz et al. , 2009] , and groundwater age distributions [ Katz , 1999; Katz et al. , 2001] . While these methods describe nitrate attenuation over large spatiotemporal scales, local scale denitrification measures in karst aquif ers are springshed behaviors. A goal of this work is to examine how spatially distributed measurements relate to aggregate responses, increasing understanding of how point measurem ents link to highly convolved flows. Distributed springshed measures may illuminate poorly described local phenomena that may explain documented variation in denitrification rates across UFA springs [ Heffernan et al. , 2012] . One potential mechanism that distributed measures may elucidate is the influence of confining units on denitrification through their role in focusing surficial flows and reactants to the con fined UFA, where conditions may be favorable for denitrification, via breaches in the confining layer. Shallow water tables and poor drainage, as observed surficial aquifers in confined regions, increase the likelihood of denitrification [ BÃ¶hlke , 2002] . Furthermore, confined aquifers have been suggested as a
22 favorable places for denitrification reactions, both in the low hydraulic conductivity confining layer and at the interface between confined and unconfined aquifers [ Rivett et al. , 2008] . Aquifer confinement generally restricts delivery of reactants to the aquifer, however flow focusing mechanisms such as sin ks and swallets that breech the confining unit may deliver more reactants to favorable conditions for denitrification. Thus areas in the aquifer where reactants (nitrate and electron donors) are delivered to the right conditions (anaerobic Â± solid phase el ectron donors) may correlate with areas of si gnificant denitrification . The goal of this research was to investigate how karst aquifer flowpaths control denitrification rates by exerting controls on local aquifer geochemistry and delivery of reactants. It was hypothesized that: (1) that groundwater redox state is a good predictor for denitrification; (2) degr ee of aquifer confinement influence s the location and magnitude of denitrification; and (3) spatial featu res that focus surficial flows and route them to the confined aquifer (e.g. sinks, swallets, closed topographic depressions and wetlands) influence delivery of reactants and denitrification hotspots. To test these hypotheses, denitrification was characterized in the UFA using synoptic spatially distr 15 18 O), excess dissolved nitrogen gas, and geospatial data describing land surface aquifer connections. In addition, spatially distributed denitrification measures were compared with pr eviously reported spring vent denitrification measures in an effort to spatially disaggregate the integrated springshed denitrification evidence that spring vents provide.
23 Methods Study Site The study site is the Ichetucknee Springshed in North Florida (F igure 2 1) which is underlain by the UFA. The study site is a typical North Florida UFA karst basin, containing a gradient of aquifer confinement and karst topographic features that influence groundwater surface water connections, and prior studies provide a wealth of spring vent measurements [ Katz , 2004; Cohen et al. , 2012; Heffernan et al. , 2012] . These factors provide a good setting for examining how the degree of UFA confinement and land surface aquifer connectivit y influence the location and magnitude of denitrification and how spatially distributed denitrification measures compare to spring vent measures. The study area lies between two physiographic provinces; the Northern Highlands where the UFA is confined and Gulf Coastal Lowlands where the UFA is semi confined to unconfined. The Northern Highlands are characterized by poorly to moderately drained sands and low permeability clay rich sediments which overlay the UFA [ Houston , 1965; Ceryak et al. , 1983] leadin g to numerous lakes, streams, and wetlands. The transition between the two provinces is defined by an erosional feature called the Cody Escarpment and is characterized by sinking streams, sinkholes and sinkhole lakes. Surface water is quickly drained to th e UFA through these karst solution features throughout in the North Highlands Province and the transitional area. In the Gulf Coast Lowlands province the UFA is unconfined, overlain by a thin layer of well to excessively drained sands [ Ceryak et al. , 1983] . This leads to much fewer lakes, wetland s , or surface water drainage features.
24 Land u se in the study area is a mix of forest (52%), agricultur e (27%), urban (12%) and su rface water and wetlands (9%). Local nitrogen loads consist of fertilizers (51%), animal waste (27%), septic tanks (12%), atmo spheric deposition (8%) and land application of reclaimed wastewater (2%) [ Katz et al. , 2009] . The local climate is subtropical with long warm summers and mild winters. Groundwater flow is predominantly from north to south. Locally, the UFA is characterized by high transmissivity values that can r ange up to several million m 2 d 1 [ Bush and Johnson , 1988] . These values reflect the thick limestone aquifer and high permeability of connected voids, conduits and caves due to dissolution of the limestone matrix. Groundwater is primarily drained by the Ichetucknee Springs Group . The Ichetucknee Springshed is approximately 960 km 2 , as estimated by Suwanee River Water Management District (SRWMD) using water year 2003 poten tiometric surfaces to determine groundwater contributing area [ Sepulveda et al. , 2006] , a nd further validated by Champion and Upchurch . A prior study by Katz et al. , suggested there is minimal inter annual change in groundwater storage or deep inter basin transfer in the springshed . These prior studies establish the Ichetucknee S pringshed as a relatively well constrained groundwater contributing area for examining nitrate transformations within the UFA. The Ichetucknee Springs Group is a collection of many springs, including 1 first magnitude (>2.8 m 3 /s) and 4 second magnitude spr ings (flows between 0.28 and 2.8 m 3 /s) , and hundreds of unnamed seeps [ Scott et al. , 2004] . Median residence times for groundwat er discharging from the springs, estimated with binary and piston flow groundwater age modeling, is approximately 10 25 years [ Katz et al. , 2001] . During the
25 period from 2002 to 2011, mean annual discharge from the two largest springs in the group , provided about hal f of the total annual mean discharge of 8.35 m 3 /s for the Ichetucknee River [ USGS , 2014] . The present day Ichetucknee River flows southwest from the springs approximately 7 km to the Santa Fe River. In the past, the Ichetucknee River connected Lake City, F L to the present day river channel [ FDEP , 2012] . All former surface tributaries are now terminal sinks (e.g. Cannon Sink, Clay Hole Sink, and Rose Sink) directly connected to the UFA through a karst conduit below the dry river channel [ Katz et al. , 2009] . Further detailed geomorphic and hydrologic descriptions of Ichetucknee Springs and Ichetucknee River can be found in [ Rosenau et al. , 1977; Scott et al. , 2004] . Data Collection Synoptic well sampling methods investigated factors that contribute to denit rification in the UFA; major element chemistry, aquifer denitrification proxies, and dissolved gases. The spatial patterns in measurements data were evaluated with geospatial data describing land surface aquif er connections common to karst terrane such as sinks, swallets, wetlands, aquifer confinement gradients and vulnerability, and closed topographic depressions. Well s election Wells were selected from the Suwannee River Water Management District (SRWMD) monitoring well network across the geologic gradien t from confined to unconfined aquifer condition. Well selection focused on capturing the spatial heterogeneity in karst features influencing hydrologic flow paths and residence times. The number of wells was limited based upon private land well access and instrumented wells. Borehole access , the ability to gain access to the top of the well, is req uired for
26 dissolved gas sampling . Eleven wells with borehole access and five wells without borehole access were used in this study. Two aquifer sampling events we re conducted; herein described as Phase I and Phase II. Phase I sampling occurred during the wet season in July August 2013. It included measurements of geochemistry and isotopes of nitrate ( 15 N and 18 O) in 16 wells (Figure 2 1 , Table 2 1 ). Phase II sampli ng occurred during the dry season (February March 2014), and included all Phase I measurements and added dissolved gas sampling f or selected wells ( which excluded sites 2, 4, 6, 7, 11) (Table 2 1) . Wells depths in the confined region were substantially dee per than wells in the unconfined region due to the thickness of the surface confining unit overlying the UFA. Geochemical s ampling All water samples were collected after removal of at least three well volumes and stabilization of field parameters. A broad suite of solutes, were used to characterize groundwater from each well and examined for covariance with denitrification evidence. Geochemistry measurements included dissolved oxygen (DO), specific conductance (SpC), oxidation reduction potential (ORP), p H, total kjeldahl nitrogen (TKN), orthophosphate (PO 4 ), nitrate + nitrite (NO x ), a mmonium (NH 4 ) and a suite of dissolved ions (i.e. sodium, calcium, potassium, magnesium, sulfate, fluoride, chloride, iron, strontium, magnesium), Dissolved Organic Carbon (D OC), Dissolved Inorg anic into prepared glass and plastic bottles, preserved, and placed on ice to prevent degradation. Dissolved oxygen, SpC, ORP, and pH were measured using a calibrated YSI 556 multi probe. Nutrients were analyzed at the University of Florida Analytical research
27 Lab; TKN was measured using EPA method 351.2, orthophosphate using EPA method 365.1, NO x using EPA method 353.2, and NH 4 using EPA method 350.1. Major element che mistry was measured with an automated dionex DX500 Ion Chromatograph with precision of better than 3% for all runs and elements based on internal check standards every fifth sample. Charge balance errors were less than 3%. Metal concentrations (iron, stron tium, manganese) were measured using a HR ICP MS Element 2 (Thermo Finnigan). The precision and accuracy of these samples were calculated by comparing measurements from the beginning and end of each run to the external standard SLRS4 (Canadian River water standard). Internal instrument drift was corrected using a Re Rh spike added to all standards and samples. DOC samples were analyzed on a Shimadzu TOC 5000A total organic carbon analyzer after sparging for 2 minutes with carbon free air to remove inorganic the injections is the reported value. Replicates had less than 5% coefficient of variance. Dissolved inorganic carbon (DIC) was measured on carbon dioxide (CO2) extracted by acidifying samples using an AutoMate Prep Device coupled with a UIC (coulometrics) 5011 carbon coulometer. Nitrate i sotopes Stable isotopes of nitrate (i.e., 14 N and 15 N) and O (i.e., 16 O and 18 O) within the nitrate ion and dissolved gases in the groundwater samples were mea sured. Nitrate rinsed plastic bottles, sealed and frozen. Samples were analyzed at the University of California, Riverside Facility for Isotope Ratio Mass Spectrometry (FIRMS) using the bacterial conversion of nitrate to nitrous oxide method [ Casciotti et al. , 2002; Coplen et al. , 2007] with maximum two 15 N 18 O .
28 Dissolved g ases Dissolved gas samples provide direct evidence for denitrification. Excess nitrogen gas (N 2 ) above concentrations that can be explained by excess air entrainment and gas solubility in water at the time of recharge indi cate subsurface generation of nitrogen gas through nitrate transformation processes such as denitrification. Excess N 2 was incorporated into geochemical and geospatial analyses because it is the most direct evidence of denitrification. Excess N 2 was estim ated using dissolved gas samples collected during Phase II using advanced diffusion sampling devices (ADS) [ Gardner and Solomon , 2009] . After removal of 3 well volumes and attainment of field parameter stability, each ADS was left to equilibrate with the aquifer at each site for 48 hours in the middle of the o pen interval before retrieval. Dissolved gas samples were extracted from each ADS at th e University of Utah Dissolved and Noble Gas Lab and examined using a Magnetic Sector Field Mass Spectrometer. Reported analytical precision is plus or minus 4% for nitrogen gas concentrations [N 2 ], 2% for Argon concentrations [Ar], 2% for Neon concentrati ons [Ne], and 3 4% for Xenon concentrations [Xe]. Data Analysis Nitrate isotopes Dual isotope measurements allow the inference of denitrification based upon isotopic changes in enrichment of heavier 15 N and 18 O isotopes that occurs during denitrification a nd can be used to distinguish between assimilation and denitrification enrichment mechanisms and allows for better resolution of nitrogen source and denitrification evidence [ Amberger and Schmidt , 1987; Kendall and Aravena , 2000; Kendall and McDonnell , 2012] . To evaluate isotopic evidence two diagnostic measures
29 were used; the slope of the r elationship between 15 18 O values and isotopically defined N source ranges from prior research [ Kendall and Aravena , 2000; Kendall and McDonnell , 2012] . Slope of the relationship between 15 18 O values in spatially distributed samples were compared to prior spr ing vent measurements [ Martin et al. , n.d.; Hefferna n et al. , 2012] . Along distinct hydrologic flow paths experiencing denitrification 15 18 O values increase exponentially with decreasing nitrate concentration, 15 18 O enrichment ratios along a flow path of approximately 2:1 [ Vogel et al. , 1981; BÃ¶ttcher et al. , 1990; Cohen et al. , 2012] . With a 1:1 isotopic enrichment ratio, assimilation and denitrification cannot be distinguished. Thus the slope of the 15 18 O values in spatially distributed samples is useful for evaluating denitrification. 15 18 O isotope values can also be compared to isotopic values in known nitrogen pools, e.g. synthetic fertilizer, soil nitrogen, and human and a nimal waste, to infer potential nitrogen sources [ Kendall and Aravena , 2000] ( Figure 2 15 N can have significant overlap between both inorganic and organic N source ranges [ Kendall and McDonnell , 2012] , how ever values that plot outside the range of nitrogen source isotopic values may have been subject to denitrification. Significant denitrification (>15%) has to occur before a nitrate isotope signature will plot outside of the source regions. Interpretation 15 N isotopic data is challenging due variation in nitrogen sources in recharge caused by mixing of soil N reservoirs in the soil zone and N transformations that fractionate the source signals
30 [ Ke ndall and Aravena , 2000] . Therefore, other evidence is needed to confirm denitrification such as excess dissolved nitrogen (N 2 ) gas. Dissolved g as Atmospheric gas solubility in groundwater is influenced by physical conditions during infiltration; most i mportantly the recharge temperature and excess air entrainment. Dissolved neon and argon gas concentrations in groundwater samples ([Ne] and [Ar]) were used to quantify variability in recharge temperatures ( ) and excess air mass ( ) dissolved in groundwater recharge, allowing estimates of potential denitrification N 2 enrichment in the absence of degassing and excess air fractionation (Plummer et al., 2001). and can be estimated sim ultaneously using [Ar] and [Ne], which have well constrained solubility temperature relationships and provides an estimate of the excess air entrained in the aquifer water. Measurements of dissolved noble gases concentrations such as Xenon [Xe] can improve estimates of and [ Aeschbach Hertig et al. , 1999a, 1999b] . Dissolved gas concentration at the time of recharge in a groundwater sample ( ) can be determined by simultaneous solution of 3 equations [ Hamme and Emerson , 2004] . The solubility constant , , as a function of recharge temperature ( : (2 1) the equilibrium solubility of other dissolved gases (e.g. , , and ) in water at temperature at the rechar ge elevation: (2 2)
31 where has units of milligrams of gas per liter of water , are temperature gas solubility constants; and t he expected total concen tration for each gas as a function of equilibrium concentration and excess air: (2 3) where 1 to mg L 1 ) and is the partial pressure for each gas co ncentration. Excess air mass is estimated by solving equations 1 3 for each measured gas, minimizing the sum of squared errors between modelled and observed data. The estimated magnitude of d enitrification is the difference between measured N 2 co ncentrations ( ) and predicted N 2 concentration at the time of recharge ( ): . (2 4) concentrations ( mg L 1 ) during recharge are a function of gas s olubility at the recharge temperature, excess air mass, and partial pressure of N 2 in the atmosphere computed as (2 5) where k is 0.028, is the partial pressure of N 2 in the atmosphere ( 0.78084), is solved using e quations 2 1 and 2 2, and and are determined using [Ne], [Ar], and [Xe]. Initial nitrate concentrations Excess N 2 ( was used to estimate initial nitrate concentrations ([NO 3 ] 0 ) and residual fraction of nitrate remaining ( f ). Excess N 2 measurements were combined with isotopic data to compute 15 N enrichment ratios ( ) and estimate initial nitrogen
32 source ( 15 N NO 3 ] init ) based on enrichment ratios and measured 15 1 5 N NO 3 ] obs ) (Table 2 6). The value for was used to determine the original nitrate concentration before denitrification at the time of recharge ( ): (2 6) where is the measured nitrate concentration in each groundwater sample. The value includes nitrate derived from nitrification in the unsaturated zone or UFA. Prior research contends that denitrification is the only significant n itrate transformation; i.e. that dissimilatory nitrate reduction to ammonium (DNRA) and direct nitrate assimilation are negligible based on concentrations of ammonium and particulate and dissolved organic nitrogen (DON) below detection limit at spring ven ts [ Heffernan et al. , 2012] 15 N NO 3 are assumed to be zero; spatially distributed measures in this study confirm ammonium concentrations below detection. Isotopic enrichment ratios by combining geochemical, isotopic, and dissolved gas measurements, and assuming Rayleigh distillation kinetics [ Mariotti , 1983; BÃ¶hlke , 20 02; Green et al. , 2008] . This assumes that areal rates of N inputs are equivalent and of identical isotopic source signatures. The land use in the springshed is relatively homogenous with fertilizer application to crops and improved pasture as the dominan t source of N [ Katz et al. , 2009] and prior research asserts this assumption is reasonable [ Heffernan et al. , 2012] . However, local variations in N source may be detected in spatially distributed measure ments leading to a larger
33 range of uncertainty in the isotopic enrichment factor. This uncertainty in enrichment was considered for the estimated N sources from samples. Original nitrogen source ( ) 15 N ( ), estimated isotopic enrichment factor ( ) and the residual fraction ( ) using the following relation: (2 8) where the residual fraction ( ) of original nitrate was determined as: (2 9) This assumes a constant fractionation in space and time and ignores potential mixing effe cts (Green et al., 2010). Source N values were compared to known ranges for [ Heaton , 1986] . While N sources may vary in space and time, these values provide a benchmark for evaluating general source distribution. Geospatial analysis To test the hypot heses that aquifer confinement and spatial distribution of karst flow focusing mechanisms influence the spatial heterogeneity in denitrification, a suite of geographic datasets was assembled. One data set focused on characterizing aquifer confinement using the generalized aquifer confinement map for Florida (Figure 2 1A), [ USGS , 1998] , intermediate confining geologic unit thickness data, 2014 (Figure 2 1B) [ SJRWMD , 2014] , and the US Environmental Protection Agency DRASTIC index (Figure 2 2A) [ Aller et al. , 1987; FDEP , 1990] . The DRASTIC index is a tool for understanding aqu ifer susceptibility to contamination. The index includes depth to water, net recharge, soil media, topography, aquifer media, and hydraulic conductivity.
34 The second data set included karst features that collect flow and reactants such as sink and swallet l ocations [ FGS , 2 003] , wetlands [ USFWS , 2014] , and closed topographic depressions (CTDs) (Figure 2 2B) [ FGS , 2004] . Metrics evaluated for each sampled well include: distance to nearest sink, swallet, CTD, or wetland, number of CTDs or wetlands within 0.5 km, total area of CTDs or wetlands within 0.5 km, confining layer thickness and DRASTIC Index. Geospatial data were analyzed to investigate the correlation between features that influence land surface aquifer connectivity and denitrification proxies ( ). Nitrate isotope data are subject to several fractionation processes between the time of infiltration and aquifer recharge that can confound inference of denitrification. Therefore, only excess dissolved nitrogen gas ( ) da ta were examined with spatial and geochemical data. Linear regression was used to examine the correlation between geospatial information, aquifer geochemistry ( redox state , nitrate availability, electron donor availability), and dissolved gas denitrificati on measurements ( ) . For further evaluation of spatial trends across the collected data, chemistry, , and geospatial data were spatially aggregated into two groups for comparative analysis; sites with and without detection a nd sites above and below the Cody Escarpment (i.e. confined and unconfined). The non parametric Mann Whitney U test was used to test the null hypothesis that the distributions of both populations were equal and whether the mean parameters of each group wer e significantly different. Statistical power was computed using G*Power [ Faul et al. , 2007] .
35 Spatially distributed and spatiotemporally averaged nitrate transformation d ata To examine h ow heterogeneous nitrate transformations at the local scale integrate to the basin scale, , initial nitrate at time of recharge ([NO3] 0 ) and residual nitrate fraction ( f ) from distributed synoptic well samples were compared to the spring vent metrics for the springs reported by Heffernan et al. , [2012 ]. Nitrate isotope 15 18 O in nitrate isotopes (Figure 2 3) , 15 N and DO (Figure 2 15 N isotopic enrichment ratios (Figure 2 15 N source values. Results Nitrate Isotopes Nitrate isotopes collected in Phase I and II indicated a significant influence of fertilizer and soil nitrogen derived nitrate wi th minimal influence of animal and human waste (Figure 2 15 18 O values ranged from 0.55 to 10.0 (Table 2 2) ; consistent with groundwater nitrate from cultivated lands [ Bohlke and Denver , 1995; Kendall and Aravena , 2000] . Values of 15 N less than 6 are consistent with fertilizer sources. Significant correlation between nitrate isotopes 15 18 O), with a slope for their association of 1.41 (Figure 2 3), which was much higher than slopes reported elsewhere, suggest 18 O enrichment from nitrification was influencing the expected association. Up to one third of the oxygen isotopes in nitrate can come from atmospheric sources [ Snider et al. , 2010] . 18 O values are approximately 24 [ Luz and Barkan , 2011] , leading to o 18 15 N. 15 N source values in this study ranged from (Table 2 2) . These values overlap with previously identified ranges for fertilizer and are
36 comparable to measured spring 15 N values ( 0.6 to Heffernan et al. (2012), and for values previously reported by [ Katz et al. , 2004, 2009] , for s prings across in North Florida Dissolved Gases Excess nitrogen gas ( ) values ranged from 1.10 to 1.37 mg L 1 . The values were above zero for half of the sampled sites, ranging from 0.08 to 1.37 mg L 1 (Table 2 2 ). V alues of (>0.55 mg L 1 ) were detected in 4 wells (Table 2 3 and Figure 2 1). In the subsequent analyses, all values gr eater than zero were considered. Selected Geochemistry of Sites with Excess Nitrogen Gas All measured geochemistry values for this study are provided in Appendix 1. S elected mean geochemical data from both sampling eve nts were compared for sites with above zero (Table 2 3). The h igh est values were detected in sites 5, 9, 10, and 12. The highest values occurred at s ite 9 which had highly reducing anoxic conditions with low mean DO concentrations and negative ORP, low me an nitrate and DOC concentrations and highest measured iron and manganese concentrations . Site 10 had the second highest value, with moderate mean DO concentrations ( 5.67 mg L 1 ), positive ORP , and highest mean nitrate, sulfate , a nd DOC concentrations . Site 12 had anoxic conditions with low mean DO concentrations ( 2.69 mg L 1 ) , substantial mean nitrate concentrations (0.73 mg L 1 ) and lowest mean DOC concentrations. Site 5 was oxic with high mean DO concentration (6.81 mg L 1 ), positive O RP, low mean nitrate concentration, and the second highest mean DOC concentrations (0.93 mg L 1 ).
37 Geospatial Analysis To investigate the role of confinement as a driver of the spatial distribution of denitrification, mean behaviors of spatially averaged g eochemistry for the confined and unconfined regions of the Ichetucknee Springshed were examined. The an alysis revealed mean ORP ( 116 mv) was significantly different from ORP in the unconfined region (116mv, p=0.0 1); iron and manganese concentrations were higher in the confined region than the unconfined region, p=0.09 and p=0.02, respectively. However, parameters which describe conditions that facilitate denitrification (DO) and reactant availability DOC, and nitr ate, had similar means (p > 0.19 ) for the c onfined and unconfined regions (Table 2 4). To gain insights into factors associated with values, mean geochemical (e.g. NO x N, SpC, ORP) and geospatial data (e.g. DRASTIC Index, confining unit thickness) were compared for sites with and without detection (Table 2 5). This comparison indicated sites with were associated with higher mean nitrate concentrations (0.95 mg L 1 , P<0.009), and lower confini ng unit thickness (14.33 ft, P=0.12 ) , higher mean ORP ( 128.70 mv, P =0.13), shallower depth to water (29.81 ft, P=0.13). However there was an insignific ant difference in DOC (P =0.58) , iron (P=0.54, manganese (P=0.79), and DO concentrations (P=0.79 ). Linear regressions between concentrations an d geospatial data were used to test the prediction that aquifer confinement (e.g., DRASTIC index, confinement unit thickness) and proximity to surface features that enhance surface to aquifer connectivity (e.g. to CTDs, Sinks, Swallets) drive observed vari ation in the spatial distribution of denitrification. Linear regression s between denitrification evidence,
38 aquifer confinement, and proximity to CTDs, sinks, swallets, and wetlands w ere conducted for all measured values and values greater than zero. Only values greater than zero showed strong spatial correlation with the number (R 2 = 0.85, p=0.006, Figure 2 4A), area (R 2 = 0.85, p=0.006, Figure 2 4B), and nearest di stance to (R 2 =0.64, p=0.04, Figure 2 4C) wetlands within a 500m buffer, all other relationships were statistically insignificant (p<0.7). Spatially Averaged and Flow Averaged Denitrification Evidence In this study, mean across all spring shed groundwater samples was 0.31 mg L 1 which was similar but less than the 0.51 mg L 1 mean from spring vent samples. Mean [NO3] 0 values for the wells (1.22 mg L 1 ) were similar to prior spring samples (1.18 mg L 1 ). Groundwater [NO3] 0 values ranged from 0.02 to 5.56 mg L 1 (Table 2 2), whereas in spring samples [ Heffernan et al. , 2012] [NO3] 0 values ranged from 0.78 to 1.68 mg L 1 . The residu al fraction of nitrate ( f ) values for well samples ranged between 0.04 and 1.0 with a mean of 0.71 (Table 2 2 ), whereas the reported springs samples ranged from 0.28 to 0.85, with a mean of 0.52 [ Heffernan et al. , 2012] . 15 18 O isotopic enrichment ratio for well samples (1.41, Figure 2 3) was much higher than the enrichment ratio (0.61) calculated from all available spring samples (including spring measur ements collected during this study and samples collected during the study by Martin et al. [in press], or the enrichment ratio (0.99) calculated from a suite of UFA spring vent samples [ Cohen et al. , 2012; Heffernan et al. , 2012] . The linear relationship between 15 N and DO in well samples was insignificant (P=0.6, Figure 2 5A), however the slope obtained for the springs using all availab le data was significant (P =7 E 8, Figure 2 5B).
39 Discussion Dat a from this study illustrated the utility of spatially distributed nitrate isotope and dissolved gas evidence to characterize the magnitude and spatial distribution of denitrification in karst springsheds, provided information about conditions (e.g. redox state and reactant availability) that facilitate denitrification in an eogenetic karst springshed across an aquifer confinement gradient, and examined how land surface aquifer connections may facilitate denitrification. These data show significant heteroge neity in N cycling, geochemistry and redox state. H igh geochemical heterogeneity and low spatial density of available wells for analysis contributed to underpowered statistical results . Statistical comparisons for geochemistry in samples obtained in the c onfined and unconfined region showed powerful (>0.65) and significant (P<0.02) differences in DIC, pH, ORP, and manganese. However, high heterogeneity (i.e. high standard deviation) and low effect size values led to low powered statistical evaluations for all other parameters. Specifically in samples obtained where the UFA is confined, mean DO and nitrate concentrations we re lower and mean DOC concentrations were higher than the unconfined region. This led to the expectation that these parameters would be s ignificantly different between the unconfined and confined aquifer groundwater samples. However, Mann Whitney U test results for these parameters were underpowered and differences were insignificant (P>0.19) (Table 2 4). Parameters which describe redox st ate (e.g. DO and ORP) and reactant availability (e.g. nitrate and DOC concentrations) were expected to be significantly different between sites with and without detection. However, all parameters had low statistical power (<0.3) and only mean nitrate concentrations were significantly
40 different (P<0.01) (Table 2 5). Together these statistical analyses illustrate the extreme heterogeneity in geochemistry and discordance between expected conditions that facilitate denitrification and observed local porewater conditions. To examine how observed geochemical heterogeneity relates to observed springshed scale denitrification evidence, spatially aggregated (i.e. averaged) point measures were used to reconstruct previously observed spring vent (i.e . flow averaged) denitrification evidence from spatially distributed well samples. Utility of Nitrate Isotopes and Dissolved Gas Measurements in Karst Aquifers: The utility of distributed isotopic measures as primary evidence to describe nitrate transform ations in karst aquifers is limited. Spatially distributed denitrification evidence was equivocal using isotopic methods. Denitrification could not be inferred with isotopic data because enrichment was not significant enough to clearly plot outside any of the N source ranges or rule out source mixing with organic waste or atmospheric nitrogen sources (Figure 2 3). This finding implies that mixing of multiple isotopic sources remains the strong est influence on local scale isotopic ratios , perhaps masking the influence of denitrification on isotopic ratios . Interpretation of these ratios depends on th e assumption that measurements we re made along distinc t hydrologic flow paths which are difficult to discern in karst aquifers. In contrast, spring vent isotope m easures unequivocally describe flow averaged springshed denitrification, as evident by 15 18 O covariance across multiple studies (Figure 2 3) [ Katz , 2004; Cohen et al. , 2012; Heffernan et al. , 2012 ] . Dissolved gas measurements are not influenced by local scale isotopic variations in nitrate source or nitrogen transformation in the unsaturated zone before recharge. Dissolved gas measures provide direct evidence of that denitrification occurred in t he
41 UFA. These measures require assumptions about gas solubility and recharge temperatures, but adding noble gases (e.g. Xe) reduces uncertainty. Dissolved gas measurements provided unequivocal evidence of denitrification ; future sampling effo rts should foc us on collecting these data as a direct measure of denitrification . Conditions That Facilitate Denitrification Redox state and reactant availability It was expected local conditions such as aquifer redox state and reactant availability would influence concentrations . In a national study, Burrow et al. , noted that redox conditions, and iron concentrations explain most of the variation in nitrate concentrations. However due to preferential flow paths in karst aquifers, the relationshi ps between redox state, reactant availability, and denitrification S ites with detection were not significantly associated with indicators of redox state such as low DO and low or negative ORP or electron don or availability (e.g. DOC, iron, manganese) (Table 2 5 ). However, t he highest value did occur at s ite 9 where conditions were the most anoxic (lowest DO concentrations, negative ORP) . Significant variations in iron and manganese concentrat ions in the unconfined region where all detection occurred, and high iron and manganese concentrations where the highest value occurred support that solid phase alternate electron donor availability may be facilitating c hemolithoautotrophic denitrification . However, across all measured sites, redox state (e.g. DO concentrations and ORP ) was not a good predictor of the magnitude of denitrification (Table 2 3), the spatial distribution of values (Table 2 5 ) , 15 N enrichment (Figure 2 4A). Often values were high even where DO
42 concentrations were unfavorable for denitrification [ Rivett et al. , 2008] . Together these data support that low redox state facilitates denitrification but is not necessary. There could be several explana tions for this unexpected discordance between redox stat e and observed denitrification facilitate oxygen metab olism to reduce DO. Also, l ocal reduction zones in microbial biofilms on aquifer substrates could create favorable conditions where transformations could occur , t his localized biogeochemistry in biofilms has been documented previously on saturated soil aggregate biofilms [ Sexstone et al. , 1985] , in biofilm reactors [ Kuenen et al. , 198 6; Revsbech et al. , 1989; Biesterfeld et al. , 2003; Satoh et al. , 2004] and sulfidic karst environments [ Macalady et al. , 2006; Gra y and Engel , 2013] . Furthermore, Biesterfeld et al. [ 2003 ] noted significant biofilm denitrification even with available DO in the bulk aqueous phase. Another mechanism that can explain the lack of concordance between redox state and denitrification is g roundwater mixing. Measured DO concentrations and ORP in wells may reflect a mixture of local porewater and higher DO concentrations and ORP from more recently recharged water conveyed though karst preferential flowpaths. Furthermore, measured may reflect contribution from upgradient denitrification. This finding suggests that future studies would benefit from quantifying water sources near wells using end member mixing models, determining the potential role of DOC in reducing DO to concentrations which facilitate denitrification, and examine the potential importance of microsite denitrification through sampling of local biofilms using recently developed sonication sampling [ Ugolini et al. , 2014] .
43 Despite unfavorabl e local DO concentrations and ORP values in some wells, evidence of denitrification was observed where reactant availability (e.g. nitrate and DOC) was high (sites 10 and 5), and deep groundwater upwelling has been documented previously (site 12) [ Champion and Upchurch , 2003a; Moore et al. , 2009, 2010] . Site 12 is located near the springs and thus may have contributed from converging older groundw ater flow paths. This is not surprising given that convergent flow paths near springs may mix water of various ages and sources [ Martin and Gordon , 2000] . Aquifer confinement Aquifer confinement was hypothesized as a potential driver of observed variation in denitrification evidence because confined aquifers have been prev iously shown to be hotspots of denitrification [ Rivett et al. , 2008] and in the springshed the conditions for denitrification are favorable in the confined region with relatively high mean NO x N and DOC , low mean DO concentrations, negative mean ORP, and evidence of available solid phase electron donors evident from high iron, manganese, concentrations (Table 2 4). In spite of this , a distinct progression of geochemical conditions from anoxic to oxic conditions was not observed in groundwater samples taken from wells across the confinement gradient and all evidence of denitrification was observed in the unconfined region. Denitrification evidence examined with proxies fo r aquifer confinement indicate that all sites with occurred below the Cody Escarpment in the unconfined region of the springshed (Figure 2 1B, Table 2 4); sites with were associated with high nitrate concentrations and shallow water table depths, and decreased confining un it
44 thickness (Table 2 5). Furthermore, ORP was significantly lower in the confined region than the unconfined region while mean DO, DOC, and nitrate concentrations, had statistically similar mean behaviors (p > 0.19) for the confined and unconfined groups (Table 2 4). This suggests that geochemical conditions imposed by aquifer confinement that favor denitrification were not driving observed nitrate transformations. In the confined region, horizontal surficial flow paths focus may flow toward depressions, s inks, and wetlands that may be connected to the UFA, possibly creating nitrate transformation hotspots through delivery of missing reactants over small spatial scales that the limited density of wells in this study did not capture. It is possible that ther e is higher spatial heterogeneity in groundwater concentrations in the confined Ichetucknee Springshed due to solute delivery through focused recharge, whereas groundwater in the unconfined region may be more spatially homogeneous due to solute delivery th rough more uniform recharge. To explore this idea, the coefficient of variation (CV) in geochemical parameters was compared between the unconfined and confined region (Table 2 6). Indeed, the variation in samples from the con fined region was higher for man y geoche mical parameters . However, nitrate, iron and manganese concentrations in the unconfined region had larger CV. This finding paired with detection only in the unconfined region and significant variation support that availability of nitrate and solid phase iron and manganese may influence the spatial distribution of denitrification, as suggested by Heffernan et al. [ 2012] . Higher spatial variability in geochemistry in the confined region indicates that sampling effor ts must be denser in order to effectively characterize average conditions over the confined area .
45 Land Surface aquifer c onnections To test the prediction that proximity to flow focusing features in karst terrain influence the spatial distribution and magni tude of denitrification, linear regression between denitrification evidence and proximity to depressions, sinks, and wetlands was evaluated for all measured values and values greater than zero. However, only values greater than zero showed strong correlation with wetland metrics (number of wetlands, wetland area, and nearest distance; Figure 2 4). This implies that UFA wetland connections may contribute to denitrification in the aquifer. Most surveyed w etlands in the region are geographically isolated (i.e. no persistent surface outflows) and coincident with local depressions that may collect water and essential reactants (e.g. DOC) from the surrounding environment. Thus wetlands may represent connection s that provide a key DOC source to an e lectron donor limited aquifer. Alternatively, the denitrification evidence observed in the UFA wells could be indicative of reactions taking place in the saturated zone in wetlands where conditions do not allow excess N 2 to be returned to the atmosphere. However many sites were near wetlands and did not have , suggesting that either wetland proximity is a necessary but not sufficient condition to facilitate nitrate transformations , or that surveyed wet lands are not in a hydrologic flowpath upgradient from the sites that did not show denitrification evidence. The underlying aquifer may be required to have the right conditions (e.g. lower DO, enriched nitrate) and be connected to a wetland delivering reac tants for the wetland to influence the magnitude of measured . Whether wetlands convey missing reactants to the aquifer or the observed occurred in
46 the wetlands and was subsequently conveyed to the aquifer is a subject for further research. Spatially aggregated and flow averaged m easurements Spring vent samples are spatiotemporally averaged (i.e. flow averaged) ; they i ntegrate mean behavior of basin scale geochemistry, reactant availability and recharge across converging hydrologic flowpaths at springs. In contrast, spatially aggre gated measurements integrate local scale point measurements of geochemistry, reactant availability and recharge variation resulting from a particular hydrologic flowpath. With sufficient sampling density, one would expect that spatially aggregated samples could reconstruct observed behaviors deduced from spring measurements. S patially aggregated measurements (0.31 mg L 1 ) and m ean [NO3] 0 (1.22 mg L 1 ) were close to flow averaged spring measurements for (0.51 mg L 1 ) and [NO3] 0 (1.18 mg L 1 ), suggesting good correspondence between the two methods. Inert excess N 2 gas records denitrification and is not influenced by local scale nitrogen cycling. On the other hand, isotopic denitrification evidence was influenced by local s cale variability in N source and cycling. Typical metrics used to evaluate denitrification using isotopic evidence 15 15 N enrichment trends) are highly influenced by local heterogeneity in nitrogen transformation mechanisms such as nitrification. I 15 for the synop tic well sampling were much lower than those measured in the springs (Figure 2 5 B). Furthermore unlike prior spring measurements , the distributed measures of NO x N 15 15 N was not corr elated with decreasing DO concentrations (Figure 2 5 A). These findings underscore
47 how significant heterogeneity in the UFA influences local scale isotopic denitrification evidence. Complex hydrologic flow paths can mix denitrified water from upgradient fl ow paths , where conditions and reactant availability are more favorable , with recently recharged water contributing uncertainty in the source of observed . Low density point sampling renders study findings subject to the confounding effects of this heterogeneity. H igh ge ochemical heterogeneity and low spatial density of available wells for analysis contributed to underpo wered statistical res ults. Thus if characterization of chemistry at the springshed scale is desired, unreasonable numbers of samples may be required to capture all of the variability. This heterogeneity presents challenges in reassembling springshed N dynamics from spatially d istributed point measurements. While distributed measures taken through time may provide insights into springshed behaviors (i.e. change detection), spring vent measures have the advantage that they are inherently spatiotemporally averaged and provide clea r relationships that are not influenced by local heterogeneity. The contrasting behavior of measurements at differen t scales further underscores the challenges of interpreting spatially distributed measurements in karst aquifers. Study Summary Many concep tual models for interpreting geochemical data rely on spatial or temporal progressions in geochemistry that are often difficult to define in karst aquifers due to variation and connectivity in preferential flow paths. Significant heterogeneity was observed in groundwater redox condition, reactant concentrations, and the distribution of land surface aquifer connections. While spatially distributed measures provided insights into local scale evidence of denitrification, extension of these behaviors to
48 larger spatial scales was plagued by uncertainty induced by process variability (i.e. different N transformation mechanisms), local recharge and geochemical variation, and poorly described upgradient and downgradient hydrologic flow paths. Denitrification was de tected even where local aquifer conditions are presumed to be unfavorable. Thus, assumptions of local denitrification potential based upon porewater chemistry alone may neglect the potentially significant influence of biofilm nitrate transformations or upg radient nitrate transformations. Direct sampling of nitrate transformation in aquifer material biofilms may provide insights into nitrate transformation where aquifer conditions are considered unfavorable for nitrate transformations. Aquifer confinement an d k arst flow focusing mechanisms (e.g. sinks) were hypothesized to be important drivers of spatial variability in denitrification . Aquifer confinement can lead to conditions that facilitate nitrate transformations (e.g. low DO) and flow focusing mechanisms can quickly route reactants into the aquifer (e.g. DOC). N either aquifer confinement, sinks, swallets, nor closed topographic depressions were correlated with denitrification proxies (i.e. excess N 2 gas). Results support that wetland aquifer connections m ay be associated with denitrification proxies. However, this relationship only held for sites with detection, indicating that wetland connectivity to local aquifer conditions may not be simply predicted by geographic proximity in karst aq uifers where extensive subsurface preferential flow paths exist. The potential wetland aquifer denitrification relationship observed in this study warrants further investigation.
49 Comparison of spatially distributed and spring vent sampling methods support that high density point measurements a re required to reproduce spring vent signals, especially if characterization of major anions and cations and isotopic evidence is desired. However, dissolved excess N 2 gas and reconstructed initial nitrate concentrati on measurements for both methods had surprisingly good concordance given the limited available spatial density for sampling. Thus, spatially distributed measurements provide information about spatial heterogeneity in geochemistry, local scale N cycling and attenuation, and when spatially aggregated, reproduce direct denitrification evidence reasonably. In summary, data support that denitrification occurred only in the unconfined region where the UFA is closely connected to the surface, where nitrate reactan t availability is high, and that low redox state (e.g. negative ORP, and low DO) was not a good predictor for denitrification. Isotopic denitrification evidence was inconclusive and highly influenced by local variations in N loading and cycling. Direct evi dence of excess dissolved N 2 gas supports that denitrification occurred throughout the Ichetucknee Springshed, was highly spatially heterogeneous, and has a suggestive correlation with wetlands. While spatially distributed dissolved gas measures were heter ogeneous, point scale aggregation resulted in reasonable correspondence between reconstructed initial nitrate concentrations and excess N 2 measured at springs.
50 Figure 2 1. Ichetucknee Springshed Sampling Site map . A) Site location in Florida with Flori dan aquifer confinement. B) Ichetucknee Springshed showing Gulf Coast Lowlands, Northern Highlands Province, Cody Escarpment, Rose, Clay Hole, and Cannon Sinks, streams, Ichetucknee Springs Group, and sampling site locations.
51 Figure 2 2 . Spatial data for Ichetucknee Springs and Cody Escarpment. A) S patially distributed DRASTIC Index sampling sites. B) L ocations of potential land surface aquifer connections; wetlands, closed topographic depressions, sinks and swallets.
52 Figure 2 3 . Measured isotopes of nitrate, showing nitrogen source ranges for synthetic (shaded), soil, and animal waste nitrogen sources, and linear regression for data collected in the springshed (black) and spring vents (gray circles) as part of this study, samples collected by Mart in et al., [ 2016] (gray triangles), as 18 O 15 15 18 O slopes for the springshed (solid line), and all springs (gray line).
53 Table 2 1. Ichetuckn ee Spring Basin w ell information. Site number Latitude Longitude Elevation (ft) Well depth (ft) Depth to water (ft) Open interval range (ft) 1 82.59 30.18 197 836 95 183 203 2 82.61 30.10 140 133 103 113 133 3 82.64 30.18 141 836 95 689 836 4 82.70 30.07 60 84 15 50 84 5 82.71 29.98 77 65 49 45 65 6 82.70 30.07 76 150 31 120 150 7 82.71 30.19 161 180 106 140 180 8 82.75 29.97 45 36 21 16 36 9 82.75 30.02 58 53.1 28 28 53 10 82.72 30.05 54 35 22 33 35 11 82.81 30.21 199 165 143 128 165 12 82.75 29.98 50 34 26 14 34 13 82.76 29.96 48 39 26 19 39 14 82.79 29.96 50 37 21 17 37 15 82.77 29.98 55 41 33 21 41 16 82.75 30.02 65 49 36 29 49 Table 2 2 . Estimates of excess N 2 (N 2 den ), initial NO 3 N concentrations [ , resi dual fraction of original NO3 (f) observed at the sampling site, and estimated 15 15 N NO 3 ] init ) based on range of fitted global enrichment factors. All nitrate, nitrogen, and excess air concentrations are in units of mg L 1 , f is unitless, temperature is in Celsius, and isotope value units are in per mil. Date Site Excess Air Recharge Temp. NO x N [N 2 ] den [NO 3 ] 0 f 15 N NO 3 ] obs 5 N NO 3 ] init 3/18/ 14 1 1.26 23.52 0.02 1.10 0.02 1.00 5.26 5.26 3/18/ 14 3 0.00 23.87 0.03 0.09 0.0 3 1.00 3.51 3.51 3/7/ 14 5 0.35 23.41 0.38 0.60 1.59 0.24 6.75 4.49 3/9/ 14 8 0.80 22.70 0.39 0.20 0.80 0.49 8.72 7.60 3/14/ 15 9 1.81 22.22 0.12 1.37 2.86 0.04 4.01 0.91 3/14/ 14 10 0.94 22.75 3.68 0.94 5.56 0.66 4.47 3.82 3/9/ 14 12 1.52 23.31 0.78 0.55 1.88 0.41 3.43 2.04 3/7/ 14 13 1.28 20.86 0.04 0.63 0.04 1.00 7.61 7.61 3/9/ 14 14 1.62 20.49 0.37 0.08 0.53 0.70 2.89 2.34 3/7/ 14 15 0.45 22.63 0.06 0.00 0.06 1.00 2.49 2.49 3/14/ 14 16 0.15 24.52 0.03 0.11 0.03 1.00 1.73 1.73
54 Table 2 3 . Phase I an d Phase II g eochemical parameters for the four sites with the highest excess N 2 gas ([N 2 ] den ) values. Site Sample DO ( mg L 1 ) ORP (mv) SO 4 2 ( mg L 1 ) NO x N ( mg L 1 ) DOC ( mg L 1 ) Fe 2+ ( mg L 1 ) Mn 2+ ( mg L 1 ) 5 Phase I 6.81 158.60 0.32 0.33 1.23 0.004 0.0 00 5 Phase II 6.81 262.40 0.48 0.38 0.63 0.012 0.001 5 Mean 6.81 210.50 0.40 0.35 0.93 0.01 0.00 9 Phase I 0.15 177.30 8.16 0.07 0.26 0.463 0.038 9 Phase II 0.53 141.00 8.53 0.12 0.70 0.526 0.035 9 Mean 0.34 159.15 8.34 0.10 0.48 0.49 0.04 10 Phas e I 5.27 41.60 28.93 4.39 1.17 0.003 0.001 10 Phase II 6.06 119.10 27.42 3.68 1.17 0.004 0.001 10 Mean 5.67 80.35 28.18 4.03 1.17 0.00 0.00 12 Phase I 2.37 118.40 5.12 0.69 0.17 0.003 0.002 12 Phase II 3.01 176.50 5.08 0.78 0.54 0.001 0.000 12 Mean 2. 69 147.45 5.10 0.73 0.36 0.00 0.00 Table 2 4 . Confined and unconfined aquifer Mann Whitney U test results for selected geochemical and geospatial data, showing m ean ( Mann Whitney U statistic, probability (P), and statistical p ower . The null hypothesis is that mean values for sites in the confined and unconfined regions are the same. Statistical power analysis assumed a significance value of 0.05. Confined n=7 Unconfined n=20 Effect Size U P Power DIC (mg L 1 ) 25.05 12.96 44.88 9.21 1.48 12.00 0.00 0.88 pH 7.76 0.51 7.08 0.56 1.10 123.00 0.00 0.65 K + (mg L 1 ) 3.83 8.20 0.41 0.36 0.81 122.00 0.00 0.41 SPC (ms cm 1 ) 294 195 493 178 1.00 20.50 0.01 0.57 ORP (mv) 116 218 100 116 1.25 23.00 0.01 0.76 Mn 2+ (mg L 1 ) 0.02 0.03 0.00 0.01 1.11 111.00 0.02 0.65 PO 4 3 ( g L 1 ) 0.03 0.02 0.06 0.04 0.83 32.00 0.04 0.42 Fe 2+ (mg L 1 ) 0.10 0.11 0.06 0.16 0.21 101.00 0.09 0.02 DOC (mg L 1 ) 1.16 1.01 0.52 0.34 1.01 94.00 0.19 0.58 DO (mg L 1 ) 2.77 3.60 4.83 2.64 0.69 47.00 0.21 0.31 NO x N (mg L 1 ) 0.39 0.54 0.76 1.26 0.33 53.00 0.36 0.10
55 Table 2 5 . Excess nitrogen gas Mann Whitney U test results for selected geochemical and geospatial data, showing m ean ( U statistic, probability (P), and statistical p ower . The null hypothesis is that mean values for sites with and without [N 2 ] den are the same. Statistical power De tected [N 2 ] den n =6 No [N 2 ] den n=5 Effect Size U P Power NO x N (mg L 1 ) 0.95 1.35 0.04 0.01 0.86 30 0.00 0.23 Confining unit t hickness (ft) 14.33 16.28 46.60 50.00 0.86 6 0.12 0.23 pH 6.82 0.18 7.48 0.86 1.01 6 0.13 0.30 ORP (mv) 128. 70 139.84 23.76 257.06 0.74 24 0.13 0.19 Depth to w ater (ft) 29.81 10.60 58.55 37.13 1.00 6 0.13 0.29 DOC (mg L 1 ) 0.66 0.27 1.28 0.97 0.85 8 0.25 0.23 Fe 2+ (mg L 1 ) 0.09 0.21 0.06 0.13 0.17 11 0.54 0.06 Mn 2+ (mg L 1 ) 0.01 0.01 0.02 0.03 0.46 13 0.79 0.10 DO (mg L 1 ) 4.47 2.60 3.82 3.31 0.23 17 0.79 0.06
56 Figure 2 4 . Linear regression for all sites (gray line) and sites with observed excess N 2 gas (black line) for three wetland metrics showing R 2 correlation coefficient and p value for all sites (left) and sites with excess N 2 (right). A) Linear regression for wetland area within 0.5 km in (km 2 ). B) Linear regression for numbe r of wetlands within 0.5 km. C) Linear regression for di stance to nearest wetland in km.
57 Figure 2 5 . Spr ingshed and spring vent linear regressions for 15 N , dissolved oxygen and enrichment ratios. A) R 15 N and dissolved oxygen for springs (gray line) and spatially distributed sites (black line) . B) The 15 N enrichment factors computed using linear regression of spatially distributed sites (black line) and all spring measurements (gray line). Table 2 6 . Comparison of coefficients of variation (CV) in measured geochemistry between the confined and unconfined regions of the Ichetucknee Springs Basin . Parameter Confined CV Unconfined CV Difference pH 0.07 0.08 0.01 DO (mg L 1 ) 1.30 0.55 0.75 ORP (mv) 1.87 1.16 0.71 SPC (ms cm 1 ) 0.66 0.36 0.30 Na + (mg L 1 ) 0.79 0.40 0.38 K + (mg L 1 ) 2.14 0.88 1.26 Mg 2+ (mg L 1 ) 0.73 0.82 0.08 F (mg L 1 ) 0.96 1.16 0.20 Cl (mg L 1 ) 0.58 0.37 0.21 SO 4 2 (mg L 1 ) 2.01 1.01 1.00
58 Table 2 6 (continued) Parameter Confined CV Unconfined CV Difference NO x N (mg L 1 ) 1.39 1.67 0.28 DOC (mg L 1 ) 0.87 0.66 0.21 DIC (mg L 1 ) 0.52 0.21 0.31 Fe 2+ (mg L 1 ) 1.19 2.45 1.26 NH 3 (mg L 1 ) 1.29 0 .94 0.35 Mn 2+ (mg L 1 ) 1.08 2.52 1.44 Sr 2+ (mg L 1 ) 1.33 0.60 0.72
59 CHAPTER 3 NITRATE TRANSFORMATION MECHANISMS AND RATES IN AN UNCONFINED EOGENETIC KARST AQUIFER ACROSS A REDOX GRADIENT Background Local biogeochemical conditions control N transformat ions in the subsurface environment [ Pinay et al. , 1993; Vidon and Hill , 2004] . N transformations requir e electron donors, and favorable environmental conditions such as low dissolved oxygen (DO), pH, temperature, and micronutrient availability [ Rivett et al. , 2008] . N transformation pathways, include assimilatory nitrate reduction to ammonia (ANRA), anaerobic ammonium oxidation (ANAMMOX), dissimilatory nitrate reduction to ammonium (DNRA) [ Tiedje , 1988] , and denitrification [ Davidson et al. , 2003; Weber et al. , 2006] . ANRA is used by heterotrophic microbes to convert nitrate into bioavailable ammonium [ Guerrero et al. , 1981; Guerrero , 1985] . This mechanism contributes to nitrogen cycling rather than removal since die off of the microbial population c ould lead to remobilization of N. ANAMMOX bacteria convert ammonium and nitrite to N 2 gas under anoxic conditions [ Kartal et al. , 2011] . Microbes capable of ANAMMOX are pervasive through terrestrial and marine ecosystems [ Humbert et al. , 2010] and this process was recently found to be significant in a freshwater aquifer study [ Smith et al. , 2015] . DNRA converts nitrate to ammonium. DNRA differs from denitrification in that its end product, ammonium, is less mobile due to potential adsorption onto aquifer substrate , highly bioavailable and is retained in th e environment [ Bu rgin and Hamilton , 2007; Rivett et al. , 2008; Giblin et al. , 2013] . Denitrification mitigates nitrogen enrichment by reducing nitrate to nitrogen gas and returning it to the atmosphere [ Ottley
60 et al. , 1997; BÃ¶hlke , 2002; Rivett et al. , 2008; Groffman et al. , 2009] and is a significant part of global nitrogen cycli ng [ Seitzinger et al. , 2006] . Denitrification rates have high uncertainty, due to issues in quantifying denitrification byproducts, spatiotemporal variability in nitrate and denitrification reactants, and multiple denitrification pathways [ Davidson and Seitzinger , 2006; Groffman et al. , 2009; Groffman , 2012] . Thus improved quantification of denitrification rates and controlling factors in the laboratory and field are needed [ Davidson and Seitzinger , 2006; Rivett et al. , 2008] . Denitrification and DNRA both occur under anoxic conditions through either organic matter oxidation (heterotrophic), or biotic or abiotic iron (Fe 2+ ), sulf ide (S 2 ), or manganese (Mn 2+ ) oxidation (chemolithoautotrophic) [ Korom , 1992; Brunet and Garcia Gil , 1 996; Otte et al. , 1999] . In the chemolithoautotrophic pathway, the end product can either be N 2 derived from denitrification, or ammonium through DNRA [ Brunet and Garcia Gil , 1996] . Metal oxidation can provide viable alternative electron donors [ Kolle et al. , 1985; Robertson et al. , 1996; Tesoriero et al. , 2000] . This may be important in low carbon aquifer systems like the UFA [ Cohen et al. , 2007; Heffernan et al. , 2012] . Microbially mediated N transformations, ANRA, and subsequent N release from biomass fermentation may lead to nutrient spiraling in aqui fers. This concept has been previously applied to river systems to account for changes in nutrient cycling as nutrients are carried downgradient [ Newbold et al. , 1981] . The nutrient spiraling analysis framework employs quantification of uptake lengths, average distance traveled by dissolved nutrients in the water column before uptake, and which describe the distance over which the nutrient exists in the abiotic and biotic phase. Ammonium and nitrate uptake lengths are a widely used parameter for
61 quantifying nutrient cycling in streams [ Covino et al. , 2010; Mulholland et al. , 2016] . However application of this concept in aquifers is limited; nitrogen studies have typically focused on nitrate transport and attenuation. Uptake length is a metric of nutrient use efficiency, quantifying nutrient uptake relative to supply. This information is critical for det ermining the assimilative capacity of aquifers and managing nitrogen loads. Study Objectives In this study, aquifer and microbial responses to nutrient additions and rates of nitrate transformation were directly measured in the UFA, an eogenetic karst aqu ifer in the SE USA (Figure 3 1). Potential nitrate reduction mechanisms in the UFA were measured across a redox gradient (e.g. ANRA, ANAMMOX, DNRA, and denitrification). Nitrate reduction processes DNRA and ANRA retain nitrogen in the environment whereas A NAMMOX and denitrification remove it. Thus quantifying nitrate transformation mechanisms and rates is essential to long term management of nitrogen enrichment. Prior research has suggested that denitrification is the most significant nitrate transformatio n mechanism where it occurs [ Heffernan et al. , 2012] ; dissimilatory reduction of nitrate to ammonium (DNRA) is negligible in the UFA and other aquifers due to l ow ammonium and dissolved organic carbon (DOC) concentrations [ Katz et al. , 2004; Brown et al. , 2008; Rivett et al. , 2008; Heffernan et al. , 2012] ; and nitrogen transformation rates are limited by carbon availability [ Katz , 1999, 2004; Katz et al. , 2001; Cohen et al. , 2007; Heffernan et al. , 2012] . This study evaluates these assumptions derived from spring vent measurements using in situ ni trate transformation measures.
62 Estimates of denitrification have previously been obtained in the study area at two spatiotemporal scales using proxies for denitrification (isotopes of nitrate) and direct evidence of denitrification such as excess nitrogen gas (N 2 ); synoptic distributed UFA monitoring well measurements (e.g. chapter 2) and spatiotemporally averaged measurements at spring vents emanating from the UFA [ Katz and Griffin , 2008; Heffernan et al. , 2012] . In this study, the fate of nitrogen was examined using single well push pull tracer tests (PPTTs) combined w ith microbial genetic sampling of key genes in N transformation pathways in order to more accurately quantify aquifer response to nitrate additions, measure in situ nitrogen transformation mechanisms and rates, and examine potential nutrient cycling in the UFA. Methods Study Site The study area is the Ichetucknee Springshed which discharges from the UFA in North Florida (Figure 3 1). The springshed is approximately 960 km 2 as estimated by water year 2003 potentiometric surfaces [ Sepulveda et al. , 2006] and further validated by geochemical methods [ Champion and Upchurch , 2003b] . While a confinement gradient exists across the spri ngshed , the sampling sites were located where the UFA is unconfined. In the unconfined region w here limestone is not exposed at the surface, a thin layer of well to excessively drained sands is present [ Houston , 1965; Ceryak et al. , 1983] . The Ocala L ime stone, which makes up the UFA locally, is highly permeable and transmissive with connected voids, conduits, and caves as a result of dissolution of the limestone matrix [ Bush and Johnson , 1988] . This leads to heterogeneous hydrologic flow paths and high surface area where biological reactions may occur. Groundwater flow is predominantly from north to south within the area and the aquifer is primarily
63 drained by the Ichetu cknee Springs . Detailed geomorphic and hydrologic descriptions of Ichetucknee Springs and local geology have been previously reported [ Rosenau et al. , 1977; Scott et al. , 2 004] . Two wells from the Suwanee River Water Management District (SRWMD) monitoring well network in the springshed were selected for the tracer tests and microbial analyses; MW14 (oxic) and MW15 (anoxic). These wells were located 40 0 m from each other (Fi gure 3 1C ) and located on protected state forest lands in the unconfined region of the springshed. Both wells have the same diameter (0.10 m), similar surface elevations ( 20.12 and 17.68 m , respectively), dep ths (14.94 m and 16.15 m ), ope n intervals in the upper Ocala L imestone ( 6.10 m and 7.62 m ) , and water levels (9.144 m and 9.45 m above mean sea level) . Table 3 1 provides measured background geochemistry of the two wells. Both have low mean nitrate concentrations (< 0.12 mg L 1 detection limit) and low mean DOC concentrations (0.34 and 0.48 mg L 1 , respectively). The oxic site has 7.75 mg L 1 mean DO concentrations and highly positive ORP. The anoxic site has mean DO concentrations equal to 0.34 mg L 1 , negative ORP values, and a strong hydrogen sulfide gas smell was observed during sampling. These sites thus represent springshed end members of redox state and present an ideal test location for examining nitrate transformat ion rates and mechanisms in unconfined eogenetic karst aquifers across a redox grad ient. Push Pull Tracer Tests Single well Push Pull Tracer Tests (PPTTs) are widely used to understand nitrate reaction kinetics [ Isto k et al. , 1997; Schroth et al. , 1998, 2000; North et al. , 2004; Kim et al. , 2005, 2011] . In a PPTT, reactive and conservative tracers are injected, reactions occur for a specific time, and continuous measurements of recovered reactive and
64 conservative tra cers provide solute breakthrough curves (BTC). The BTC analysis evaluates nitrate transformation rates that occurred while the solution was resident in the aquifer. For each PPTT, an injectate solution was created using local groundwater combined with rea ctive and conservative tracers. At each site, one PPTT experiment examined aquifer response with nitrate only tracer addition (NO 3 only) and a subsequent PPTT examined response with nitrate and fumarate addition (NO 3 DOC) with fumarate as a DOC source. For each solution, lab grade potassium nitrate, potassium chloride, and sodium fumarate were measured and mixed with 250ml of twice deionized water to create a concentrated solution for mixing with local groundwater. A 1.22 m interval of open borehole was is olated using an inflatable borehole packer to reduce tracer solution volume required for the PPTT. Vertical distributions of nitrate and nitrite (NO x N) and field parameters (DO, temperature, pH, and conductivity) were evaluated for each 1.22 m interval. N O x N concentrations from grab samples were <0.12 mg L 1 . Field parameters measured by a 556 MPS multi probe (YSI, USA) continuously lowered down the borehole and NO x N grab sample concentrations at 1.22 m intervals were vertically homogenous in each well ( <10% variation), so the packer was placed in the center of the open interval in each well, at 13.41 m depth. The equivalent of six interval volumes (40 L) of local groundwater were extracted and split into two carboys with one way air release vales. One ca rboy had pre prepared concentrated injectate solution to be diluted with groundwater, and the other carboy was used to hold groundwater to push the injected solution further into the formation. Once
65 each sealed carboy was filled with site water, nitrogen g as was introduced through an aquarium bubbler to displace dissolved oxygen and mix the solution. Starting with the test solution, the water from the two carboys were gravity injected into the open interval. A summary of the injectate solution concentration s and resting time ( t* ) values for each PPTT experiment is provided with experiment results in Table 3 3. A calibrated Submersible Ultraviolet Nitrate Analyzer probe (SUNAV2; Satlantic, USA) and 6920V2 conductivity and chloride probe (YSI, USA) was used to measure background concentrations before tracer injection and monitor tracer recovery during the experiment. However due to potential interference between fumarate and nitrate in the SUNAV2 sensor readings, grab samples were collected, analyzed in the l ab and used for all BTC analyses. Injected solution was extracted using a Rediflo2 pump (Grundfos, USA) at a steady rate verified periodically throughout each experiment using a graduated cylinder and stopwatch. Grab samples were collected, filtered to 0. acid, and immediately placed on ice to prevent further microbial reactions. Nitrate and nitrite nitrogen (NO x N) in grab samples was measured using EPA method 353.2, ammonium (NH 4 N) using EPA method 350.1, and chloride using EPA Method 325.2 in the University of Florida (UF) Analytical Research Lab. Fumarate was measured using a Shimadzu TOC 5000A total organic carbon analyzer after sparging for 2 minutes with carbon free air to remove inorganic C in the UF Department of Geol ogical Sciences Organic Geochemistry Analytical Lab. Three to five 60 ÂµL samples were measured and the mean of the injections is the reported fumarate value.
66 Break through curve (BTC) analysis was used to determine zeroth order nitrogen loss rates [mass v ol 1 time 1 ] using the following relation [ Istok et al. , 1997] :, (3 1) where is mass of reactive tracer injected, is recovered reactive tracer mass (area under the recovered reactive tracer BTC), is conservative tracer mass recovery fraction (mass of conservative tracer injected divided by the mass of the conse rvative tracer recovered), is injected solution volume, and is time between injection and midpoint of mass recovery of the tracer test pull phase. Zero order denitrification rates (i.e. rates independent of nitrate concentration) have been reported for NO 3 N concentrations >1mg N /l [ Morris et al. , 1988; Smith and Duff , 1988; Korom , 1992] , suggesting that the supply of electron donors controls denitrification rates at these NO 3 N concentrations. Background analyte concentrations were measured by grab sampling prior to PPTTs, and nitrate concentrations were confirmed before tracer injection at each well using the SUNAV 2 probe. Background concentrations were subtracted from all measured values before using BTC integration to estimate mass recovery. For the tracer analysis and computation of the zero order NO 3 N loss rate, injected tracer mass for chloride, NO x N, and fum arate were calculated based upon the mean of triplicate samples of injectate solution (Table 3 3). Loss rates were computed using the difference between injected NO x N and recovered NO x N + NH 4 N. Reported lab accuracy for each analyte (10%) was used to p ropagate uncertainty through BTC integration and rate calculatio ns using Monte Carlo analysis. For each
67 analyte injection concentration and BTC grab sample concentration, a random error from a normal distribution with zero mean and standard deviation equal to the lab reported accuracy was added and Monte Carlo replicates were generated until the ensemble mean of cumulative mass recovery estimate for each measurement time converged. Reported values are the Monte Carlo ensemble mean and standard deviation for injected mass, recovered mass, and resulting zero order rate calculations. Microbial Analyses Microbial DNA and RNA were analyzed in groundwater samples to assess baseline and PPTT gene numbers, and gene expression using qPCR (quantification of DNA) and R T qPCR (quantification of mRNA) . Previous studies have looked at changes in microbial community structure in PPTTs [ Kim et al. , 2011] . Here, the focus was quantifying changes in N transformation functional gene copy numbers and expression including: total 16S rRNA genes (measure of abundance ), DNRA gene ( nrfA ) [ Song et al. , 2014] , ANAMMOX gene ( hzoA) [ Sun et al. , 2014] , and denitrification genes ( nirK and nirS ) [ Henry et al. , 2004; Kandeler et al. , 2006] . collected volume through the sampler (at least 5 L) was measured. Samplers wer e dried and Mobio, USA). DNA and RNA sample filters were sealed and kept at 4 Â°C for preservation and immediate analysis. DNA and RNA were extracted from filters using Powerw aterÂ® Sterivex DNA Isolation Kit (Mobio, USA) [ Hollibaugh et al. , 2014] . RNA was converted to cDNA using cDNA synthesis kits ( Quanta BioSciences , USA). Table
68 3 2 lists all primers of N related functional genes. Duplicate sample and standard DNAs time PCR system (Applied Biosystems, USA) [ Bae et al. , 2014] . Abundance of all genes was calculated as copies of genes per liter of water sample. Standard curves were generated by serial dilution (10 3 to 10 8 ) of reference plasmids containing target nirK / nirS , nrfA , hzoA gene fragments as t he insert to the plasmids extracted from E. coli hosts using TOPO TA cloning kits (life technologies, Thermo Fisher Scientific Inc., USA). Concentrations and purity of plasmid DNAs were measured with absorbance from 230 nm to 320nm using Eppendorf biophoto meter plus (VWR International LLC, USA), and plasmid DNAs with ratio of A260/A280 equal to 1.7 2.0 were used for qPCR analysis. DNA data evaluated total numbers of microbes (approximated by 16S rRNA copy numbers) and specific genes related to subsurface pr ocesses including DNRA ( nrfA ), denitrification ( nirK and nirS ). Functional gene responses to each tracer addition relative to baseline were evaluated by dividing the number of genes by 16S rRNA gene copy numbers (normalization). Normalized gene expression ( c DNA/16S rRNA) data examine gene expression for active processes within sampled water and assess aquifer response to nutrient addition. Gene expression rates were computed by dividing the functional gene mRNA copy numbers (as cDNA) by the total gene copy numbers (i.e. nrfA c DNA/ nrfA DNA). Microbial samplers collected microbes only from the water column so genetic information about potential nitrate transformations in aquifer biofilms was unavailable. However, pore water biomass N (PBN) in the recovered in jectate solution was
69 estimated to evaluate ANRA in the pore water and estimate the proportion of NO 3 N loss explained by PBN changes. Microbial population sizes were estimated using 16S rRNA copy numbers, approximately 3.6 copies of per cell [ Klappenbach et al. , 2001; Schippers et al. , 2005] 3 ) [ Phillips et al. , 2010] and carbon content of bacterial biomass (5.6 x 10 13 g C/ 3 ) [ Bratbak and Dundas , 1984; Bratbak , 1985] were used to estimate carbon mass. Using ranges for N:C weight ratios for freshwater microbes (0.17 to 0.22) [ Fagerbakke et al. , 1996] and solving for N mass in terms of carbon mass yields a N:C molar conversion factor of 0.198 to 0.257. N biomass per liter was estimated by multiplying the estimated carbon mass by this factor. Since DNA sampling occurred throughout the pull phase, genetic data were assumed to represent the entire volume of water extracted. Total PBN est imates were obtained by extrapolating values for N biomass per liter to total number of liters extracted during the pull phase. Results and Discussion Push Pull Tracer Tests Background NO 3 N concentrations were low at both wells, below practical quantifica tion limit for grab sample analysis methods (0.12 mg L 1 ) but within NO 3 N precision and accuracy range for the SUNAV2 probe, 0.004 mg L 1 and 0.028 mg L 1 respectively. Background NO 3 N, chloride, and DOC concentrations were low and similar for the anoxi c and oxic site (Table 3 1). The resting time (t*), percent recovery, measured values and uncertainty for mass injected, mass extracted, and zero order NO 3 N loss rates are provided for each PPTT (Table 3 3). Cumulative mass recovery data (Figure 3 2, Tab le 3 3) show that for all PPTTs, conservative tracer recovery ranged from 98 to 101%. Values slightly above
70 100% were due to reported uncertainty in background and grab sample measurements. Cumulative NO x N tracer recovery ranged between 84 to 93%, indicat ing that nitrate transformations occurred in the aquifer under both anoxic and oxic conditions, with and without fumarate addition. Fumarate tracer recovery percentages were 84% and 58% at the oxic and anoxic sites, respectively. Unrecovered molar mass rat ios of fumarate and NO3 N (C:N) were similar for both sites; 2.08 and 2.20 for the oxic and anoxic sites, respectively. This suggests some similar mechanisms in carbon utilization for both sites. The 2:1 ratio supports fermentive DNRA reaction stoichiometr y [ Robertson et al. , 1996] . However, low ammonium recover y with carbon addition implies that the DNRA reaction may be coupled to direct ammonium removal via nitrifi cation and/or direct microbial assimilation and supports other N transformation mechanisms such as ANRA and denitrification may contribute to some observed carbon losses. The oxic site ha d dissolved oxygen near saturation, positive oxidation reduction pote ntial, and relatively low background nitrate, DOC, iron, and manganese (Table 3 1). These conditions would suggest low potential for local denitrification because high DO impedes denitrification [ Rivett et al. , 2008] and low C:N ratios impedes heterotrophic DNRA [ Kelso et al. , 1997] . However, cumulative mass recovery show e d only 94% of injected NO 3 N was recovered (Table 3 3, Figure 3 2). At the oxic site, NO 3 only injectate solution resulted in similar NO 3 N loss rates for NO 3 only and NO 3 DOC (Figure 3 2A, Figure 3 2A ). Measured zero order NO 3 N loss rates ranged from 0. 03(+/ 0.01) for NO 3 only to 0.05(+/ 0.01) mmol L 1 hr 1 for NO 3 DOC (Table 3 3). Zero order fumarate reduction rates were 0.81 mmol L 1 hr 1 (Table 3 3). A significant NH 4 N increa se was not measured (Figure 3 2D , Table 3 3), suggesting net ammonium
71 pro duction is low. However, gross ammonium production could be significant and masked by microbial assimilation . Observed nitrate losses were unexpected in the oxic zone; however Bengtsson and Annadotter  also noted significant denitrification in aquife r substrate microcosms under aerobic conditions. The anoxic site has low dissolved oxygen, negative oxidation reduction potential, high iron, high manganese, and evidence of hyd rogen sulfide gas (Table 3 1). This evidence supports a local environment wher e heterotrophic and chemolithoautotrophic DNRA and denitrification may be favored. Furthermore, the presence of high iron and manganese and hydrogen sulfide gas indicates that geochemistry near the site favors reactions lower on the redox chain [ Rivett et al. , 2008] . Cumulative mass recovery for all experiments indicate 84 92% of injected NO 3 N was recovered (Table 3 3). With NO 3 only, nitrate loss rates were 0.03(+/ 0.01) mmol L 1 hr 1 and 14 mmol of injected NO 3 N was converted to NH 4 N (Figure 3 2C , Table 3 3). The addition of NO 3 DOC quadrupled the zero order NO 3 N loss rate to 0.13 (+/ 0.02) ( Figure 3 2D , Table 3 3), indicating that other nitrate transformation mechanisms such as denitrification may have been stimulated by DOC. At the anoxic site, zero order fumarate reduction rates were 0.78 mmol L 1 hr 1 (Table 3 3). Microbial Analyses Microb ial DNA and RNA were sampled in groundwater samples taken from each site before the PPTT and during the pull phase of each PPTT. Microbial data for both sites indicate changes from baseline conditions in response to PPTT injectate solutions (Figure 3 3). G enes for denitrification ( nirK and nirS ; Figure 3 3C and 3E ) and genes for DNRA were detected ( nrfA , Figure 3 3E ). However, genes used in ANAMMOX nitrate transformation pathway ( hzoA ) were not detected in any samples. These data suggest
72 that denitrificati on, ANRA, and DNRA were the likely significant mechanisms in observed NO 3 N losses. At the oxic site, microbial samplers for the NO 3 only PPTT were contaminated by filter failure and were excluded from analysis. However, baseline data and microbial respons e to NO 3 DOC were measured. Microbial evidence shows baseline nrfA expression activity, suggesting under background conditions the aquifer is primed for DNRA. With NO 3 DOC, microbial 16S rRNA gene copies (abundance) increased 30 fold relative to baseline ( Figure 3 3 A ) in response to carbon availability. However, DNRA was not a dominant process since nr f A copy numbers and nrfA expression and expression rate did not increase and significant excess NH 4 N was not recovered. Normalized nrfA copies were close to baseline (Fig ure 3 3B ), and nfrA expression (cDNA/16S rRNA gene; Figure 3 3B ) and expression rate ( cDNA/DNA) decreased (Figure 3 3B ). Denitrification competes with DNRA, which may explain the decrease in nrfA expression rate with NO 3 DOC. At the oxic site with NO 3 DOC, normalized nirK copies were similar to baseline conditions (approxim ately 1% difference, Figure 3 3C ) and normalized nirS act ivity increased 67% (Figure 3 3E ) relative to the increase in 16S rRNA gene copies suggesting that denitrification g ene activity increased. However, nirK and nirS expression (cDNA/16S rRNA) was not detected so the genes were not sufficiently active to be detected in the sampled water, further supporting the hypothesis that observed nitrate losses occurred in biofilms on aquifer material which has been found to be a significant location for transformations [ Holm et al. , 1992; Hancock et al. , 2005; Macalady et al. , 2006] even under high DO conditions [ Chang et al. , 2006] .
73 At the anoxic site with NO 3 only, the number of 16S rRNA gene copies decreased 36 fold c ompared to baseline (Figure 3 3A ). Approximately 5% of the injected nitrate was converte d into ammonium, thus DNRA represents about 58% of the unrecovered nitrate mass. The positive response in DNRA microbial indicators ( nrfA cDNA and DNA ; Figure 3 3B, Figure 3 3D ), and decrease for denitrification microbial indicators ( nirS , nirK cDNA and DN A; Fig.3C , Fi gure 3 3E ) with NO 3 only tracer addition suggest that most DNRA groups were chemolithoautotrophic. Furthermore, DNRA generated 14 mmol NH 4 N (Table 3 3). This result was unexpected; DNRA is typically favored when C:N ratios are very high [ Korom , 1992; Kelso et al. , 1997; Burgin and Hamilton , 2007; Rivett et al. , 2008] . Research contends that microbes that facilitate DNRA are limited to anoxic zones of the aquifer [ Hill , 1996] . This finding paired with low ammonium concentration measurements has been the reason DNRA is often assumed to be negligible in many aquifers. Relatively few studies have confirmed DNRA in aquifers [ Bengtsson and Annadotter , 1989; Thayalakumaran et al. , 2008; Hsu et al. , 2010] . Furthermore, these studies assume DNRA vis Ã vis increased ammonium concentrations alone and do not pro vide direct evidence. With NO 3 DOC injections at the anoxic site, the number of 16S rRNA gene copies increased 14 fold r elative to baseline (Figure 3 3A ), and normalized nirK copies (Figure 3 3C ) and normalized nirS copies (Figure 3 3E ) were close to base line values, suggesting that denitrification gene activity in the sampled water increased at the same rate as population. Normalized nrfA copy numbers (Figure 3 3B) and expression (Figure 3 3D ) decreased below baseline and nrfA expression rates (Figure 3 3 F ) increased. These data support the premise that nrfA
7 4 as 16S rRNA gene copies increased, indicating additional DOC and nitrate addition will not increase microbial groups with nrfA . These microbial data were associated with smaller mass of recovered NH 4 N (<0.001 mmol). Reduced NH 4 N recovery and microbial data suggest the relative importance of DNRA nitrate transformation decreased with NO 3 DOC and o ther nitrate loss mechanisms may be more important such as A NRA and denitrification in aquifer matrix biofilms. Microbial biomass change was evaluated as a mechanism for unrecovered nitrate (Table 3 4). These data show minimal PBN mass for baseline conditions and NO 3 only at the oxic and anoxic site less than 0.08 ( +/ 0.01) mmol, this is expected as significant ANRA is only expected with increased carbon availability [ Guerrero , 1985] . Increased PBN occurred with NO 3 DOC; 4.27(+/ 1.86) mmol and 2.93(+/ 1.28) mmol for the oxic and anoxic sites, respectively. Thes e values suggest that some NO 3 N was converted into biomass with fumarate addition, however pore water ANRA was much lower than the amount of observed nitrate reduction for all experiments (> 19 mmol). Thus, ANRA in the recovered injectate cannot account f or missing NO 3 N mass. Biofilm ANRA, denitrification, or DNRA may explain the unrecovered nitrate. Microbial biomass estimates in the recovered injectate can account for some of the observed nitrate mass loss but not enough to account for a significant por tion; all pore water biomass nitrogen (PBN) values were less than 4.27(+/ 1.86) mmol (Table 3 4) while unrecovered NO 3 N was greater than 19 mmol for all PPTTs (Table 3 3) . However, microbial sampling only examined pore water microbes, neglecting microbia l communities in biofilms attached to aquifer material surfaces. Biofilms may be where many N transformations are occurring. Biofilms have been extensively documented in
75 sulfidic karst caves and aquifers [ Macalady et al. , 2006; Gray and Engel , 2013] and may occur in anoxic and oxic karst aquifers and biofilms create zones of significantly different conditions relative to their surroun ding environment [ Kuenen et al. , 1986; Biesterfeld et al. , 2003] . Most microbial biomass is attached to aquifer materials [ Holm et al. , 1992] and may represent the regions of greatest metabolic activity [ Blakey and Towler , 1988] . Thus, local scale denitrification, DNRA combined with ammonium adsorption, or ANRA in biofilms may explain observed nitrate losses. Study Summary This study examined two geochemical end members in the Ichetucknee springshed which discharges groundwater from the karst UFA. Observed NO 3 N losses cannot be completely explained by estimated as similatory reduction of nitrate to ammonium (ANRA) biomass estimates, or recovered ammonium, suggesting denitrification occurred. Zero order NO 3 N loss rates support the premise that nitrate losses are occurring in the aquifer across a redox gradient with and without carbon addition. Measured rates were consistent with reported denitrification rates in other basin fill aquifer PPTTs, which range from 0.01 to 0.36 mmol L 1 hr 1 [ Istok et al. , 1997; Kim et al. , 2005] . Computed zero order rates were independent of concentration and are an indicator of the local aquifer response to nutrient addition, however it is highly unlikely that these rates are constant in the aquifer throughout space and time. If that were the case, eogenetic karst aquifers would not be susceptible to persistent nitrate contamination. In the oxic well test , zero orde r NO 3 N loss rates were low but not insignificant, 0.03 (+/ 0.01) mmol L 1 hr 1 with NO 3 only to 0.05(+/ 0.01) mmol L 1 hr 1 with NO 3 DOC. Dissimilatory reduction of nitrate to ammonium was minimal (<0.001 mmol) with NO 3 -
76 only and less than 1 mmol with NO 3 DOC. The nirS gene activity showed the most dramatic response to tracer addition, indicating an increase in denitrification genes. Measurement of nitrate loss where the aquifer is oxic suggests that more research is needed to investigate differences in N loss mechanisms in groundwater, the porous aquifer matrix and within aquifer biofilms. Similarly, nitrate losses were observed in the anoxic well test for both tracer experiments. With NO 3 only, dissimilatory reduction of nitrate to ammonium dominated N l oss; 58% of nitrate reduction was conversion to ammonium (14 mmol) and zero order NO 3 N loss rates were 0.03 (+/ 0.02) mmol L 1 hr 1 . In addition, i ncreased normalized nrfA DNA and cDNA, increased ammonium concentrations, presence of hydrogen sulfide gas, and increased iron concentrations suggest that iron or sulfide driven chemolithoautorophic dissimilatory reduction of nitrate to ammonium occurred. Significant ammonium was not recovered after NO 3 DOC addition (<0.001 mmol). However, zero order NO3 N loss rates quadrupled to 0.13 (+/ 0.02) mmol L 1 hr 1 . Negligible ammonium mass, microbial data, and significant increase in NO 3 N loss rate support that with NO 3 DOC another N transformation mechanism, such as denitrification, may be occurring in biofilms of the aquifer matrix. Results of this study indicate that nitrogen transformations in karst aquifers are more diverse than previously reported, which has implications for future study and management of nitrate enrichment. Nitrate losses occurred in both oxi c and anoxic zone with rates that varied with fumarate addition. While zero order NO 3 N loss rates were greater at the anoxic site with NO 3 DOC, nitrate transformation in oxic areas of the aquifer was not negligible, and evidence of DNRA was found. Thus, f uture efforts must
77 extend conceptual models to account for a broader suite of nitrate transformations and conditions. Observed nitrate losses without corresponding increases in nirS and nirK expression or expression rate in recovered injectate supports t he hypothesis that N transformations are likely occurring in porous media biofilms. Future work should focus on sampling of these biofilms and conducting experiments which directly investigate their influence on macroscale N transformations. Unexpected significant chemolithoautotrophic DNRA was observed, suggesting that low DOC concentrations may not limit nitrogen transformations in the UFA, alternate electron donors may facilitate both DNRA and denitrification. While significant DNRA has been documente d in a sandy Swedish aquifer microcosm study [ Bengtsson and Annadotter , 1989] and coastal marsh and sediment studies [ Burgin and Hamilton , 2007] , to the author s knowledge, this study represents the first time that chemolithoautotrophic DNRA has been documented in a karst aquifer. The fate of nitrate converted to ammonium through DNRA is not well studied and may be as important as respiratory denitrification [ Burgin and Hamilton , 2007] . However am monium will likely readily be oxidized back to nitrate in aerobic zones, and ammonium sorption onto aquifer materials is expected to be significant in many aquifer systems [ Davidson et al. , 2003; Buss et al. , 2004; Rivett et al. , 2007] Thus, DNRA is just a temporary sink for nitrate. Interactions among DNRA, sorption, and nitrification will influence the fate and transport of nitrate and complicates conceptual models for nitrate transport in the karst aquifers. Furthermore, the presence of ANRA and DNRA as N retention processes influences how microbial c ommunities take advantage of and
78 respond to nutrient addition and supports the occurrence of nutrient spiraling in the aquifer. Conceptual models of N transformation in aquifers may be improved by interpreting N transformations in terms of nutrient spirali ng.
79 Figure 3 1. Push pull tracer test site map. A) Site location in Florida illustrating Upper Floridan Aquifer confinement. B) Ichetuc knee Springshed showing streams and Ichetucknee Springs Group. C) P ush pull tracer test site locations.
80 Table 3 1 . Selected geochemical data for characterization of push pull tracer test sites. MW14 MW15 Sample date 8/2/2013 3/14/2014 8/2/2013 3/14/2015 pH 7.26 6.7 8.42 6.91 DO 7.49 8.02 0.15 0.53 ORP 89.9 154 177 141 Mg 2+ 2.67 1.99 3.16 3.04 Cl 2.78 2.7 3 .84 4.07 SO4 2 3.55 3.47 8.16 8.53 NO X N 0.11 0.03 0.07 0.12 NH 4 N 0.07 0.00 0.08 0.00 PO4 3 0.06 0.06 0.01 0.01 DOC 0.12 0.56 0.26 0.7 DIC 50.29 52.43 36.78 39.47 Fe 2+ 0.00 0.00 0.46 0.53 Mn 2+ 0.00 0.00 0.04 0.04 Table 3 2. Primers for Nitroge n Functional Genes showing functional genes hazA , nrfA , nirK , and nirS , associated nitrate transformation reaction name and process, gene sequence and length in base pairs (bp). Sequence Length (bp) hazA ANAMMOX (NO 2 + NH 4 + N 2 ) hzsA_1857R AABGGYGAAT CATARTGGC 1838 1857 nrfA DNRA (NO 2 NH 3 ) nrfA F2aw nrfA7R1 250 nirK Nitrite reductase (NO 2 NO) nirK876F ATYGGCGGVAYGGCGA nirK1040R GCCTCGATCAGRTTRTGGTT 165 nirS Ni trite reductase (NO 2 NO) nirSCd3aF AACGYSAAGGARACSGG nirSR3cd GASTTCGGRTGSGTCTTSAYGAA 400
81 Table 3 3 . Push pull tracer test summary, showing injected and recovered mass, recovered mass fraction, zero order nitrate loss rate, and zero o rder fumarate loss rate with one standard deviation errors, and resting time in the aquifer (t*). Site Test t ype Quantities Cl NO x N NH 4 N Fumarate C Oxic Nitrate Mass i njected (mmol) 259 (+/ 1) 316 (+/ 3) Mass e xtracted (mmol) 252 (+/ 2) 297 (+/ 1 ) 0 Recovered m ass 0.97 0.94 Zero order r ate ( mmol L 1 hr 1 ) 0.03 (+/ 0.01) t* (hr) 21.5 Oxic Nitrate + Fumarate Mass i njected (mmol) 240 (+/ 1) 381 (+/ 3) 314 (+/ 0) Mass e xtracted (mmol) 239 (+/ 3) 356 (+/ 2) 1 (+/ 0) 262 (+ / 8) Recovered m ass 1 0.93 0.84 Zero order r ate ( mmol L 1 hr 1 ) 0.05 (+/ 0.01) 0.81 (+/ 0.0) t* (hr) 20 Anoxic Nitrate Mass i njected (mmol) 260 (+/ 0) 306 (+/ 0) Mass e xtracted (mmol) 261 (+/ 7) 282 (+/ 1) 14 (+/ 0) Recovered m ass 1.01 0.92 Zero order r ate ( mmol L 1 hr 1 ) 0.03 (+/ 0.02) t* (hr) 10.25 Anoxic Nitrate + Fumarate Mass i njected (mmol) 262 (+/ 0) 379 (+/ 6) 311 (+/ 0) Mass e xtracted (mmol) 258 (+/ 7) 320 (+/ 1) 0 181 (+/ 1 5) Recovered m ass 0.99 0.84 0.58 Zero order r ate ( mmol L 1 hr 1 ) 0.13 (+/ 0.02) 0.78 (+/ 0.0) t* (hr) 20.5
82 Figure 3 2. Push pull tracer test cum ulative mass recovery results. Chloride (Cl), nitrate + nitrite nitrogen (NO x N), ammonium nitrogen (NH 4 N) and fumarate recovery shown with shaded gray envelopes that represent standard deviation in cumulative mass recovery uncertainty. A) O xic aquifer response to nitrate. B) N itrate and fumarate (DOC) additions. C) A noxic aquifer response to nitrate . D) N itrate and fumarate additions.
83 Figure 3 3. Microbial genetic data showing baseline conditions and response to tracer additions. Error bars show standard deviation in measured data. A) M icrob ial DNA copies (16S rRNA genes). B) N ormal ized nrfA gene copies. C) N ormalized nirK copies. D) normalized nrfA gene expression (cDNA) . E) N ormalized nirS gene copies. F) N ormalized nrfA cDNA gene expression rate. Table 3 4. Estimated pore water biomass nitrogen from assimilatory reduction o f nitrate to ammonium in recovered injectate for all push pull tracer tests. PPTT Biomass N (mmol) Oxic Baseline 0.08 (+/ 0.01) Oxic Nitrate + Fumarate 4.27 (+/ 1.86) Anoxic Baseline 0.02 (+/ 0.01) Anoxic Nitrate 0.01 (+/ 0) Anoxic Nitrate + Fumar ate 2.93 (+/ 1.28)
84 CHAPTER 4 EXAMINING DRIVERS OF KARST AQUIFER HYDROLOGIC RESPONSE USING STOCHASTICALLY GENERATED PREFERENTIAL FLOW PATH TEMP LATES AND GLOBAL SENSITIVITY ANALYSIS Background Karst Aquifers Karst hydrologic flow paths include d rainage of surface runoff through sinkholes, temporary storage of water within a near surface epikarst zone, storage and laminar flow through the porous limestone matrix, turbulent flow though conduits, and discharge from conduit networks to perennial spri ngs [ Taylor and Greene , 2005] . Hydrologic flowpaths are controlled by the structure and organization of sinks and conduits, which collect surface drainage and infiltrating water and provide subsurface flowpaths with lower resistance and hydraulic gradients than the surrounding porous media [ White , 1993] . This aquifer behavior results in convergent flow to large springs, troughs in potentiometric surfaces, decreased hydraulic gradients , and increased hydraulic conductivity in the downgradient direction [ Worthington , 2009] . Karst dissolution processes lead to preferential subsurface flow, conduits, anisotropic groundwater flow, sinkholes and sinking streams or swallets [ Ewers , 2006] . Hydraulics in karst aquifers are controlled by high permeability features (e.g. fractures, conduits, and sinks) embedded in a lower permeability porous limestone matrix [ Atkinson , 1977] . While conduit networks typically represent 0.1 1% of the aquifer porosity, they often convey much of the flow in karst aquifers [ Bonacci , 1987; Worthington , 1999; Ford and Williams , 2013] . Karst flow heterogeneity is a function of conduit induced hydraulic conductivity contrasts and conduit densit y [ KovÃ¡cs and Perrochet , 2008] .
85 Karst Aquifer Development Research suggests that UFA karst has developed over multi episodic dissolution events that reinforce preferential pathways when karst is exposed during sea level retreat [ Hanshaw and Back , 1979; Hunn and Slack , 1983; Randazzo and Bloom , 1985; Randazzo , 1997; Florea et al. , 2007; Gu lley et al. , 2013a, 2013b] . Karst aquifer formation (e.g. speleogenesis and conduit network evolution) is governed by boundary conditions (such as recharge rates, topography and sea level), porous matrix properties, horizontal preferential flow pathways ( HPFs; e.g. fractures, protoconduits, bedding planes) and vertical preferential flowpaths (VPFs; e.g. root networks, fractures, sinkholes) geometry, topology, and location [ Palmer , 1991; Scott et al. , 2004; Taylor and Greene , 2005; Ford and Williams , 2013] . HPF and VPF properties define a preferential flow path temp late and exert primary influence on karst development [ Palmer , 1991] . These pathways deliver recharge which is undersaturated with respect to calcite that contributes to karst development by driving dissolution of lime stone [ Pa lmer , 2003] . During karst development, dissolution along these preferential flow paths leads to self organized conduit networks [ Worthington and Ford , 2009] . Dissolution leads to a positive feedback where larger con duits capture more flow and are preferentially enlarged. As preferential paths are dissolved undersaturated water penetrates further along the conduit flow path ultimately resulting in cal cite undersaturated water exits the conduit into the spring . Often, at this point flow along the length of an entire conduit is turbulent [ Dreybrodt , 1990; Palmer , 1991] . Undersaturated water throug hout the conduit allows for rapid enlargement of conduits within the network, leading to decreases in hydraulic heads,
86 flow field reorientation, development of tributary conduits, and increased springflows [ Ewers , 2006; Ford and Williams , 2013] . Characterization of hydrologic flow paths in karst aquifers relies on well described conduit geometry, however few conduits can be directly observed [ Palmer , 2006; Pardo IgÃºzquiza et al. , 2012] . Traditional methods of aquifer characterization such as aquifer tests, borehole log s, and geophysical observations provide little information about the spatial configuration and properties of conduit networks [ Teutsch and Sauter , 1991] . Indirect geochemistry methods [ Birk et al. , 2004; Raeisi et al. , 2007] and spring hydrograph analysis [ Bonacci , 1993; Baedke and Krothe , 2001; Geyer et al. , 2008a; Covington et al. , 2009; Bailly Comte et al. , 2011; Fiorillo , 2014; Katsanou et al. , 2015] have been widely used to obtain information about macro scale properties of karst aquifers. Karst Aquifer Modeling Much work has been done desc ribing karst formation processes [ KirÃ¡ly and Morel , 1976; Kaufmann and Braun , 2000; Liedl , 2003; Kaufmann , 2009] and modelling flow and transport in karst systems [ Kiraly , 1985; Dreybrodt , 1990; Palmer , 1991; Eisenlohr et al. , 1997a, 1997b; de Rooij et al. , 2013] . Recent advances in integrated karst conduit and matrix flow models, discrete continuum models, have made simulation of karst networks in regional scale models possible (e.g. MODFLOW CFP, [ Shoemaker et al. , 2005] ; DISCO, [ de Rooij et al. , 2013] , FEFLOW [ Diersch , 2005] ). However, accurate karst aquifer groundwater flow and transport models requ ire well defined hydraulic and geometric parameters of conduits and the porous matrix [ KovÃ¡cs et al. , 2005] at fine spatial scales. In the few applications of basin to regional scale karst hydrologic models to date [ Shoemaker et al. , 2008; Gallegos et al. , 2013; Saller et al. ,
87 2013; Xu et al. , 2015] conduit properties (e.g. conduit geometry, tortuosity, radius, roughness, and conduit matrix exchange) were typically specified a priori or , in some cases , calibrated. Karst Conduit Evolution Algorithms There have been significant efforts to simulate the evolutio n of realistic karst networks that honor boundary conditions and hydrogeology over geologic timescales. These network generating processes can be separated into two approaches, structure imitating and process imitating. Structure imitating approaches use stochastic concepts to reproduce conduit networks. These methods include stochastically generated fracture networks [ Andersson et al. , 1984; Cacas et al. , 1990; Reeves et al. , 2008] , geostatistical methods [ Jaqu et and Jeannin , 1994; Fournillon et al. , 2010] fractals [ Maramathas and Boudouvis , 2006] , and lattice gas cellular automa [ Jaquet et al. , 2004] , However these methods typically do not consider the underlying processes driving karst conduit developme nt and may lead to poor network connectivity. Pseudogenetic algorithms [ Borghi et al. , 2010, 2012, 2016] , and inverse modelling [ KovÃ¡cs , 2003] address this issue by incorporating some of the physical concepts driving conduit network formation. Process imitating methods use dissolution models to reproduce conduit evolution and speleogenesis [ Lauritzen et al. , 1992; Groves and Howard , 1994; Howard and Groves , 1995; Dreybrodt , 1996; Kaufmann , 2009; Worthington and Ford , 2009; Kaufmann et al. , 2010] . These methods are limited b y lack of information about preferential flow path temp lates, and initial and boundary conditions. Process imitating models have made significant progress in describing karst dissolution in small scale box models or in networks of regularly ordered and co nnected
88 fractures [ Dreybrodt , 1990; Kaufmann , 2003, 2009; Reimann et al. , 2011a, 2011b] . However, no work has us ed process based speleogenesis to examine how a conduit network evolves at the basin scale under assumptions of random numbers and locations of VPFs, random HPF fracture orientation and length scales, and random porous media properties Evaluati on of flow and transport response in these types of evolved networks is need ed to provide insight into what preferential flow path template attributes drive springshed behavior. Sensitivity Analyses Sensitivity analysis (SA) provides a mechanism to identif y and prioritize influential model inputs, identify non influential parameters so that they can be fixed to representative values, and calibrate model inputs to achieve desired model behaviors [ Saltelli et al. , 2004; Pappenberger et al. , 2008] . SA methods are highly useful for understanding models of complex non linear systems. Much research has documented the benefits an d limitations of SA methods [ Foglia et al. , 2009; Nossent and Bauwens , 2012; Herman et al. , 2013; Li et al. , 2013; Shin et al. , 2013] . Sensitivity analysis techniques have been widely applied to hydrological models [ van Griensven et al. , 2006; MuÃ±oz Carpena et al. , 2007; Pappenberger et al. , 2008 ; Srivastava et al. , 2014] and contaminant transport models [ James and Oldenburg , 1997; Pan et al. , 2011] . The complexity of SA analyses ranges from one at a time sensitivity analysis (OAT SA) to global sensitivity analysis (GSA). In OAT SA only one variable is changed and a response metric is evaluated relative to that change. The main difference between OAT SA and GSA methods is that GSA methods account for non linear interactions among parameters. This increased information from GSA methods comes
89 at a computational cost, often requiring thousands of model executions , l eading to limited application of variance based GSA in hydrologic models [ Shin et al. , 2013] . An alternative to OAT SA and full variance based GSA are a subset of GSA methods called global screening methods (e.g. the Morris method, [ Morris , 1991; Campolongo et al. , 2007] ). The Morris method provides a more rigorous evaluation of model parameter importance than single variable OAT SA methods. Global screening methods provide qualitative information about the effective parameters and interactions with significantly fewer model executions. While variance decomposition GSA methods like the method of Sobol [ Sobol , 2001] and the extended Fourier Amplitud e Sensitivity Test (FAST) [ Saltelli et al. , 1999] provide detailed quantitative information about parameter contributions to uncertainty, recent studies have illustrated that the Morris method accurately identified the non influential and most sensitive mode l parameters with significantly fewer model executions [ Confalonieri et al. , 2010; Herman et al. , 2013; Li et al. , 2013; Srivastava et al. , 2014] . To date, SA methods employed in many karst hydrologic and dissolution studies have focused on OAT SA [ KirÃ¡ly and Morel , 1976; Kaufmann and Braun , 2000; Birk , 2003; Scanlon et al. , 2003; KovÃ¡cs et al. , 2005; Perrin et al. , 2007; Doummar et al. , 2012] and local sensitivity measures [ Dafny et al. , 2010; Mazzilli et al. , 2012] . Karst modelling using OAT SA has been used extensively to evaluate the effect of properties of the porous matrix (e.g. hydraulic conductivity) and properties of well defined conduit geometries (e.g. roughness, exchange parameters) on spring hydrograph baseflow recession [ KirÃ¡ly and Morel , 1976; Eisenlohr , 1996; Eisenlohr et al. , 1997a; Cornaton , 1999] . Recently, OAT SA of conduit conductivity and frequency was used to define two
90 end members of karst domains : a conduit influenced flow regime (CIFR) where recession is con trolled by the conductive capacity of the conduit system, and matrix restricted flow domain (MRFD) where recession is controlled by the low permeability porous matrix [ KovÃ¡cs et al. , 2005] . Study Objectives Coupled conduit network evolution, flow, and transport has ye t to be evaluated in a rigorous global sensitivity analysis (GSA) framework. In this study, the Morris GSA method wa s used to examine how stochastically generated preferential flow path temp lates evolve into karst conduit networks that create first magnitu de springs (flow > 2.8 m 3 s 1 ), gain insight into conduit network evolution processes, and assess hydrologic and transport pulse response in evolved networks. A discrete continuum flow and transport model that simulates carbonate dissolution was incorporat ed into an iterative karst development framework [ Borghi et al. , 2010] . Conduit networks and resulting s pring systems were evolved through dissolution of stochastically generated preferential flow path temp lates. In this framework , sink to spring connectivity is not specified a priori an d resulting network development and spring generation is not guaranteed. This study addressed three questions. First , what preferential flow path temp late properties lead to development of a first magnitude spring system ? To answer this question, an ensemb le of preferential flow path temp lates were evolved until a steady configuration was achieved. Evaluation of ensemble behaviors during dissolution, such as conduit development along HPFs, plots of network radii, Reynolds numbers which describe the distribu tion and development of turbulent flow through the domain, and changes in hydraulic heads and springflows were used to address this research question.
91 Karst networks exist across a gradient of poorly developed to mature and develop over multiple events ov er geologic timescales. Thus a second research question ad dressed in this stu dy was how much time t o dissolve networks to allow direct comparisons of hydraulic and transport pulse response among networks ? Similarities in network evolution among replicates in Monte Carlo ensembles were used to evaluate when a replicate had evolved into a mature karst drainage network. Conduit properties are uncertain and water managers need to understand how this uncertainty propagates into uncertainty in predicting springsh ed hydrologic and transport pulse response . A third question this research addressed wa s how do initial preferential flow p ath assumptions influence flow and transport response in evolved networks? To evaluate this question, ensemble statistics were used t o evaluate the mean and variance of springshed hydrologic and transport pulse response. Furthermore, GSA me thods examined which preferential flow path parameters mo st strongly influence hydrologic and transport metrics. It was hypothesized that (1) a limi ted range of parameters exi s ts which result in a development of a first magnitude spring system that is the main outlet for recharge fluxes in a hypothetical springshed. Parameter interactions that increase connectivity (e.g. fracture length, fracture dens ity, VPF density) and decrease matrix permeability (e.g. hydraulic conductivity) will be the most influential; (2) spring hydrologic and transport pulse response uncertainty in systems where a first magnitude spring develops will be reasonably similar acro ss different conduit networks, but distributed aquifer behavior (head fields, regions vulnerable to surface solute applications) will be more highly variable and depend ent on conduit configuration ; ( 3) first magnitude spring
92 conduit networks subject to the same boundary conditions will have conduit development in similar regions of the springshed, allowing for the mapping of vulnerable areas based on climatic, topographic and hydrologic features. The conceptual geologic model and hydrologic model for this s tudy were based upon the Silver Springshed in North Florida (Figure 4 2) . Methods Overview of Flow Transport and Conduit Dissolution Model The conduit generation algorithm employed in this study adapts a pseudogenetic scheme [ Borghi et al. , 2012] to include a process based flow, transport , and conduit dissolution model, stocha stically generated preferential flow path temp lates and boundary conditions that honor what is known about the geology, hydrology , and climate of the system being modeled (Figure 4 1). The physics based discrete continuum flow and transport model DisCo [ de Rooij et al. , 2013] was used to ev olve karst networks from initial preferential flow path temp lates rather than the more empirical using fast marching algorithm used by Borghi et al.  . This approach is advantageous because it incorporates the physics of speleogenesis into dev eloping conduit networks from the initial preferential flow path temp lates and uses the advection dispersion equation to represent calcium transport. Details underlying the DisC o model are included in Appendix B. The process start s with a geologic model of the Silver S prings springs hed to specify topography, geologic model structure, springshed boundaries and boundary conditions, and the spring outlet (Figure 4 2) . Then stochastic preferential flow path temp lates are generated with random vertical preferent ial pathways (VPFs), random horizontal preferential pathways (HPFs), and a fixed outlet for Silver Springs. The
93 conduit evolution algorithm is run iteratively to evolve conduit netwo rks using coupled flow, calcite dissolution, and transport equations in th e DisC o model . At each successive iteration , conduits are enlarged based on the mass of carbonate dissolved, limestone density, and fracture surface area. An advantage to this algorithm is that conduit properties evolve from initial conditions using physic s describing dissolution and transport of calcite undersaturated water. Furthermore, exchange between the conduit network and porous matrix is a function of model discretization and porous media properties via the Peaceman well index instead of a calibrate d parameter [ Peaceman , 1978, 1983 ] . Hydrologic Model The conceptual DisCo hydrologic model was derived from the Silver Springs local scale model (Figure 4 2), a subdomain of the Northern District Model (NDM) [ Hydrogeologic , 2013] . All surface topography, layering, a nd total Floridan Aquifer depths were taken from the NDM but porous media properties were homogenized and boundary conditions were simplified to facilitate analysis of heterogeneity induced by variations in preferential flow path temp lates. The model conta ined 7 layers, 161 rows and 92 columns with 762m s quare model cells (Figure 4 3). The active domain wa s defined by the springshed boundary determined from the NDM 1000 year capture zone [ Hydrogeologic , 2013] . Model layer 1 represents the surficial aquifer and epikarst, which collects and conveys rec harge to the karst aquifer through both focused (VPF) and diffuse recharge. Layers 2 through 7 represent the limestone aquifer . Despite, homogenous property assumptions, NDM model layer discretization was maintained to facilitate incorporation of conduit n etwork replicates in future local scale hydrologic decision models.
94 Silver Springs is primarily sourced from Ocala Limestone [ Faulkner , 1970 ] which is represented by layer 4 in the original NDM. Most reported cave systems in the UFA form at the contacts between sha llow stratigraphic units (e.g. Ocala Limestone and Avon Park formation contact) [ Beck , 1986] . The interface between the Ocala Limestone and Avon Park Formation is presumed to be the location of extensive karst development because of the lower relative permeability of the Avon Park Formation [ Phelps and Survey , 2004; Flore a and Vacher , 2007b] . Prior studies in the UFA in north Florida [ Florea et al. , 2007; Langston et al. , 2012a; Gulley et al. , 2013a, 2013b] suggest that conduits are located 1 0 15m below sea level, the top and bottom median elevation values for model layer 4 are 4.23 and 37.61, respectively. Accordingly, HPFs were placed throughout layer 4 to generate preferential flow path temp lates. Boundary conditions Lateral no flux bo undaries were placed at along the defined regulatory springshed boundaries, resulting in 4 179 active cells in each layer. A fixed head boundary represented the spring with the Silver Spring pool elevation, 12.198m above sea level [ Hydrog eologic , 2013] . A specified flux boundary represented steady state surface recharge of 1.182E 8m s 1 which is based upon the mean areal recharge for the NDM model within the Silver Spring subdomain. The resulting total steady state recharge flux is 28.65m 3 s 1 which is greater than average springflow for the Silver Spring of 20.9m 3 s 1 calculated from 1935 to 2014 [ USGS , 2014] . B oundary conditions during the conduit evolution process are generally unknown. However, applied recharge during conduit evolut ion must be carefully considered to avoid unrealistic scenarios such as unrealistic high hyd raulic heads above land surface . Simulating surface flow routing and storage using rigorous mass balance methods is
95 computationally demanding and not feasible over the time frames requ ired for conduit development. Therefore a computationally efficient surface boundary condition was developed to remove surface water from the model domain so that unrealistic hydraulic heads do n ot occur before the network has been able to develop and capture recharge . For each cell, a spill elevation was defined using the surface topography [ Wang and Liu , 2006] . Spill elevations were computed based upon adjacent cell elevations . If the cell was located in a local depression, the spill elevation was set to the height of the lowest cell elevation along the closed topo graphic depression boundary. This boundary limited the hydraulic heads to the spill elevation. The surface flux was calculated as: (4 1) where is the surface flux, is the drain conductance, is the hydraulic head in the topmost cell, and is the spill elevation associated with the drain. As the conduit network evolve d through time heads throughout the domain we re lowered and flow wa s captur ed from the surface bounda ries and routed to the spring. Preferential flow path temp late generation Karstification requires initial apertures of some nominal size to facilitate dissolution . These are referred to here as the preferential flow path temp late and characterized by the random length, density and orientation of HPFs and random location of VPFs . An ensemble of possible preferential flow path temp lates were created and implemented in the hydrologic model to obtain insights into hydrologic response uncertainty. The conceptua l geologic model for local HPFs in this study is based upon knowledge of HPF development along bedding planes of the Ocala Limestone and northwest and northeast trending fractures and joints developed since the Oligocene
96 epoch approximately 35 m.y.a [ Vernon , 1951; Ceryak et al. , 1983] (Figure 4 2) . The network of HPFs were conceptualized as intersecting fracture segments within a horizontal plane that were mapped to layer 4 of the hydrologic model, which corresponds to the Ocala Limestone . A limitation of this approach is that karst net works only evolve within the initial HPFs. Each HPF was generated as a line segment using typical stochastic fracture generation techniques [ Bauer et al. , 2005] . In order to facilitate connectiv ity , several geometric restrictions were applied to the HPF generation algorithm; (1) to prevent superposition of HPFs on one another a minimum distance between HPFs of 10 m was established, (2) to prevent very small HPF segments a minimum distance of 10 m was specified between each HPF intersection point, and (3) to prevent slightly subparallel HPFs a minimum angle of intersection of 5Â° between HPFs in the same fracture set was implemented. Length distribution of HPFs was determined based upon a gamma dist ribution with two parameters, k and theta. For the gamma distribution; k values less than 1 approach power law behavior, k values equal to 1 approach exponential behavior and k values greater than one represent distributions with long tail behaviors. The t heta parameter is a scale value. The expected mean HPF length for each dis tribution is equal to k*theta. All HPFs were set to an initial diameter of 2 mm. Sinks and other VPFs have been observed near photolinear features, are common in valley floors [ Littlefield et al. , 1984] , and are often associated with the oldest karst depressions [ Upchurch and Littlefield Jr , 1988] represented by high hydraulic conductivity cells connecting the land surface of the model to the conduit layer (laye r 4). In a previous UFA study near the Suwannee River,
97 mean sinkhole densities were approximately 6 per km 2 [ Denizman , 1998] , and have been reported to range between 3.5 to 19.9 per km 2 in other temperate subtropical settings [ Ford and Williams , 2013] . In the Silver Springs springshed, density of sinkholes may be estimated based upon Florida Geological Survey data (Figure 4 2) whic h was based upon geographic analysis of surface topography, however it is unknown how many of the mapped sinks or depressions extend to the aquifer. Due to the uncertainty in the locations and spatial density of HPFs and VPFs in the Silver Springshed, vari ous densities and configurations were evaluated in the preferential flow path temp late algorithm. Each preferential flow path temp late contained a specified number of stochastically generated gamma distributed HPFs, a specified number of randomly located V PFs, and homogeneous matrix porosity, hydraulic conductivity, specific storage coef ficient. Values for the parameters of the gamma distribution, HPF and VPF densities, and porous media properties varied across the realizations using values drawn from a spe cified probability distribution (Table 4 1). Across the ensemble, HPF fracture sets had mean orientations of 45Â° and 315Â°, northeast and northwest, following orientations of a prior photolineament study [ Vernon , 1 951] , however the orientation of individual HPFs were allowed to randomly vary around this mean with a tolerance drawn from a uniform distribution . This random tolerance is referred to as the spread in following analyses. The l ocation of each HPF and VPF was randomly selected from a uniform distribution of x and y coordinates within the domain . One HPF was deterministically placed at the location of Silver Springs.
98 Dissolution model DisCo was run for each preferential flow path temp late to simulate the co nduit network evolution subject to the specified climatic, topographic, and hydrologic boundary conditions of the Silver Springshed model. Calcite undersaturated water was applied at land surface, and moved through the epikarst and VPFs to the HPFs where d issolution could occur, enlarging the HPF radii. Recharge entering the HPFs from the VPFs was assumed to 100% undersaturated with respect to calcite. Groundwater entering the HPFs from the limestone porous media was assume d to be 100% saturated with respec t to calcite. Morris Method Global Sensitivity Analysis Morris Method Global Sensitivity Analysis (MM GSA) was used to sample the wide range of possible preferential flow path and porous matrix properties influencing spring network genesis and hydrologic a nd transport response. Furthermore, insights from the MM GSA were used to evaluate properties driving network development, observed behaviors, and examine sensitivity to input parameters. The MM GSA used estimated ranges in values for HPF density and orie ntation, VPF location and density, and porous matrix properties from literature and previous studies in the Silver Springshed (Ta ble 4 1). The analysis evaluated dissolution and hydraulic and transport pulse response. Not all preferential flow path temp lat es have the potential to create a first magnitude spring, thus behavioral Monte Carlo filtering (BMCF) of the MM GSA results was used to determine parameter ranges most likely t o create hat is, preferential flow path temp lates that develop into a first magnitude spring.
99 The Morris method classifies input parameters into three groups, parameters that have negligible effects, those with large linear effects without interactions, and those with large non linear and interactive effects [ I ooss and LemaÃ®tre , 2014] . The Morris a k dimensional p level grid [ Campolongo et al. , 2007] ; where k is the number of independent model input parameters (i.e. input parame ter space is a vector with k ) and p is the number of levels which span the specified value range for each parameter. For each model run , a change in a single parameter is generated and it s effect on the output ( x ) is examined. The change in x in response to a change in the i th parameter ( ) is the elementary effect ( ) of the i th input on output x calculated as: (4 2) where represents a non dimensional step change of parameter within the specified range, i.e . i /( p 1) where i p 2. For this study, p =4 was selected [ Khare and MuÃ±oz carpena , 2014] . Differen t sample of . These effects are used to derive sensitivity indices based upon the mean ( Âµ ) and the standard deviation ( ) of the over the ensemble. Large Âµ values indicate a large influe nce of the corresponding parameter on model output. A modified Âµ value, Âµ * was proposed by [ Campolongo et al. , 2007] , which averages the absolute values of to eliminate effects of opposite signs in non monotonic model responses. This Âµ * index was shown in the same study to be an acceptable approximation for the variance based total index provided by more rigorous GSA analyses. The signs of elementary effects were shown to vary during analysis so both Âµ and Âµ * were evaluated.
100 Large v alues of indic ate that the corresponding input parameter exerts a non linear effect on model output or has non linear interactions among other input parameters. The number of simulations in the Morris analysis ( ) is: (4 3) w here is the sampling size for each trajectory and is the number of parameters being investigated. Parameters were sampled using the eSU method [ Khare and MuÃ±oz Carpena , 2014] . A value of =10 was selected based on previous research [ Saltelli et al. , 2004; MuÃ±oz Carpena et al. , 2007] and = 8 for the MM GSA analysis (Table 4 1) resulting in = 90 simulations. To better control random behaviors, the random seed for each replicate was changed for each gamma dis tribution. Morris results were analyzed using the Elementary Measures and Plots Matlab software [ Khare an d MuÃ±oz Carpena , 2014] . The blue lines in the plot reflect a 1:1 line 2 . Values inside the red plotted lines are highly interactive, those within the blue lines are moderately interactive and those outside the lines are non interactive parameters. Dissolution Response Metrics The spring hydrograph as a function of dissolution time was used to determine when the conduit network reached a stable configuration. Hydraulic and transport pulse exper iments were evaluated on the stable network. Differences in conduit radii statistics (median , m ax ), and network plots were used to compare conduit networks. Effective porosity for the evolved karst aquifer was computed as: (4 4)
101 where is the conduit volume, is the matrix volume, and is the porosity. The HPF density was defined as: (4 5) where is the total length of HPFs and is the spri ngshed area and has units of [km 1 ]. Dissolved HPF density was used to evaluate magnitude of dissolution across replicates. For that metric, is the total HPF length for HPFs with radii greater than 0.5 m was used in equation 4 5. Hydrograph Res ponse Metrics Hydrograph response metrics were derived directly from simulated spring hydrograph storm pulse responses. Metrics include steady state flow magnitude, timing and magnitude of peak flow, and early and late hydrograph recession coefficients. Hy drograph recession can be separated into at least two periods, an early period which is controlled by drainage of water in the conduit network, epikarst, and matrix cells containing conduits, and a later period that is controlled by slow draining of storag e in the lower permeability porous matrix often referred to as baseflow recession [ Atkinson , 1977; Padilla et al. , 1994; Kiraly , 2002; KovÃ¡cs et a l. , 2005; Geyer et al. , 2008b] . Due to significant difference in response timescales for focused and diffuse recharge, early hydrograph behavior is dominated by focused recharge while late recession is dominated by diffuse recharge [ Covington et al. , 2009] . Recession coefficients and inflection points were calculated based upon linear regression of the natural logarithm springflow through time (lnQ vs t plots). Defining inflection points is subjective, however the point where the hydrograph recession becomes linear on the
102 lnQ vs t plot (i.e. exponential deca y) has been widely used to describe late recession behavior [ Baedke and Krothe , 2001; Dewandel et al. , 2003; KovÃ¡cs and Perrochet , 2008; Chang et al. , 2015] . This was the approach used to separate early and late reces sion behavior (Figure 4 4). The effects of conduit networks on flow toward the spring were evaluated in using fraction of conduit flow through 3 concentric circular control planes through the aquifer located 2, 8, and 12km away from a spring for steady sta te and peak flow. Transport Response Metrics Transport response metrics were derived from surface tracer pulse chemographs ( b reakthrough curves ) . A uniform instantaneous pulse of 1 kg m 2, was applied to the surface in a steady flow field. Surface and spr ing chemograph analysis included peak magnitude and timing, and moment analysis including total mass (0 th moment), mean travel time (normalized 1 st moment), and travel time standard deviation (square root of normalized 2 nd central moment). Since the experi ments allowed two boundaries for tracer mass to exit, the fraction of mass delivered to the spring or to the surface boundary varied among replicates . Combined t racer mass recovery from the surface and spring boundary conditions was required to be at least 90% before moments were calculated . The concentrations at each boundary at the end of the simulation were examined to get a relative sense of whether the remaining solute was contained in near surface storage or aquifer storage . Behavioral Monte Carlo Fil tering Behavioral Monte Carlo filtering (BMCF) was used to evaluate parameter ranges that were more likely to generate a spring conduit network. If the preferential flow path temp late evolved into a first magnitude springshed (>2.8 m 3 s 1 ), it was consider ed
103 behavioral. Sensitivity indices of behavioral classification s to MM GSA parameters w ere evaluated to determine which parameters were driving the behavior . A probability distribution was developed for each parameter for the behavioral and non behavioral groups. Kolmogorov Smirnov (KS) two sample test was used to test the null hypothesis that the distributions were identical . Results Dissolution Beha viors a nd Sensitivity The MM GSA provided insights into paleo karst template properties that lead to deve lopment of a spring conduit network. Steady state springflow was used to determine if networks generated 1 st magnitude springs. As networks evolved, they began to capture more of the surface recharge, consequently spring flux increased through time. Networ ks reached a steady spring flux over varying timescales (Figure 4 5A), however every conduit system reached a maxima after which steady state flux did not change significantly (<0.01% change in springflow). Of the 90 simulations in the MM GSA only 61 were behavioral. To test that non behavioral networks continue to be non behavioral with longer dissolution, dissolution was allowed to proceed, and springflow evolve, for 1.5 times the time when dissolution steady flux behavior was observed for each replicate . Results confirmed that non behavioral networks continue to be non behavioral (Figure 4 5B). Although a steady spring flux was achieved at a certain point during dissolution, the HPFs continue d to dissolve and radii increase d . This ha d minimal effect on s teady state spring fluxes but a significant effect on peak response to the hydraulic pulse (Table 4 2). As the conduit diameter in the network increase d , the head gradient across the conduit d ecreased and
104 the conduit be haved like a fixed head boundary equa l to the spring head, with reduced resistance to flow. Prior research of single conduit evolution in box models with fixed boundary conditions at each end has shown similar dissolution progression, where initial conduit development is sl ow but reaches a s teady value [ Kaufmann , 2003, 2009; Kaufmann et al. , 2010] . Figure 4 5 illustrates this behavior a t the macro scale for a conduit network. It is posited that steady flow behavior represents karst system maturity, and the majority of flow paths that will connect to the network have become stable. To test this hypo thesis, plots of conduit flow at steady state (15,000 years), 20,000 years, and 25,000 years were examined for one behavioral case during network development (Figure 4 6). While some very small paths continue to develop near the distal edges of the fractur e network over time , the majority of the network stays the same. These data support that when the steady flux is obtained that the configuration of the network has become stable ; this state could be described as self organized. Table 4 3 summarizes differ ence between behavioral and non behavioral groups for dissolution. Minimum, maximum, mean and standard deviation in t emplate HPF density wa s similar between the groups, while dissolved HPF distributions have different behaviors. This congruence between beh avioral and non behavioral initial network HPF densities indicates that the statistical parameters dictating the probability distribution preferential flow path temp late parameters may not be the only important driver of behavioral networks; random placeme nt of HPF and network connectivity within the domain may assert significant control over patterns in connectivity and eventual network development.
105 As expected, radii we re more developed , with higher ensemble median and mean radii in the behavioral replic ates. The small influence of conduit networks on effective porosity wa s consi stent with prior karst modelling studies, conduits make up a maximum of 1 % of the aquifer volume [ Bonacci , 1987; Worthington , 1999; Ford and Williams , 2013] . The MM GSA dissolution results (Figure 4 7) showed that hydraulic conductivity was the most sensitive parameter and it had relatively low interactions. Higher values of hydraulic conductivity contributed to lower dissolved HPF density (i.e. density of HPFs that have a radii greater than 0.5 m) (Figure 4 7A ), and lower media n and maximum radii (Figure 4 7B, Figure 4 7C ). I n this experimen t , distributed recharge was applied across the top surface and flow focusing though the VPFs to the conduits occurred if conduit heads decreased and as a result gradients between the epikarst and conduits increased. H igher hydraulic conductivity values all ow ed distributed recharge to be distributed easily within the porous matrix, reducing flow focusing in VPFs and producing lower head gradients between the matrix and conduits and less impetus for c onduits to develop. Thus in this experiment higher hydrauli c conductivity was observed to impede network development. Network development was also sensitive to parameters that specified connectivity of the preferential flow path s . Dissolved HPF density was sensitive to the HPF l ength distribution parameters (k and theta ), and fracture spread. Median radii were sensitive to length scale (theta) , which influences network connectivity . Maximum radii was sensitive to the number of VP F s . Thus higher connectivity in the preferential
106 flow path temp late due longer HPFs ( hi gher K*theta), higher variation in HPF orientation and more VPFs led to more network development. Steady State a nd Hydraulic Pulse Response Sensitivity Steady state flow (Figure 4 8A) was most sensitive to hydraulic conductivity, and the fracture length s cale (theta ). Increased hydraulic conductivity was associated with lower steady state fluxes, lower peak flows, and slower recession , which is somewhat counter intuitive. However because hydraulic conductivity had such a strong influence on conduit develop ment, the observed sensitivity wa s dominated by the indirect influence on conduit development rather than the direct influence of porous media flow behavior. This explanation is supported by the fact that high est values of matrix hydraulic conductivity (1e 4 m s 1 ) often result in non behavioral systems. Of the 20 simulations with the high est hydraulic conductivity values, only 4 generated a first magnitude spring. Non behavioral networks had steady state fluxes <0.1 m 3 s 1 and no hydraulic pulse response o r recession. After hydraulic conductivity, hydraulic response were most sensitive to the fracture length scale parameter theta and the number of VPFs. Longer fractures and more VPFs lead to increased peak flow and faster recession wh ich wa s expected becau se more developed conduits will lead to rapid response and rapid recession of hydraulic pulses. Transport Pulse Response Sensitivity Behavioral Morris transport pulse response results also showed that hydraulic conductivity was the most sensitive least int eractive parameter (Figure 4 9A). Higher hydraulic conductivity decrease d the magnitude of peak solute concentration and increase d the peak travel time, first normalized temporal moment (mean travel time)
107 and the second centralized moment ( travel time stan dard deviation ). Again this result wa s likely attributed to the fact that high hydraulic conductivities promote more flow through the matrix and inhibit development of the conduit network , thus more solute wa s transported more slowly through the porous med ia matrix. Transport pulse response was also sensitive to the fracture length scale parameter and spread . In general , increased network connectivity from longer fractures and more variable orientations increased the magnitude of peak solute concentration , and reduce d the travel time peak, mean , and standard deviation . Thus as expected, parameter values that increased the conduit network connectivity led to faster transport by focusing flow and solute into the network (Figure 4 9D). Hydraulic Pulse Respons e Metrics Ensemble hydrograph plots illustrate the variability in spring flow over the entire ensemble (Figure 4 10) and behavioral replicates (Figure 4 11). Behavioral replicate spring flows ranged from approximately 83 % to 100% of recharge. The mean stea dy state spring flow was 27 . 19 m 3 s 1 and the standard deviation around this mean was 3.70 m 3 s 1 . Hydrograph metrics for behavioral and non behavioral replicates (Table 4 4) illustrate expected behavio rs; steady state and peak flow wa s substantially higher in behavioral cases. Hydrograph recession values are in the range of values estimated from the Silver Springs data , where the best fit of transfer function hydrologic response modeling indicated exponential function with a response time ranging between 1 .5 4.5 years [Jim Jawitz, personal communication, September 1, 2015] . In this study, early recession hydrologic response time (1/alpha) ranged from 0.5 5 to 4.5 years with a mean of 0. 25 years. Late baseflow recession response time (1/alpha) ranged from 3 t o 274 years with
108 a mean of 20 years . The lower numbers are closer to the estimated values from the exponential response function lag times , but the larger values suggest that storage release during baseflow recession can extend for long periods of time. F low through concentric control planes toward the spring provides further insight into causes of behavioral and non behavioral groups (Table 4 5). In behavioral cases, flow through the control planes w as dominated by conduit flow with nearly all of the flow occurring in the conduits and small decreases in percentage of conduit flow with distance from the spring. In contrast, non behavioral case mean behavior shows small flow fractions near the spring and increasing mean flow fractions as distance from the sp ring increases. This suggests that in these non behavioral cases local flow systems developed and routed water to the surface boundary instead of the spring ; t his finding wa s further supported by observed non behavioral networks discussed below . Transport Pulse Response Metrics Entire ensemble chemographs (i.e. b reakthrough curves) show substantial variation in transport pulse response (Figure 4 12). Entire ensemble results show high variability due to low spring mass flux for non behavioral cases. However, the comparison of the entire ensemble and behavioral group chemographs reveal highly consistent peak travel time of approximately 14 years (Figure 4 13), which wa s in good agreement with prior measurements of apparent groundwater age at Silver Springs tha t range from 9.9 to 27.5 years with a mean of 18.73 years [ Phelps and Survey , 2004] . Transport pulse response statistics provided further insights into differences b etween behavioral and non behavioral groups (Table 4 6). Overall mass recovery was good with a minimum of approximately 73 % and mean of 86 % mass recovery through both the surface and spring boundary for non behavioral cases and a minimum of 81 %
109 with a mean of 97 % for behavioral cases (Table 4 6 ) . Unrecovered mass likely remains in the domain at the end of the simulation, although model precision could account for up to 4% error at each boundary over the 1250 year simulation. Transport BTCs showed expected b ehaviors, behavioral replicates infiltrated much more mass into the domain illustrating the role that connected conduit networks play in aquifer vulnerability. Maximum, minimum and mean peak travel times , mean travel times , and travel time standard deviati on s were s maller in the behavioral cases. Behavioral and Non Behavioral Networks The MM GSA results (Figure 4 14) , show that hydraulic conductivity and HPF length scale (theta) are the most sensitive parameters for behavioral cases. The BMCF analysis exa mined differences between behavio ral and non behavioral groups. The KS test D parameter is the distance parameter, an indicator of differences between the empirical cumulative distribution functions between behavioral and non behavioral groups for each inp ut parameter. The magnitude of the D parameter depends on the number of samples; a smaller value indicates a higher probability that the samples in each group are drawn from the same statistical distribution, while larger values indicate a lower probabilit y that samples are drawn from the same distribution. The P value represents the probability that the two distributions are the same. This analysis show ed that hydraulic conductivity and the fractur e length scale parameter theta were the most sensitive para meters that distinguish between behavioral an d non behavioral realizations. Histograms of behavioral versus non behavioral paramet er s show that hydrau lic conductivity (p<0.001) and num ber of VPFs (p=0.006) were significantly different among behavioral and non behavioral groups, further underscoring the importance of these parame ters on network development. While distributions for
110 behavioral and non behavioral groups were significantly different for the number of VPFs, the number of VPFs did not show up as a sensitive parameter in the behavioral Morris sensitivity plot (Figure 4 14). This may b e due to the influence of non behavioral groups on steady state flow sensitivity . Results previously discussed indicated that the fracture length scale ( theta ) was also an important parameter in determining conduit network development, however it did not produce histograms that were different for behavi oral and non behavioral groups, likely due to high interactions of thi s parameter shown in Figure 4 15 . Differences in n etwork connectivity , such as the number of connections between HPFs were not explicitly described by the parameters so that may be driving some observed behaviors. Parameters with small Âµ values that are highly interactive such as spread m ay be exhibiting significant interactions with conduit geometry rather than the Morris parameter values. Network plots for non behavioral (Figure 4 16 ) an d behavioral groups (Figure 4 17 ) were used to further investigate causes of non behavioral groups. As expected, behav ioral groups had much more connectivity with the spring, whereas non behavioral networks tended to develop localized flow systems that discharge to the surface boundary. Figure 4 17 shows three behavioral conduit networks with the same peak solute travel time but increasing steady state flux from left to right. S teady state flux for Figure 4 17 A wa s 10 m 3 s 1 , Figure 4 17 B wa s 18 m 3 s 1 , and Figure 4 17 C wa s 24 m 3 s 1 . Th e se figures illustrate expected drivers of behavioral networks; intense development of conduits near the spring and fluxes which increase with network extent and number of fractures oriented along the dominant north south springshed head gradient s . These
111 networks also illustrate that in behavioral cases conduits tend to develop in the lo w lying areas of the model domain (Figure 4 17 , Figure 4 3 ). Behavioral steady state flow hydraulic conductivity end members show the influence of evolved conduit networks on steady state hydraulic heads . Figure 4 18 shows the lowest matrix hydraulic cond uctivity (1E 8 m s 1 ) . The low hydraulic conductivity leads to a highly spatially heterogeneous potentiometric surface. H ydraulic heads in Figure 4 18A we re significantly higher than observed in the Silver Springshed (Figure 5 1) owing to low matrix hydrau lic conductivity, however the patterns of high heads in the northern and southern part of the domain and wide flat area near the spring wa s consistent. The dark circular features represent local low head regions around the VPFs with similar hydraulic heads as the spring boundary condition. Lower hydraulic conductivity results in more pronounced differences between conduit influenced porous matrix cells and porous matrix cells more distant from conduits. Figure 4 19 shows the highest matrix hydraulic condu ctivity (1E 4 m s 1 ) . The hydraulic head values, pattern, and smooth potentiometric surface in the high hydraulic conductivity case 14 (Figure 4 19A) correspond well with observations from the Silver Springshed and illustrate the smoothing influence of hig h hydraulic conductivity on the potentiometric surface. Discussion The MM GSA results and BMC F indicate hydraulic conductivity strongly influences network development and that initial HPFs must have significant mean length scales, spread, and density to f acilitate the co nnectivity requisite to create a s ufficient drainage network to create a first magnitude spring . Th e relatively small number of behavioral replicates and the observed parameter sensitivities confirms the first
112 hypothesis that a limited rang e of parameters will result in a development of a first magnitude spring system and that parameter s that increase connectivity and decrease matrix permeability will be the most influential. The number and length of HPFs influence network drainage area, and connectivity to the spring . T he number of VPFs had to be significant in order to route calcite undersaturated water to the HPFs to facilitate network development. Behavioral Monte Carlo Filtering of MM GSA results show that the distributions between beha vioral and non behavioral parameters were highly similar. The HPFs and VPFs were randomly placed in the domain and interactions between network connectivity, HPF length scales and orientations, VPF locations, and porous media properties such as hydraulic c onductivity led to significant hydraulic and solute response differences between behavioral and non behavioral networks, with similar initial properties. Thus in addition to hydraulic conductivity and statistics of the HPF and VPF distributions, the actual connectivity of randomly placed HPFs and VPFs appears to have a large influence on whether a behavioral network evolves (Figure 4 14). This hypothesis will be further examined in the next chapter . Parameters zero we re likely influenced by non behavioral networks . For instance, sensitivity calculated for steady state flow was highly influenced by the significant difference in steady state flow between behavioral and non behav ioral groups where steady state flows were equal to almost all of the applied recharge for behavioral cases and almost none of the recharge for non behavioral cases.
113 The MM GSA results, steady flow behaviors, and BMCF results were used to examine what prop erties lead to development of a spring network and what timescales are appropriate to evolve each network for meaningful comparison. These results conf irmed that networks evolve until they reach a steady flow whe n connected pathways and head gradients stab ilize (Figure 4 5). After this point, conduit radii continue to grow but the springflow at t he outlet wa s steady . This steady flow regime represents when the contributing drainage areas to each HPF have reached a stable configuration. Comparisons between n etworks at the time of steady flow facilitated fair comparisons among networks . The MM GSA evaluated interactive and influential parameters of preferential flow path templates for hydraulic and solute pulse response in a hypothetical springshed. Hydrologic and transport pulse response metrics showed significant sensitivity to hydraulic conductivity and preferential flow path temp late properties (e.g. number of HPF, HPF length scales, spread, and numbers of VPFs). Porosity was expected to influence transport behavior by influencing velocity distributions in the porous matrix, however sensitivity to porosity was minimal. Behavioral and non behavioral groups had similar preferential flow path temp late pdf parameters, but substantially different hydrologic and t ransport response . This wa s another indication that the random configuration of these features and resulting connectivity and interactivity among parameters wa s a significant driver of response. Analysis of mean and standard deviation for hydraulic and t ransport pulse for behavioral netw orks showed mean steady state spring flows of 27.19 m 3 s 1 , a standard deviation of 3.70 m 3 s 1 (Table 4 4) and coefficient of variation (CV) of 0. 14 . Similarly
114 peak spring flows had a mean of 34.80 m 3 s 1 with a similar s tandard deviation of 6.21 m 3 s 1 with a CV of 0.18 . Transport response (Table 4 6), had a mean tracer peak travel time of 13.86 years with a CV of 0.28 , mean peak mass flux of 1.36 kg s 1 , with a CV of 0. 38 . Low CV values for hydraulic and transport metric s confirm the hypothesis that spring hydraulic and transp ort pulse response uncertainty wa s relatively low among behavioral replica tes. However, Figure 4 18, and Figure 4 19 illustrate that distributed aquifer behavior such as head fields were more highly variable and dependent on conduit configuration. To gain further insight into properties associated with behavioral networks, network plots were used to visualize differences between behavioral and non behavioral systems. Non behavioral plots showed poor c onnectivity with the spring and disconnected local flow systems (Figure 4 16 ). Large variability in resulting potentiometric surfaces were observed, even in behavioral replicates (e.g. Figure 4 18) with large variations in maximum hydraulic head values, pr onounced locally low head regions near VPFs, and pronounced undulations in the surface with lower matrix hyd raulic conductivity . Behavioral networks showed increases in steady state springflow with increasing spring connectivity and an increased distributi on of connected HPFs . I ncreased development of HPFs appeared to occur in the topographic low areas of the model domain that drain nearby high regions like th e western boundary where heads we re elevated relative to the spring (Figure 4 17 , Figure 4 18, Figu re 4 19). This last result has implications for ma nagement of springsheds because areas of higher conduit network de velopment correspond to areas of higher aquifer vulnerability. This support s the hypothesis that topographic and hydrologic features may be good predictors of
115 conduit development and thus may provide valuable information for mapping of vulnerable areas . Study Summary This study presented a framework for evolving random conduit networks in a hypothetical springshed representative of Silver Spri ngs (Figure 4 1). A generalized geologic model for local karst properties and boundary conditions was created using available literature, geologic reports, and insights from meetings with local hydrogeologists. These data were used to formulate parameter r anges and horizons for conduit development that were incorporated into a simplified regional scale hydrologic model derived from a widely used equivalent porous media regulatory hydrologic model. A preferential flowpath generation algorithm was developed t o create preferential flow path temp lates. These preferential flow path temp lates were subjected to dissolution by spatially distributed calcite undersaturated recharge to create karst aquifer realizations. Hydrologic and transport response in each replica te was used to determine whether the realization resulted in generating a significant spring network. Morris Method Global Sensitivity Analyses coupled with Behavioral Monte Carlo Filtering was used in an effort to identify porous media and preferential f low path parameter ranges more likely to result in a first magnitude spring. Furthermore, these analyses were used to determine what assumptions about initial paleo karst templates influence hydrologic and transport response. Variation in ensemble hydrolog ic and transport response was used to evaluate how uncertainty in conduit location and properties contribute to uncertainty in hydrographs, chemographs and ensemble statistics. This analysis resulted in the following insights:
116 A limited range of preferenti al flow path temp late parameters result in evolution of a conduit network that generates a first magnitude spring in the Silver Springshed During dissolution steady springflow occurs in the network and the network has reached a stable configu ration. After this time steady state sp ring flow does not change but conduit radii continue to grow, causing higher peak hydrograph response and faster recession to hydraulic pulses. L ower hydraulic conductivity and longer HPF length scales (theta) often resulted in mor e significant conduit development. Conduits tended to develop in topographic lows that drained nearby high regions. Steady spring flux was a function of number of conduit connections near the spring, overall network extent , and matrix hydraulic conductivit y . B ehavioral replicates occurred across all hydraulic conductivity ranges suggesting that connectivity is a limiting factor. HPF length scales, h ydraulic conductivity and number of VPFs and HPFs were generally the most sensitive parameters for hydrologic and transport pulse experiments. Behavioral and non behavioral networks often had similar preferential flow path temp late pdf parameters suggesting that random placement and resulting connectivity variation exerted a large control on behaviors of evolved preferential flow path temp lates. Parameter interactions among minimally sensitive parameters in the MM GSA results further support that the random behavior of networks influenced interpretation of input parameters. Hydrologic response for behavioral netw orks at t he time of steady flow resulted in mean steady state flow values of 23.78 m 3 s 1 with a standard deviation of 3.70 m 3 s 1 . Peak spring flows had a mean of 34.80 m 3 s 1 with a standard deviation of 6.21 m 3 s 1 . Transport response for behavioral net works at t he time of steady flow resulted in peak travel time with a mean of 13.86 year s and standard deviation of 3.83 years, peak mass flux values with a mean of 1.36 kg s 1 and standard deviation of 0. 51 kg s 1 , and mean travel times with a mean of 46 y ears and a standard deviation of 9.22 years. All of these results we re close to previously estimated mean groundwater age measurements reported in the literature [ Phel ps and Survey , 2004] .
117 Figure 4 1. Conduit Evolution Algorithm, adapted from Borghi et al., .
118 Figure 4 2. Silver springshed site location. A) Springshed location in Florida . B) Silver S prings, 1000 year spring groundwater capture zone, pho tolineament defined regional fracture sets, extensive sinks (black dots) and model domain.
119 Figure 4 3. Silver Spring model domain showing model grid and model surface topography elevations. Table 4 1. Horizontal and vertical preferential flowpath (HP F and VPF) and matrix properties for Morris Method Global Sensitivity Analysis. Parameter Units Probability Distribution * Source HPF Number C ount U(2500, 4000) [ Vernon , 1951] Northeast HPF Set 1 Orientation D e grees Fixed 45Â° [ Vernon , 1951] Northwest HPF Set 2 Orientation D egrees Fixed 315Â° [ Vernon , 1951] HPF Spread D egrees U(0Â°, 60Â°) Estimated HPF k U nitless U(0.75, 2.00) Estimated HPF theta M U(4000, 6000) Estimated VPF Number C ount U(250, 500) [ Denizman , 1998; Ford and Williams , 2013] Matrix Porosity U nitless U(0.25, 0.40) [ Langston et al. , 2012b] Matrix Hydraulic Cond. m s 1 LU(1.00E 4, 1.00E 8) [ Heath , 1983] Matrix Specific Storage m 1 U(1.00E 4, 1.00E 6) [ Batu , 1998] Epikarst Hydraulic Cond. m s 1 Fixed 1.00E 3 [ Heath , 1983] Epikarst Porosity U nitless 0.3 [ Langston et al. , 2012b] Epik arst Specific Storage m 1 Fixed 1.00E 4 [ Batu , 1998] *U(minimum, maximum) uniform distribution probability range for GSA LU(minimum, maximum) is log uniform dis tribution probability range for GSA
120 Figure 4 4. Example lnQ vs t plot showing hydrograph (black), log transformed flow (blue) and computed inflection point separating early and late hydrograph recession periods. Figure 4 5. Simulated springflow du ring dissolution . A) Steady flux dissolution time. B) One and a half times the steady flux dissolution time. Table 4 2. Comparison of flows and radii statistics at steady flux dissolution time and one and a half times the steady flux dissolution time. Steady state flow (m 3 s 1 ) Peak flow (m 3 s 1 ) Max radii (m) Mean radii (m) Median radii (m) Steady flux dissolution time Minimum 0.00 0.00 0.32 0.03 0.01 Maximum 28.67 43.88 14.26 3.48 4.29 Mean 18.45 23.60 3.91 1.24 1.26 Standard Deviation 13.0 5 17.02 2.71 0.97 1.18
121 Table 4 2 (Continued) Steady state flow (m 3 s 1 ) Peak flow (m 3 s 1 ) Max radii (m) Mean radii (m) Median radii (m) 1 .5 steady flux dissolution time Minimum 0.00 0.00 0.05 0.01 0.01 Maximum 28.68 0.00 0.35 0.03 0.01 Mean 18.63 46.95 15.26 3.72 4.59 Standard Deviation 13.18 25.25 4.19 1.32 1.35 Figure 4 6. Conduit flow during evolution. A) N etwork at 15,000 years . B) Network at 20,000 years. C) N etwork at 25,000 years. Table 4 3. Morris Method Global Sensitivity Analysis results showing behavioral and non behavioral group statistics for dissolution. Metric All HPF Density (km 1 ) Dissolved HPF Density (km 1 ) Median Radii (m) Max Radii (m) Effective Porosity Conduit Volume Fraction Behavioral Min imum 1.53 0.83 0.00 0.23 2.82 0.00 Maximum 6.11 5.72 0.00 4.29 14.26 0.01 Mean 3.75 2.55 0.00 1.83 5.30 0.00 Standard Deviation 1.13 0.94 0.00 1.01 2.03 0.00 Coefficient of variation 0.30 0.37 0.67 0.55 0.38 0 .17 Non Behavioral Minimum 1.53 0.00 0.01 0.32 0.25 0.00 Maximum 6.81 2.17 0.44 3.72 0.40 0.00 Mean 3.88 0.38 0.06 1.00 0.32 0.00
122 Table 4 3 continued Metric All HPF Density (km 1 ) Dissolved HPF Density (km 1 ) Median Radii (m) Max Radii (m) Effective Porosity Conduit Volume Fraction Standard Deviation 1.33 0.70 0.11 1.22 0.06 0.00 Coefficient of variation 0.34 1.82 1.78 1.23 0.18 1.75 Figure 4 7. Morris Method Global Sensitivity Analysis dissolution statistics showing Âµ, line) and Âµ to +/ 2 (blue lines) which aide in interpretation of parameter interactivity. A) D issolved HPF density . B) Median Radii. C ) M ax radii.
123 Figure 4 8. Morris Method Global Sensitivity Analysis hydrologic pulse response metrics . A) S teady state flow. B) Peak flow. C) Early recession coefficients. D) L ate recession coefficients.
124 Figure 4 9. Morris Method Global Sensitivity Analysis transport pulse response metrics . A) T ransport peak value . B) T ransport peak time . C) N ormalized first moment (mean breakthrough time) . D) N ormalized second central temporal moment (tracer variance) . Figure 4 10. Ensemble hydrograph showing mean (black line) and standard deviation (gray shaded region).
125 Figure 4 11. Behavioral ensemble hydrograp h showing mean (black line) and standard deviation (gray shaded region). Table 4 4. Morris Method Global Sensitivity Analysis results showing behavioral and non behavioral group hydrograph statistics. Metric Steady state f low (m 3 s 1 ) Peak f low (m 3 s 1 ) A lpha early recession c oefficient ( year 1 ) Inf lection t ime (days) Alpha late recession c oefficient ( year 1 ) Behavioral Minimum 23.78 33.73 0.23 2.03 0.00 Maximum 28.67 43.88 19.66 24.34 0.33 Mean 27.19 34.80 4.63 9.26 0.05 Standard Deviation 3.70 6.21 4.19 8.06 0.07 Coefficient of variation 0.14 0.18 0.90 0.87 1.42 Non Behavioral Minimum 0.00 0.00 0.00 0.00 0.00 Maximum 0.09 0.09 0.00 0.00 0.00 Mean 0.05 0.05 0.00 0.00 0.00 Standard Deviation 0.04 0.04 0.00 0.00 0.00 Coefficient of vari ation 0.77 0.77 0.00 0.00 0.00
126 Table 4 5. Behavioral Morris results showing behavioral and non behavioral group statistics for flow fractions through concentric control planes located 2, 8 and 12 km from the spring at steady state and peak flow. Met ric 2 km Steady State Flow Fraction 2 km Peak Flow Fraction 8 km Steady State Flow Fraction 8 km Peak Flow Fraction 12 km Steady State Flow Fraction 12 km Peak Flow Fraction Behavioral Minimum 0.99 0.99 0.98 0.98 0.97 0.97 Maximum 1.00 1.00 1.00 1.00 1.00 1.00 Mean 1.00 1.00 0.99 0.99 0.99 0.99 Standard Deviation 0.00 0.00 0.00 0.00 0.01 0.01 Coefficient of variation 0.00 0.00 0.00 0.00 0.01 0.01 Non Behavioral Minimum 0.00 0.00 0.00 0.00 0.03 0.03 Maximum 0.96 0.96 0.97 0 .97 1.00 1.00 Mean 0.34 0.34 0.45 0.45 0.61 0.61 Standard Deviation 0.40 0.40 0.33 0.33 0.38 0.38 Coefficient of variation 1.17 1.17 0.73 0.73 0.62 0.62 Figure 4 12. Ensemble chemograph showing mean (black line) and standard deviation (gray shaded region).
127 Figure 4 13. Behavioral ensemble chemograph showing mean (black line) and standard deviation (gray shaded region). Table 4 6. Behavioral Morris results showing transport metric summary for behavioral and non behavioral cases. Metric Percent Tr acer Mass Recovery Percent Tracer Infiltrated Transport Peak Time (years) Transport Peak Value (kg s 1 ) Mean Tracer Breakthrough (years) Tracer Standard Deviation (years) Behavioral Minimum 80.70% 27.28% 9.01 0.38 28.37 37.74 Maximum 99.70% 97.71 % 22.60 2.41 71.39 87.95 Mean 97.23% 76.76% 13.86 1.36 46.01 54.15 Standard Deviation 4.15% 20.78% 3.83 0.51 9.22 10.09 Coefficient of variation 0.04 0.27 0.28 0.38 0.20 0.19 Non Behavioral Minimum 72.92% 0.01% 5.96 0.00 60.94 85.30 Max imum 95.13% 0.25% 84.08 0.01 177.49 130.07 Mean 85.52% 0.12% 29.33 0.00 104.13 109.99 Standard Deviation 7.80% 0.09% 23.68 0.00 32.75 13.09 Coefficient of variation 0.09 0.73 0.81 0.83 0.31 0.12
128 Figure 4 14. Morris Method Global Sensitivity Analys is results showing sensitive and interactive parameters for behavioral networks.
129 Figure 4 15 . Behavioral Monte Carlo Filtering results for input parameters showing parameter distributions for behavioral (blue) and non behavioral (green), networks and st atistically significant parameters (red boxes).
130 Figure 4 16 . Two non behavioral conduit networks, showing development of localized flow systems. Figure 4 17 . Three behavioral conduit networks with the same transport peak time but increasing steady s tate flux from left to r ight. A) Steady state flux wa s 10 m 3 s 1 . B) Steady state flux wa s 18 m 3 s 1 . C) Steady state flux wa s 24 m 3 s 1 .
131 Figure 4 18 . B ehavioral example with spatially heterogeneous potentiometric surface. A) P otentiometric sur face. B) Conduit network.
132 Figure 4 19 . B ehavioral example with relatively smooth potentiometric surface. A) Potentiometric surface. B) Conduit network.
133 CHAPTER 5 EVALUATING THE IMP ORTANCE OF CONDUIT AND SINK LOCATION ON FLOW AND TRANSPO RT IN AN IDEALIZED SILVER SPRINGS MODEL Background Locations of karst conduits and conduit properties are unknown and poorly constrained . U ncertainty in conduit location, density, network geometry and connectivity on hydrologic response is not well quantif ied, leading to limited use of discrete continuum models that in corporate conduit networks for regional scale hydrolog ic regulatory models [ Reimann et al. , 2011b] . Nevertheless it is well documented that conduit networ ks influence flow and transport in karst aquifers considerably [ Bonacci , 1987; Worthington , 1999; Ford and Williams , 2013] . To date equivalent porous media (EPM) models have been typically used for management of water resources. These models increase matrix hydraulic conductivity to account for the influence of conduits. While this approach can provide reasonable water balances for calibrated models hydrologic flowpath heterogeneity, which is particularly important for predicting solute transport, is neglected . Furthermore, as shown in the previous chapter and by [ Bonacci , 1987] the majority of the flow in karst aquifers can occur in the conduit system . Therefore it is of prime importance to evaluate the uncertainty of conduit location, orientation and connectivity on flow and transport in karst aquifers. Study Obj ectives The objective of this study was to perform a Monte Carlo analysis of vertical preferential pathway and conduit network geometries to gain insight into potential flow and transport uncertainty in the idealized Silver Springshed. While conduit networ ks and flow and transport responses are simulated in an idealized system, interpretation of flow and transport phenomena in the idealized system should provide much needed
134 information about the importance of incorporating vertical prefe re ntial flow path an d conduit network uncertainty into hydrologic decision models. This study used insights obtained from the Morris Method Global Sensitivity Analysis from the previous chapter to examine hydrologic and transport pulse response in networks that generate sprin gflows and potentiometric surfaces that are representative of the Silver Springs system. Of the 90 simulations included in the Mor ris Global Sensitivity Analysis , only two networks had potentiometric surfaces, steady state flows [ Munch et al. , 2006] and peak transport times similar to prior published data for the Silver Springshed ( [ Phelps and Survey , 2004] , and St. Johns River Water Management District (F igure 5 1A) . While chapter 4 examined the influence of parameters governing the probability density function for horizontal preferential flowpaths (HPFs), density and angles of orientations for HPFs, density of vertical preferential flow paths (VPFs) and properties of the porous matrix, this study examined the uncertainty in hydrologic and transport pulse response to the actual spati al configuration of HPFs and VPFs in preferential flow path temp lates in a Monte Carlo framework. Chapter 4 behavioral Monte Carlo filtering results indicated that the length scale factor and hydraulic conductivity were highly influential in determining conduit network development and first order magnitude spring generation, but the mechanisms explaining differences between non behavioral and behavioral replicates were not fully uncovered by the Morris Method. It is hypothesized that variations in network connectivity imposed by the specific random placement and orientation of preferential flow paths, even for the same governing Morris Method parameters, led to differences
135 that contributed to behavioral network formation. Simply put, the random placement of preferential flow paths matters. Th is study asked two questions (1) for specified statistical properties of the preferential flow path temp late (HPF length, density, and orientation and number of VPFs) and porous media properties that are representative of the Silver Springshed, does random placement of individual HPFs and VPFs within the domain create differences in connectivity among replicates that drives significant differences in resulting conduit networks and the tendency to creat e first magnitude springs; and (2) what is the uncertain ty in head fields, and hydraulic and transport pulse response contributed by randomly placed VPFs and HPFs that honor statistical properties of the preferential flow path temp late and porous media properties representative of the Silver Springshed? To ans wer these questions, the statistical properties of the preferential flow path network and the porous media properties were fixed in a Monte Carlo experiment, so the only random variables are the placement and orientation of HPFs and the locations of VPFs. Differences between behavioral and non behavioral replicates were used to test the hypothesis that, in addition to the statistical parameters of the HPF length distribution and the porous media hydraulic conductivity (shown in chapter 4), the actual random placement of the preferential flow paths in the do main wa s an also imp ortant factor in determining whether a first magnitude spring develops, as well as the resulting behavior of head fields, spring hydrographs and chemographs. Methods Only t wo networks from the previous chapter showed head fields, spring flows and solute chemographs ( break though curves ) that were behavioral for the Silver
136 Spring system. Those networks were used as a basis to develop parameters and an additional behavioral constraint for defining whether each Monte Carlo replicate behaved like Silver Spring. To more rigorously define the additional hydraulic head behavior condition for a replicate to be categorized as behavioral for Silver Springs , the hydraulic heads of all transects for the previous chapter were examined . Simulations with hydraulic head values within 2 m of the 2010 potentiometric surface values where the four transects intersected the springshed boundary (Figure 5 1) were considered behavioral . Combined use of the first magnitude springflow and head metric s to the screen replicates included in the Morris Method ensemble in chapter 4 resulted in 2 Silver Springs behavioral replicates , 14 and 81. Th e se networks had very similar hydraulic properties (Table 5 1). However, ma ny of the preferential flow path temp lates with similar properties in the Morris Method ensemble were not behavioral, indicating the potential importance of VPF and HPF connectivity on eventual development of a spring network. Representati ve values near th e mean of the two behavioral networks were selected and fixed (Table 5 2), while the actual placement and orientation of HPFs and the location of VPFs was randomly varied . Monte Carlo replicates were generated and run through the dissolution, hydraulic pu lse and solute pulse experiments until the ensemble mean values for the flow and transport outputs of interest converged for the entire ensemble. Outputs of interest included: steady state flow, steady state hydraulic head transects from the spring though the model domain in the four cardinal directions (i.e. north, south, east, and west), hydraulic pulse peak flow, solute pulse peak travel time, solute pulse peak concentration, first normalized temporal moment of the solute breakthrough curve
137 (mean travel time), square root of the normalized second central temporal moment of the solute breakthrough curve (travel time standard deviation). Mean convergence was tested using a t test to test the null hypothesis that the ensemble means were equal. As discussed in chapter two, springs integrate springshed behaviors through converging flowpaths, and therefore solute chemographs from insta ntaneous surface solute pulses were typically smooth and did not reveal vulnerable areas of the springshed that contribute to e arly arrival of solute at the spring. Therefore in this study , a backwards transport pulse Monte Carlo experiment at the spring outlet was conducted to illustrate the spatially variability of solute mass loading, peak solute travel times, mean solute trave l times and thus vulnerable locations within the springshed. For the equivalent porous media (EPM) model and each behavioral conduit model replicate, a backwards pulse from the spring equal to the 1 kg m 2 over the entire springshed area was applied at t he spring. The velocity field was reversed, the simulation was run for 1250 years, and the breakthrough of solute at the model land surface was observed . Distributed plots of the 0 th moment (mass), peak travel time, and mean travel t ime at the land surface provide information on the spatial variability of vulnerable areas with in the domain. F or the behavioral ensemble , mean hydraulic heads and moments were computed and plotted for comparison. Results Convergence of ensemble output means occurred after 400 r eplicates. All P values for the hydraulic and transport pulse output metrics (Table 5 3) and head transects in the cardinal directions (Table 5 4) were greater than 0.7, allowing the acceptance of the null hypothesis that the means were equal. For the head transects,
138 only the means of the individual values along the transect and t statistics are provided in the table, however each location along each transect met the above criteria. Based on the results shown in chapter 4 it was expected that, because of th e high value for hydraulic conductivity required to simulate behavioral head fields, the number of behavioral replicates would be low, and in fact only 73 out of 400 replicates were behavioral for the first magnitude spring criterion . Only 37 of these were behavioral for the additional hydraulic head criterion. This confirms the hypothesis that in addition to HPF and VPF statistics prop erties and the hydraulic conduc tivity of the porous media (shown in the previous chapter), the actual random placement of H PF and VPF elements exerts significant control on behavioral conduit network development, and resulting flow and transport behaviors for behavioral and non behavioral networks . Hydrologic and Transport Response Hydrologic and transport metrics (Table 5 4) , hyd rographs (Figure 5 2), head transects (Figure 5 3), and chemographs (Figure 5 4 ) for the full ensemble showed high coefficients of variation (CV) , indicating substantial variability in dissolution, hydrologic, and transport behavior . For example, the CV for steady state flow was 1.86 , peak flow was 1.86 , peak solute travel time was 0.46 , and mean tracer breakthrough was 0.31 . Hydraulic pulse response statistics and hydrographs for the 37 behavioral replicates (Table 5 6 and Figure 5 5 ) showed that var iation among behavioral networks was very small , with a mean value of 28.66 m 3 s 1 with low standard deviation (<0.001 m 3 s 1 ), and extremely low CV (<0.001). Similarly, mean peak flow for behavioral replicates was 31.46 m 3 s 1 with a CV of 0.01. If a f irst magnitude spring developed in any
139 particular replicate, it tended to capture all of the applied recharge and springflows tended to the value of the areally integrated recharge . Mean behavioral head transects also showed small standard deviatio n with v alues of approximately less than approximately 1 m in each cardinal direction near the spring ( Figure 5 6 ). Non linear behavior of mean head gradient near the spring was observed , which wa s a diagnostic behavior for conduit dominated, rather than equivalen t porous media dominated flow systems at the regional scale [ Worthin gton , 2009] . Behavioral ensemble chemographs (Figure 5 7 ) also showed limited variability in solute in pulse response. Mean behavioral ensembl e peak travel time (19.95 years) was similar to previously measured mean groundwater ages of approximately 20 yea rs [ Phelps and Survey , 2004] . However, the highly skewed mean solute breakthrough curve resulted in much larger values for mean travel time (fi rst normalized temporal moment), approximately 57 years, and mean travel time standard deviation (square root of second centralized moment), approximately 73 years . Coefficients of variation of transport solute metrics we re very small; the CV of peak solute arrival time was 0.06 , the peak solute mass flux was 0.20, mean tracer br eak through was 0.05 , tracer spread ( standard deviation ) was 0. 08 . This low variability in steady flow, heads, hydraulic pulse breakthrough curves and solute pulse breakthrough curves across behavioral rep licates indicates that the actual spatial configuration o f the network wa s not extremely important to accurately predicting integrated spring measures within the set of behavioral replicates.
140 Backward Transport Pulse Monte Carlo A backwards transport pulse experiment at the spring outlet was conducted for the EPM model and the behavioral ensemble and to illustrate the spatial variability of solute mass loading, peak solute travel times, mean solute travel times and thus vulnerable locations within the sprin gshed. The backwards transport pulse experiment on the EPM model showed expected behaviors with smooth heads and relatively uniform contours (Figure 5 8 A) and relatively uniform 0 th moment at the surface (Figure 5 8B ) where mass from the spring backward tr ansport pulse arrives equally throughout the domain. The n otable exception wa s mass exiting the surface boundary in the northwest area of the domain. Peak arrival time and first normalized moment spatial maps (Figure 5 8 C and Figure 5 8 D) show concentric c ontours which illustrate that the EPM model, which neglects the influence of preferential flow paths, leads to uniform behaviors in backwards transport pulse response. The behavioral Monte Carlo ensemble mean hydraulic heads (Figure 5 9 A), 0 th moment (Figu re 5 9B ), peak arrival times (Figure 5 9C ), and 1 st normalized moment (Figure 5 9D ) show much more spatial variability than the EPM model. Backwards transport Monte Carlo e nsemble mean hydraulic heads reasonably reproduce the 2010 Silver Springshed potenti ometric surface. The 0 th moment shows most of the mass from the spring arrives in the relatively flat hydraulic head area near the spring with significantly less contribution from the high head areas at the north and south ends of the domain. The peak arri val time shows significant spatial variability; first arrival less than 10 years in areas spatially distant from the springs on the western edge of the domain and greater than 100 year peak arrival times at the distal edges of the north and
141 southern areas of the domain. The 1 st normalized moment shows that mean tracer breakthrough wa s less than 50 years near the spring, greater than 500 years in all of the high topographic areas, and 150 to 250 yea r mean breakthrough times that we re spatially associated wit h topographic lows. Discussion Differences in mean behaviors between the entire Monte Carlo ensemble and the behavioral ensemble were used to test the hypothesis that actual connectivity of preferential paths within a particular realization of the Monte Ca rlo experiment exert s important control on eventual development of a first magnitude spring. Since the ensemble statistical properties of every realization were the same, if random placement results in a significant number of both behavioral and non beh avi oral groups this hypothesis wa s supported. The small number of behavioral samples ( 37 out of 400) supports that the statistical parameters of the VPF and HPF distributions and the hydraulic con ductivity of the porous matrix we re not the only drivers , and p erhaps not the most significant driver s of first magnitude spring development. Connectivity of individua l realizations was found to be a significant driver of spring development and the consistent mean and standard deviation behavior among behavioral repl parameter values explored here. Recharge was either lost completely to the surface boundary or directed co mpletely to the spring outlet; the steady state spring flow was either <0.1 m3 s 1 (f or non behavioral replicates) or equal to total recharge applied (for behavioral replicates). Standard deviations for steady state flow, peak flow, peak arrival time, and peak mass flux were an o rder of magnitude higher for the entire ensemble than for the behavioral ensemble. Low standard deviations among behavioral replicate
142 metrics indicate that for the case s investigated here connectivity wa s extremely important in determining whether a spring w ill form, but once it formed knowledge of the actual locati on of preferential paths wa s not important for accurately predicting hydraulic and solute pulse response at the spring vent. The backwards tracer pulse experiment on the EPM model illustra tes that even when springflows we re equal to the recharge (28.65 m 3 s 1 ) and hydraulic head distributions we re representative of the Silver Springshed (Figure 5 8 A), the transport behaviors we re significantly different. The EPM model had relatively uniform mass arrival at the surface (Figure 5 8B ) and the peak arrival time (Figure 5 8C ) and first normalized moment (mean tracer breakthrough) (Figure 5 8D ) had concentric contours away from the spring. This finding illustrates the importance of including preferential flow paths for prediction of spatial patterns in springshed vulnerability. The backwards tracer pulse experiments illustrated that although the prediction of integrated spring vent hy draulic and transport response wa s not particularly uncertain, spatially variable differences vulnerability across the springshed exi st. The behavioral Monte Carlo ensemble hydraulic heads, peak arrival times, and 0 th and 1 st moments show significant variability (Figure 5 9). The mass distribution at the surface as a result of the backwards pulse (Figures 5 9B) indicated that all areas of the springshed were not equally vulnerable and did not respond equally to the solute injected at the spring. A higher percentage of solute mass exited the domain near the spring and in low topographic areas throughout the domain, while less solute mass exited at distant high topographic areas. Quick peak arrival times at areas spatially distant from the spring on the western edge of the domain further underscore how the
143 EPM model inadequately represented the distribution of vulnerability. Furthermore ove r larger timescales represented by the mean tracer breakthrough (1 st normalized moment), preferential flow paths in topographic lows exert significant influence on the breakthrough times, a finding not observable with EPM models or forward transport experi ments. In localized portions of the springshed, where the netwo rk wa s well developed, peak sol ute arrival times we re on the order of 0 40 years whereas in less well connected regions, typically in distant high topographic location s, peak solute arrival tim es were greater than 100 years . However total mass, peak travel time and mean travel tim e we re not well predicted by distance from the spring alone, topography and conduit network geometry interact ed to im pact springshed vulnerability. This emphasizes the importance of including conduit networks in vulnerability mapping exercise s , and s uggests that simulating transport using equivalent porous media flow and transport models is insufficient to characterize springshed vulnerability or plan effective land mana gement strategies for reducing contaminant concentrations in spring flows. Nevertheless the correspondence of the backward moment zeroth and peak travel time plots with the potentiometric surface plot indicate that if extremely high resolution maps of pot entiometric surface map were available they could be used as first order screening tools for vulnerable regions. To further examine the variability in behavioral replicates, two behavioral replicates were compared, cases 79 and 83. The HPF density for cas es 79 and 83 was representative of both non behavioral and behavioral networks (Figure 5 10 A, Figure 5 11 A). The hydraulic heads differed but showed the same spatial patterns with a low, flat
144 area in the center of the domain near the spring and high head a reas in the northern and southern portions of the domain (Figure 5 10 B, Figure 5 11 B). The resulting fracture networks after dissolution were different but have similar characteristics as previously observed behavioral replicates; extensive development of networks near the spring and long branches that extend up through the domain to the northern and southern boundaries (Figure 5 10 C, Figure 5 11 C). Figures 5 1 2 and 5 1 3 show results of the backward pulse responses for case 79 and 83 respectively. Figure 5 1 4 shows the breakthrough curve at the spring for forward solute pulse experiments on each of these replicates. The forward pulses were smooth and resulted in different solute mass contributions but showed similar peak travel times to the spring. Peak sol ute travel times at the surface as a result of the backward pulse (Figure 5 1 2 A and 5 1 3 A) illustrate the importance of VPFs and the conduit network configuration to determining vulnerable areas of the springshed that quickly contribute solute to the sprin g outlet. For both replicates, the western side of the domain, spatially distant from the spring, showed the fastest peak pulse arrival times (i.e. less than 10 years). In contrast, high topographic areas in the northern and southern part of the domain sho wed much longer peak travel times even where extensive conduit network development occurred beneath them. The spatial map for distribution of mass exiting the surface of the domain (Figures 5 1 2 B and 5 1 3 B) showed that areas contributing significant mass t o the spring were similar for the two replications and corresponded quite well on to the piezometric surface of the individual replicates (Figures 5 10 A and 5 11 A). While there was spatial variability within and across the two replicates the highest mass w as recovered near the spring in both cases, and less mass exited the
145 domain in areas with higher topography distant from the spring in both cases. The spatial map for distribution of mean arrival time (Figure 5 1 2 C, Figure 5 1 3 C) showed that the mean trave l time wa s highly influenced by conduit geometry with exit times less than 50 years occurring in very close proximity to the developed conduit networks. Together these data show that for behavioral network s vulnerable regions we re spatially variable but qu ite similar between the replicates and the ensemble mean behaviors . Study Summary Hydrologic models used for karst aquifer management typically neglect the influence of preferential flow pathways on hydrologic and transport response due to lack of informa tion about the location and properties of preferential flow pathways. This study has shown that in addition to statistical properties of vertical and horizontal preferential pathways and the hydrologic properties of the porous media, the actual location of VPFs and HPFs in relation to each other and the spring outlet determines whether a spring will develop in a karst landscape. However , for the Monte Carlo simulation conducted in this study , if a network developed the uncertainty in hydraulic and solute pu lse response at the spring vent due to unknown locations of network elements (e .g . VPFs and HPFs) was minimal. This indicates that incorporating preferential flow processes in to the model may be more important than exact knowledge of preferential flowpath locations and orientations for understanding spatiotemporally integrated behaviors at springs. However backwards transport pulse experiments confirmed that conduit networks do contribute significant spatial heterogeneity in tracer travel times, surface vul nerability, and long term breakthrough tailing behavior.
146 Figure 5 1 . Potentiometric surface maps for Silver Springs. A) C ontoured potenti ometric surface in ft. for 2010. B) A dapted equivalent porous media (EPM) NDM model hydraulic heads for layer 4 . Table 5 1. Hydrologic and transport pulse metrics for behavioral cases from Chapter 4 with values similar to published values for Silver Springs. Peak f low (m 3 s 1 ) Steady state f low (m 3 s 1 ) Transport peak t ime (years) Transport peak v alue (kg s 1 ) C ase 14 30.00 28.90 18.44 0.57 Case 81 34.97 28.94 17.54 1.00 Table 5 2. Behavioral network properties from cases similar to Silver Springs from c hapter 4 and selected parame ters for c hapter 5 Monte Carlo e xperiments. K Theta Number of VPF Number of H PF Hydraulic c onductivity (m s 1 ) Spread (Â°) Porosity Specific s torage (m 1 ) Case 14 1.17 6000 500 3000 0.0001 20 0.25 3.4 E 05 Case 81 1.58 4000 200 3500 0.0001 40 0.35 6.7 E 05 Average 1.375 5000 350 3250 0.0001 30 0.3 5.1E 05 Selected p arame ters 1.5 5000 350 3250 0.0001 30 0.3 1.0 E 05
147 Table 5 3. Output statistic convergence results for Monte Carlo simulations. Output Statistic Mean 300 r eplicates Me an 400 r eplicates T statistic P value Steady state f low (m 3 s 1 ) 6.47 6.36 0.17 0.90 Peak f low (m 3 s 1 ) 6.67 6.56 0.39 0.73 Transport peak t ime (years) 18.88 18.56 0.31 0.79 Transport peak v alue (kg s 1 ) 0.13 0.15 0.26 0.76 1st normalized m oment (years) 71 .57 70.37 0.11 0.93 2nd normalized m oment 2 (years) 94 .24 92.66 0.39 0.73 Table 5 4 . Head t ransect output statistic conve rgence results for Monte Carlo s imulations Table 5 5. Complete ensemble m ean and standard deviation for h ydrologic and transport metrics. Statistic Steady state f low (m 3 s 1 ) Peak f low (m 3 s 1 ) Transport peak t ime (years) Trans port peak v alue (kg s 1 ) Mean tracer b reakthrough (years) Tracer s tandard deviation (years) Min 0.00 0.00 5.52 0.00 47.31 61.78 Max 28.67 31.46 82.27 1.03 209.68 143.01 Mean 6.36 6.56 18.56 0.15 70.37 92.66 Standard d eviation 11.85 12.22 8.49 0.28 21. 53 23.69 Coefficient of variation 1.86 1.86 0.46 1.89 0.31 0.26 Output s tatistic Mean 300 r eplicates Mean 400 r eplicates Mean T statistic Mean P value East transect 19.07 18.96 0.37 0.75 North transect 21.10 24.08 0.43 0.71 South tra nsect 19.88 19.39 0.39 0.70 West transect 19.62 19.45 0.31 0.73
148 Figure 5 2 . Complete ensemble hydrograph showing mean (black line) and standard deviation (gray shaded region). Figure 5 3. Complete ensemble head transects showing mean (black line) and standard deviation (gray shaded region). Spring is at zero on x axis and distances are in the direction away from the spring. A) North transect. B) West transect. C) South transect. D) E ast transect.
149 Figure 5 4 . Complete ensemble chemograph show ing mean (black line) and standard deviation (gray shaded region). Table 5 6. Behavioral ensemble mean and standard deviation for Hydrologic and transport metrics. Statistic Steady state f low (m 3 s 1 ) Peak f low (m 3 s 1 ) Transport peak t ime (years) Trans port peak v alue (kg s 1 ) Mean tracer b reakthrough (years) Tracer standard deviation (years) Min 28.66 29.02 16.09 0.35 51.63 61.92 Max 28.67 31.46 22.88 1.03 67.61 93.43 Mean 28.66 29.57 19.95 0.66 56.83 73.15 Standard Deviation 0.00 0.43 1.18 0.14 3. 01 5.97 Coefficient of variation 0.00 0.01 0.06 0.20 0.05 0.08
150 Figure 5 5 . Behavioral ensemble hydrograph showing mean (black line) and standard deviation (gray shaded region). Figure 5 6. Behavioral ensemble head transects showing mean (black lin e) and standard deviation (gray shaded region). Spring is at zero on x axis and distances are in the direction away from the spring. A) north transect, B) west transect, C) south transect, and D) east transect.
151 Figure 5 7 . Behavioral ensemble chemogra ph showing mean (black line) and standard deviation (gray shaded region).
152 Figure 5 8 . Equivalent porous media (EPM) model hydraulic heads and moments . A) Hydraulic heads. B) 0 th moment. C) Peak arrival time. D) 1 st normalized moment.
153 Figure 5 9 . Behavioral Monte Carlo ensemble model hydraulic heads and moments . A) Hydraulic heads. B) 0 th moment. C) Peak arrival time. D) 1 st normalized moment.
154 Figure 5 10 . Behavioral case 79 . A) Initial HPF network. B) Steady state hydraulic heads. C) D eveloped conduit network. Figure 5 11 . Behavioral case 83. A) Initial HPF network. B) S teady state hydraulic heads . C) D eveloped conduit network.
155 Figure 5 1 2 . Backwar ds transport pulse for case 79 . A ) M ass recov ery at the surface (0 th momen t). B ) P eak arrival time . C) M ean breakthrough time (1 st normalized moment ). D) M odel surface topography. Figure 5 1 3 . Backwar ds transport pulse for case 83 . A) M ass recov ery at the surface (0 th moment). B) P eak arrival time . C) Mean breakthrough time (1 st normalized moment). D) M odel surface topography.
156 Figure 5 1 4 . Forward transport pulse for case 79 and 83 showing tracer mass flux.
157 CHAPTER 6 SUMMARY AND CONCLUSIONS This research utilized a multi disciplinary approach to examine spatially distributed evidence of nitrate transformations, in situ nitrogen transformations , and the influences of geologic heterogeneity on conduit development and flow and solute transport fluxes and flowpaths in karst terrain. Spatially Distributed Nitrate Tran sformations The first element of the research examined the spatial distribution of denitrification proxies to uncover relationships between local aquifer geochemistry, aquifer confinement, and spatial flow focusing karst features associated with nitrat e re moval via denitrification. This study found that aquifer confinement was not associated with observed denitrification evidence despite local conditions which we re favorable for denitrification (e.g. low DO, negative ORP, organic (DOC) and solid phase (iron and manganese) electron donor availability. All of the measured evidence of denitrification occurred in the unconfined region of the springshed. This does not necessarily mean that denitrification is not occurring in the conf ined region but rather that de nitrification was not observed with the density of sampling of the confined region . In the confined region, surficial flow paths towards depressions, sinks and wetlands that may b e connected to the UFA may create nitrate transformation hotspots over smal l spatial scales not captured by the limited density of available wells to sample in the unconfined region. Therefore, the limited number of available sites in the confined region may not have captured enough of the variability in geochemistry to make the assertion that denitrification is not occurring there.
158 In general, local geochemical conditions were a poor predictor for denitrification evidence. It is likely that local scale flow path variability inherent to karst terrain that can may confound interpr etation of data at individual sites. Observed excess nitrogen gas may have been generated upgradient where conditions were move favorable and delivered to the site without significant dilution from other flow paths. Another possible reason for the poor cor relation between local conditions and denitrification is that denitrification may be occurring in localized microsites on biofilms within aquifer material. Significant evidence of biofilms in previous cave studies suggest that this biofilm denitrification mechanism may be a more probable explanation. Relationships between excess nitrogen gas (i.e. direct evidence of denitrification) and flow focusing mechanisms (e.g. sinks, swallets, closed topographic depression, and wetlands) were evaluated. Of all featur es examined, only proximity to wetlands had a significant correlation with excess nitrogen gas. However this relationship only held for sites w here excess nitrogen gas was detected suggesting that confounding factors such as wetland aquifer connectivity an d upgradient location within the hydrologic flow path to the sampling site may contribute to the lack of relationship between excess N 2 denitrification evidence and distance to wetland metrics . This finding presents an interesting research topic for future investigation : more rigorous evaluati on of denitrification proxies within the aquifer beneath wetlands, can be used to confirm or refute this potential relationship . Lastly, the correspondence between spatially distributed denitrification measures and sp ring vent measurements was evaluated. Spring vents measurements have long been used to characterize karst springsheds however how those measures relate to
159 processes distributed throughout springsheds is still poorly understood. Groundwater samples from wel ls represent only a spatial average of the well contributing area while spring vent measures provide flow averaged values across the entire springshed due to the inherent integration of flow paths at spring vents. This study found that local scale heteroge neity in N transformations (e.g. nitrification) confounds expected isotopic signatures. Spatially averaged excess nitrogen gas measures were close to spring vent measurements , 0.31 mg L 1 and 0.51 mg L 1 respectively, however the expected linear relationsh ips between DO and increasing 15 N were not observed in the groundwater samples from wells as they have been in integrated spring measures. Together these data suggest that spring vent measurements remain the most cost effect i ve way to characterize springshed denitrification because ob served heterogeneity in spatially distributed N tra nsformations would require significant sampling densit ies . In Situ Denitrification Measurements in Geochemical End Members The second study directly measured rates and mechanisms of N transformation in th e aquifer across a redox gradient : at a completely oxic and highly anoxic site representing end members of redox state in wells sampled for the first study. Push pull tracer tests with nitrate and nitrate plus carbon addition were combined with microbial g enetic data analyses to provide information about aquifer response to nitrate and carbon addition and to identify important transformation mechanisms. D ifference s between reactive and conservative tracer mass recovered in push pull tracer tests were used t o compute zero order loss rates. B aseline conditions and response to each tracer injection were examined using microbial DNA and rRNA which identified both the genes in the microbial communities as well as the genes that were actively being used.
160 Significa nt nitrate mass was still missing after an N balance was computed for each experiment which measured all aqueous N forms, and estimated maximum N assimilation by microbes using microbial population data , indicating potential loss by denitrification . Denitr ification genes were detected as being present in the communities but genes in sampled water were not active. This further supports the assertion that denitrification may be occurring in biofilms on aquifer material that were not sampled in this study. How ever the reported rates represent zero order (mass time 1 ) nitrate losses, instead of denitrification rates because measurements represent unrecovered nitrate and the portion of nitrate loss attributed specifically to denitrification could not be determine d. Nitrate losses were measured both with and without carbon addition. At the anoxic site without carbon addition, increased ammonium concentrations and genetic data support that dissimilatory reduction of nitrate to ammonium occurred. This is significant because this mechanism of nitrate reduction has not been documen ted at all in the UFA and only one study has reported it in freshwater aquifers [ Bengtsson and Annadotter , 1989] . With carbon addition, the mechanism of nitrate loss shifted and DNRA became less significant as evidenced with genetic data and reduced ammonium concentrations . T he mechanism in this case was likely denitrification due to enhanced carbon availability, although porewater genetic data could not confirm it was occurring the p orewater. This supports a recurrent finding in this study : unexplained nitrogen losses o ccur s in zones where conditions a re not considered to be favorable (e.g. high dissolved oxygen) , leading to the inference that biofilm processes on aquifer materials ma y be significant.
161 Zero order nitrate loss rates with carbon addition were much higher at the anoxic site ( 0.13 +/ 0.02 mmol L 1 hr 1 ) than the oxic site (0.05 +/ 0.01 mmol L 1 hr 1 ). Together th ese data indicate that nitrate wa s not conservatively transp orted in the UFA even when local conditions are considered un favorable for denitrification. While the measured rates we re consistent with literature values from similar experiments in basin fill aquifers, these rates we re high relative to the documented ni trogen loads to the springshed which ranged from 262,000 to 1.3 million kg year 1 [ Katz et al. , 2009] . Assuming that the nitrate loads for each year are completely mixed in only the top 1 % of the aquifer volume estimated from the thickness of the Ocala Limestone ( 0.58 m), a matrix porosity of 0.3, a nd the springshed area of 960 km 2 results in nitrate load 1. 32 e 5 to 6 .71e 5 mmol L 1 hr 1 , which is 3 orders of magnitude less than the measured nitrate loss rates. If this were the case nitrate contamination would not be an issue. These rates may include other short term sinks for nitrate such as assimilation and may be influenced by the high nitrate concentrations in injected tracers. Regardless, they provide a relative measure of denitrification potential and these experiments confirmed that denitrifica tion is not the only mechanism of nitrate transformation in the UFA. Observed nitrate reduction to ammonium, unexplained nitrate losses, and observed assimilation, s upport that nutrient spiraling is occurring in the UFA as water passes through areas of the aquifer with different geochemical conditions. Effects of Karst Aquifer Complexity on Fluxes and Flow Paths A stochastic framework was used to generate preferential flow path temp lates which were evolved into conduit networks for a hypothetical springshe d. Calcite dissolution and evolution of these preferential flow path temp lates was evaluated to gain insight into parameters that lead to development of a spring network. These networks
162 were placed in a hydrologic flow and transport model and used to evalu ate what preferential flow path temp late properties influence hydrologic and transport response metrics. T his study addressed three questions; what initial assumptions about preferential flow path temp lates lead to development of first magnitude springs , h ow long do you evolve networks for to provide useful comparisons across replicates, and what influence do these evolved preferential flow path temp lates have on springshed behaviors such as hydraulic head, springflow and solute transport. The most importa nt preferential flow path temp late parameters were the length scale (theta) and density of the horizontal preferential pathways (HPFs) and vertical preferential pathways (VPFs), and the variation in orientation. All of these factors contribute to network c onnectivity. The most important porous matrix property for conduit evolution, and hydraulic and transport pulse response was hydraulic conductivity. H igher porous matrix hydraulic conductivity may impede network development by decreasing head gradients, an d allowing diffuse recharge to be distributed throughout the domain in the matrix instead of focusing undersaturated water into preferential pathways connected to springs. Network development is therefore a competition between hydraulic conductivity of the matrix and connectivity of the network. Hydraulic and transport pulse response experiments showed high variability due to the large contrast between paleo karst templates that developed into a first magnitude spring and those that did not. However variati on in hydrologic and transport spring response among net works where springs develop was much lower . Evaluation of par ameter ranges using behavioral M onte Carlo filtering showed that the distributions of
163 parameters for behavioral and non behavioral samples were similar suggesting that the specified properties may not be the only driver s of non behavior. Effects of Randomly Distributed Preferential Flow P aths. Selected parameter ranges in the Morris Method study described above only generated behavioral netw orks approximately 6 6 % of the time indicating that there was another process influencing network development whose influence was not well described . This raised the question as to whether non behavioral networks were non behavioral only because of the inte raction among paleo karst template statistical properties and porous media properties, or whether the actual random placement and orientation of prefer ential flow paths mattered. Two replicates from the prior study which behaved similar to the Silver Sprin g system were selected. Preferential flow path temp late statistical properties and porous media properties were fixed and Monte Carlo simulations were conducted in which the only random behavior was the placement of HPFs and VPFs within the domain. O nly 83 of the 400 Monte Carlo realizations resulted in spring s, and these springs essentially captured all the applied recharge suggesting that specifi c network connectivity details we re an important driver of network development. Moreover, only 37 of 400 also h ad reasonable potentiometric surfaces comparable to Silver Springs. The low varia bility in spring flow and transport behavior among behavioral replicates indicated that , if only sprin g vent predictions are required, knowledge of exact conduit location is p erhaps not as important as long as preferenti al flow paths are represented, h owever a backwards tracer pulse injection Monte Carlo expe riment revealed that different conduit networks lead to spatially vulnerable regions in the springshed, an outcome not di scernable from only integrated signals measures at spring vents.
164 Insights for Further Research This study examined the role of hydrogeological heterogeneity, assumptions about preferential flow paths, and drivers of denitrification and other nitrogen trans formation mechanisms at local and springshed scales . Insights for further research gained from this study include : Denitrification was often observed where local conditions were unfavorable, in both studies where direct measurements were made. This indicat es that local scale nitrate transformations may be occurring in aquifer biofilms. Incorporation of biofilm sampling in further research may help to ascertain its potential importance. Denitrification evidence showed a suggestive correlation with nearby wet land area, number of wetlands and distance to wetlands. Further nitrogen transfo rmation research may benefit fro m a more rigorous evaluation of this wetland aquifer denitrification connection. This study uncovered evidence of nutrient cycling that should b e evaluated in future studies of nitrogen transformations in the UFA. Quantifying the spatial distribution of zones of similar character to the site where DNRA occurred in the aquifer will provide deeper insights into nitrate cycling and denitrification. F urthermore , sampling of aquifer biofilms will help determine their potential influence on regional scale nutrient spiraling in the UFA.
165 APPENDIX A SELECTED GOECHEMICAL DATA FOR THE ICHETUCKNEE SPRINGSHED Table A 1. Measured geochemical parameters for Ph ase I and Phase II sampling in the Ichetucknee Springs Basin. Site number Sample date pH DO ( mg L 1 ) ORP (mv) SPC (ms cm 1 ) K ( mg L 1 ) Na ( mg L 1 ) Mg ( mg L 1 ) F ( mg L 1 ) Cl ( mg L 1 ) SO 4 ( mg L 1 ) NO x N ( mg L 1 ) 2 8/9/2013 7.34 7.34 116.9 276 0.52 4 .45 4.62 0.00 7.25 2.07 1.14 3 8/9/2013 7.94 0.36 379.3 89 1.02 0.74 1.34 0.00 1.06 0.00 0.10 4 8/2/2013 7.84 7.82 75 854 0.64 2.07 1.22 0.00 3.16 4.18 2.63 5 7/31/2013 6.67 6.81 158.6 485 0.09 3.67 1.22 0.00 5.27 0.32 0.33 6 8/2/2013 8.46 0.44 119.6 440 0.50 4.10 5.59 0.00 6.08 1.93 0.31 7 8/8/2013 7.61 0.94 92.7 393 0.52 3.32 5.83 0.08 4.75 4.96 0.12 8 8/1/2013 6.98 1.82 144.7 409 0.53 5.62 8.99 0.04 7.21 19.98 0.41 9 8/2/2013 8.42 0.15 177.3 761 0.42 2.71 3.16 0.00 3.84 8.16 0.07 10 8/8/2013 6.32 5.27 41.6 669 1.42 2.12 1.94 0.00 4.01 28.93 4.39 11 8/8/2013 7.39 1.12 16.1 655 0.87 3.78 10.70 0.12 5.03 51.22 0.10 12 8/1/2013 6.83 2.37 118.4 409 0.20 2.60 5.59 0.00 4.45 5.12 0.69 13 8/1/2013 7.13 5.74 155.5 420 0.50 5.00 9.48 0.07 7.60 30.82 0.05 14 8/2/2013 7.18 7.72 153.6 298 0.12 1.78 1.22 0.00 3.67 4.04 0.20 15 8/2/2013 6.89 6.05 110.8 567.5 0.25 2.25 1.94 0.00 3.41 3.75 0.07 16 8/2/2013 7.26 7.49 89.9 928 0.31 2.23 2.67 0.00 2.78 3.55 0.11 1 3/18/2014 8.69 0.66 191.4 291 22.43 9.87 12.70 0.36 5.70 4.72 0.02 2 3/9/2014 7.25 8.63 143.4 273 0.50 4.33 4.38 0.21 7.57 1.85 1.20 3 3/18/2014 8.07 0.32 396 84 0.95 0.75 1.64 0.16 0.99 0.00 0.03 5 3/7/2014 6.64 6.81 262.4 524 0.06 3.70 0.82 0.06 6.05 0.48 0.38 8 3/9/2014 6.86 3.33 176.6 46 2 0.47 4.93 9.13 0.18 8.07 17.57 0.39 9 3/14/2015 6.91 0.53 141 328 0.23 2.54 3.04 0.09 4.07 8.53 0.12 10 3/14/2014 6.63 6.06 119.1 446 1.24 1.68 1.80 0.07 3.26 27.42 3.68 12 3/9/2014 6.76 3.01 176.5 404 0.22 2.95 5.54 0.13 4.97 5.08 0.78 13 3/7/2014 6.97 4.4 173 424 0.41 4.42 9.72 0.15 7.47 25.33 0.04 14 3/9/2014 7.1 7.09 178.6 297 0.14 1.78 0.85 0.07 3.52 3.97 0.37 15 3/7/2014 6.96 5.7 141.6 344 0.22 2.15 1.61 0.04 3.22 3.46 0.06 16 3/14/2014 6.7 8.02 154 393 0.12 2.12 1.99 0.06 2.70 3.47 0.03
166 T able A 2. Nutrients and dis solved carbon measurements for P hase I and Phase II sampling in the Ichetucknee Springs Basin Site number NH 3 ( mg L 1 ) TKN ( mg L 1 ) PO4 (ug L 1 ) DOC ( mg L 1 ) DIC (ugC g 1 ) 2 0.08 0.08 0.03 0.00 27.69 3 0.19 0.21 0.00 1.53 8 .05 4 0.08 0.08 0.03 0.29 40.15 5 0.10 0.01 0.05 1.23 66.70 6 0.08 0.04 0.00 0.35 35.33 7 0.09 0.04 0.04 0.58 29.31 8 0.07 0.02 0.05 0.41 44.27 9 0.08 0.05 0.01 0.26 36.78 10 0.08 0.05 0.07 1.17 44.31 11 0.36 0.47 0.06 1.42 44.45 12 0.07 0.02 0.06 0.17 48.49 13 0.09 0.07 0.04 0.30 36.74 14 0.08 0.05 0.14 0.06 34.33 15 0.07 0.08 0.10 0.20 41.75 16 0.07 0.06 0.06 0.12 50.29 1 0.44 0.74 0.04 2.71 31.28 2 0.03 0.10 0.02 0.00 26.18 3 0.03 0.14 0.00 1.87 8.39 5 0.03 0.12 0.04 0.63 66.70 8 0.04 0 .18 0.05 0.41 50.26 9 0.03 0.11 0.01 0.70 39.47 10 0.05 0.20 0.06 1.17 46.85 12 0.03 0.14 0.04 0.54 44.91 13 0.03 0.04 0.04 0.55 41.74 14 0.05 0.10 0.12 0.54 34.33 15 0.03 0.32 0.09 0.70 41.75 16 0.30 0.05 0.06 0.56 52.43
167 Table A 3. Dissolved met al concentrations for Phase I and Phase II sampling in the Ichetuck n ee Springs Basin. Site number Fe ( mg L 1 ) Mn ( mg L 1 ) Sr ( mg L 1 ) V ( mg L 1 ) Cr ( mg L 1 ) Co ( mg L 1 ) Ni ( mg L 1 ) Cu ( mg L 1 ) Zn ( mg L 1 ) 2 0.0011 0.0002 0.0702 0.0022 0.0019 0.00 00 0.0000 0.0011 0.0052 3 0.2200 0.0606 0.9234 0.0000 0.0005 0.0000 0.0000 0.0000 0.0007 4 0.0003 0.0000 0.0631 0.0020 0.0015 0.0000 0.0000 0.0005 0.0388 5 0.0039 0.0005 0.0857 0.0016 0.0009 0.0008 0.0004 0.0029 0.0059 6 0.2628 0.0078 0.0589 0.0020 0.0 000 0.0000 0.0000 0.0001 0.0602 7 0.0791 0.0084 0.0554 0.0015 0.0001 0.0001 0.0002 0.0003 0.1044 8 0.0002 0.0001 0.2393 0.0021 0.0004 0.0001 0.0000 0.0018 0.0085 9 0.4632 0.0377 0.1218 0.0000 0.0003 0.0000 0.0002 0.0001 0.0026 10 0.0025 0.0014 0.1468 0 .0075 0.0009 0.0001 0.0002 0.0004 0.0066 11 0.0781 0.0277 0.0678 0.0002 0.0000 0.0002 0.0008 0.0005 2.5932 12 0.0028 0.0020 0.1308 0.0016 0.0010 0.0001 0.0001 0.0039 0.0185 13 0.0000 0.0001 0.3148 0.0005 0.0002 0.0000 0.0000 0.0011 0.0069 14 0.0002 0.0 001 0.0676 0.0022 0.0022 0.0001 0.0000 0.0036 0.0080 15 0.0020 0.0001 0.0882 0.0038 0.0013 0.0001 0.0002 0.0052 0.0196 16 0.0009 0.0002 0.1360 0.0011 0.0016 0.0002 0.0001 0.0042 0.0063 1 0.0068 0.0142 0.0904 0.0002 0.0009 0.0004 0.0006 0.0002 0.0012 2 0.0013 0.0002 0.0828 0.0027 0.0023 0.0000 0.0001 0.0016 0.0005 3 0.2852 0.0624 0.9941 0.0000 0.0004 0.0001 0.0002 0.0000 0.0020 5 0.0115 0.0006 0.1032 0.0019 0.0010 0.0007 0.0014 0.0041 0.0200 8 0.0006 0.0001 0.2408 0.0019 0.0005 0.0001 0.0001 0.0020 0. 0058 9 0.5262 0.0353 0.1273 0.0000 0.0000 0.0000 0.0001 0.0000 0.0017 10 0.0042 0.0007 0.1432 0.0061 0.0011 0.0001 0.0008 0.0001 0.0014 12 0.0009 0.0004 0.1426 0.0018 0.0013 0.0000 0.0001 0.0030 0.0161 13 0.0015 0.0003 0.4017 0.0008 0.0005 0.0000 0.000 1 0.0017 0.0104 14 0.0006 0.0002 0.0750 0.0025 0.0021 0.0001 0.0002 0.0039 0.0105 15 0.0097 0.0006 0.1086 0.0048 0.0015 0.0001 0.0012 0.0035 0.0193 16 0.0012 0.0001 0.1426 0.0011 0.0017 0.0001 0.0001 0.0027 0.0076
1 68 Table A 4 . Geospatial statistics for Ichetucknee Springs Basin spatial analysis. Site number Confining unit thickness (ft) Number of CTD within 0.5km Nearest CTD (km) CTD Area 0.5km Number of Wetlands Within 0.5km Nearest wetland (km) Wetland area 0.5km Number of swallets within 0.5km Neares t swallet (km) 1 125 0 0.69 0.00 2 0.06 0.03 0 4.54 2 65 4 0.27 0.30 0 0.93 0.00 0 4.46 3 68 0 0.87 0.00 3 0.40 1.38 1 0.87 4 24 3 0.05 0.81 3 0.13 0.01 1 0.16 5 28 6 0.04 0.27 1 0.49 0.00 0 4.10 6 18 5 0.12 1.00 4 0.38 0.02 2 0.55 7 48 4 0.13 0.34 2 0.32 0.01 0 1.78 8 2 5 0.10 0.13 0 1.61 0.00 1 1.01 9 7 6 0.06 0.03 2 0.20 0.02 1 0.06 10 2 3 0.04 0.63 2 0.22 0.01 0 2.74 11 79 0 0.68 0.00 4 0.14 0.03 0 4.05 12 6 3 0.14 0.02 0 1.77 0.00 1 0.50 13 20 6 0.00 0.08 0 1.93 0.00 0 2.47 14 41 1 0.43 0 .00 0 2.53 0.00 0 3.79 15 12 9 0.06 0.51 0 0.90 0.00 0 1.97 16 8 5 0.27 0.09 1 0.25 0.00 1 0.25
169 Table A 5. Si te names and nitrate isotope data for Phase I and Phase II sampling in the Ichetucknee Springs Basin , showing s ites with non detect ion (ND) where insufficient nitrate was available for analysis and s ites where isotope samples were not measured (NM) . Site number Site name Phase Date 15 N NO3 18 O NO 3 1 JOHN FOLKS DOF LAKE CITY DEEP 1 3/18/2014 5.26 16.07 2 JOHN FOLKS DOF ROSE CREEK 1 3/9/2014 5.62 9.01 3 DOT LAKE CITY 1 3/18/2014 3.51 6.73 4 BOB BRENNAN BH WOODWORKS 1 8/2/2013 4.99 2.45 5 DOT 47 1 3/7/2014 6.75 4.86 6 MT SALEM CHURCH 1 8/2/2013 7.64 8.02 7 HUNTS ALUMINUM 1 8/8/2013 4.68 5.65 8 ICHETUCKNEE MW4 1 3/9/2014 8.72 10.09 9 ICHETUCKNEE MW15 1 3/14/2015 4.01 7.86 10 ANDERSON MINE 1 3/14/2014 4.47 3.53 11 BETHEL BAPTIST 1 8/8/2013 0.55 4.17 12 ICHETUCKNEE MW1 1 3 /9/2014 3.43 5.79 13 ICHETUCKNEE MW6 1 3/7/2014 7.61 6.16 14 ICHETUCKNEE MW7 1 3/9/2014 2.89 0.26 15 ICHETUCKNEE MW9 1 3/7/2014 2.49 0.63 16 ICHETUCKNEE MW14 1 3/14/2014 1.73 4.47 1 JOHN FOLKS DOF LAKE CITY DEEP 2 8/8/2013 NM NM 2 JOHN FOLKS DOF RO SE CREEK 2 8/9/2013 4.51 5.26 3 DOT LAKE CITY 2 8/9/2013 ND ND 5 DOT 47 2 7/31/2013 6.9 4.84 8 ICHETUCKNEE MW4 2 8/1/2013 9.66 12.6 9 ICHETUCKNEE MW15 2 8/2/2013 ND ND 10 ANDERSON MINE 2 8/8/2013 4.26 3.91 12 ICHETUCKNEE MW1 2 8/1/2013 3.96 6.18 13 ICHETUCKNEE MW6 2 8/1/2013 10.02 5.06 14 ICHETUCKNEE MW7 2 8/2/2013 3.66 1.6 15 ICHETUCKNEE MW9 2 8/2/2013 3.91 0.33 16 ICHETUCKNEE MW14 2 8/2/2013 2.42 5.64
170 APPENDIX B KARST DISSOLUTION MODEL OVERVIEW Dissolution Model The following text describes the dissolution model theory which has been documented in Graham et al.  for unpublished work completed for the St. Johns Water Management District. Permission to reprint this passage has been provided by the other coauthors. Theory We assume that karst aquifers can be represented by one dimensional conduits with a circular cross section embedded in a porous limestone matrix. We consider physical and chemical processes within the conduits and the matrix. Some previous studies karst evolution studies have represented conduits as ducts having a rectangular cross section. In these studies the conduits are often referred to as fractures even when the fractures are represented by one dimensional discrete features. Other studies have considered karst disso lution in a single two dimensional fracture [ Hanna and Rajaram , 1998; Andre and Rajaram , 2005; Detwiler and Rajaram , 2007; Szymczak and Ladd , 2009, 2011; Chaudhuri et al. , 2013; Pandey et al. , 2014] . Flow a nd Reactive Solute Tra nsport Conduit flow is governed by the following mass balance equation: (B 1) where C c is capacity term for conduit flow [m], p the pressure head [m], v the velocity [m/s], A the cross sectional area of flow [m 2 ], s the spat ial coordinate in the direction parallel to the conduit [m], a sink term associated with exchange from the conduit
171 to the matrix [m 2 /s] and q c,O and q c,I are conduit sink and source terms [m 2 /s], respectively . The flow velocity v in the conduit equals Q / A . The mass balance equation for matrix flow is given by: (B 2) where C is a capacity term for matrix flow [L 1 ], p the pressure head, q the darcy flux [m/s], a sink term ass ociated with exchange from the matrix to the conduit [1/s] and q m,O and q m,I are matrix sink and source terms, r respectively [1/s]. Reactive solute transport of calcium in the conduits is governed by the following advection dispersion reaction equation: (B 3) where c is the concentration [mol/m 3 ], D is the hydrodynamic dispersion coefficient for conduit flow [m 2 /s], c I the concentration at inflow boundaries and P c a calcium production term [mol/m/s] . Reactive transport in the ma trix is governed by the following advection dispersion reaction equation: (B 4) where is the water content [ ], D the hydrodynamic dispersion tensor [m/s] and P m a calcium production term [mol/m 3 /s]. Equations ( B 3) and ( B 4) are based on the assumption that solute transport between the conduits and the matrix is solely governed by advect ion.
172 Calcite Dissolution The change in conduit radius r [m] due to a dissolution rate R [mol/m 2 /s] follows from a mass balance at the conduit wall: (B 5) where is the number of moles of calcite per unit mass of calcite [mol/k g] and the density of calcite [kg/m 3 ]. Similarly, within the porous matrix, change of porosity due to a dissolution rate R [mol/m 2 /s] is given by: (B 6) where S is the specific reaction surface [1/m] of porous limestone. Ty pically, the reaction surface per unit volume of porous material is very large and calcite undersaturated water (i.e. undersaturated with respect to calcite ) entering the porous matrix quickly becomes saturated with respect to calcite. As a result, dissolu tion of the porous matrix is effectively limited to a small region with a sharp reaction front where the calcite undersaturated water is introduced. Within the bulk of the matrix continuum the calcium concentration simply equals the saturation equilibrium concentration for calcium c eq : (B 7) The dissolution of limestone at the conduit matrix interface is governed by surface controlled and transport controlled reaction rates. The first order surface controlled dissolution rate R s [mol/m 2 /s] is given by ( B 31): (B 8)
173 where k s is the surface controlled rate coefficient [mol/m 2 /s], c s the calcium concentration at the interface [mol/m 3 ] and c eq the calcium saturation equilibrium concentration [mol/m 3 ]. T he transport controlled dissolution rate accounts for transport through the diffusion boundary layer and is given by: (B 9) where k t is the transport controlled rate coefficient and c the bulk calcium concentration within the water. Equating equation ( B 8) and equation ( B 9) gives an expression for c s which can be inserted in either one of these two equations to find the following expression for the effective first order dissolution rate R [ Szymczak and Ladd , 2009; Perne et al. , 2014] : (B 10) with: (B 11) The transport controlled rate coefficient k t is given by: (B 12) where D m is the d iffusion coefficient [m 2 /s] and the thickness of the boundary layer. The thickness of the boundary layer is defined by the Sherwood number: (B 13) The Sherwood number for laminar conduit flow is 3.66 [ Goode , 1996] . For turbulent flow the Sherwood number is derived using:
174 (B 14) where N Re and N sc are the Reynolds number and the Schmidt number, respectively. The Reynolds number is given by: (B 15) with v the velocity [m/s], w the density of water [kg/m 3 water [kg/m/s]. The Schmidt number is given by: (B 16) It has been observed that as calcium concentrations approach saturation the reaction rate decreases due to impurities within the limestone which inhibit dissolution [ Svensson an d Dreybrodt , 1992; Cornaton and Perrochet , 2006] . This phenomenon is known as the kinetic trigger effect [ White , 1977] and has been modeled by switching dissolution from first order to higher order kinetics when the calcium concentration exceeds a certain value c * . The higher order effective rate is typically given by [ Lichtner , 1988] : (B 17 ) with k n defined as: (B 18) such that R 1 = R n at c = c * . A general expression for the reaction rate be written as: (B 19)
175 Figure B 1 illustrates the effect of the kinetic trigger on reaction rate R . The decrease in dissolution rates close to saturation allows calcite undersaturated water to penetrate further into the aquifer than would otherwise be possible. Figure B 1. Effect of the kinetic trigger on reaction rate as calcium concentration a pproaches saturation (c*=1.6 mol/m 2 )). Solid line corresponds to equations (B 10) and (B 17). Dashed line corres ponds to equations (B 26) and (B 27). Model Design The quasi steady state approximation Combining equation ( B 3) and ( B 5) and using P c rR results in the following reactive transport equation: ( B 20) Research has shown [ Lichtner , 1988; Hanna and Rajaram , 1998] that because the density of the limestone rock is much larger than the maximum calcium concentration, the rate of change in conduit radius is much slower than the rate of
176 change in concentration a nd the rate of change in the flow field. Thus the flow and stationary state equations: ( B 21) Within the mat rix, we assume that the porosity remains constant and that the concentration of calcium equals the equilibrium concentration. Therefore the flow and reactive transport equations in the matrix are simplified: ( B 22) Equations ( B 21) and ( B 22) allow simulation of conduit generation processes through a sequence of steady states [ Lichtner , 1988; Hanna and Rajaram , 1998] . To begin, steady state flow and concentrations fields and corresponding dissolution rates are computed based on initial conduit diameters. The dissolution rate, in turn, determines the rate of conduit r adius enlargement. The quasi steady rate of conduit diameters. The process is then repeated using the modified conduit diameters. A sequence of these dissolution time step s can be applied to simulate the dissolution process over the desired geologic timescale. Equations ( B 21) and ( B 22) constitute a speleogenesis model that solves advective dispersive reactive transport within the conduits. This is different from many
177 ot her speleogenesis models in which advective reactive transport is solved within the conduits and fractures [ Siemers and Dreybrodt , 1998; Kaufmann and Braun , 2000; Perne et al. , 2014] . The advection reaction equation is independent of downstream conditions and may be solved from upstream to downstream. Equations are solved for each conduit cell separately with complete mix ing assumed at conduit junctions. The strength of our scheme lies in the fact that it can be easily implemented in any model code capable of handling ad vective dispersive transport. Numerical Solution o f Flow The numerical solution of flow in the conduits and the porous limestone matrix follows the approach described by [ de Rooij et al. , 2013] . This solution is based on a discrete continuum approach and a finite difference scheme. The coupling of conduit matrix flow is governed by a Peaceman well index which depends on the conduit radius. Thus, after each dissolution timestep the Peaceman well indices are updated. To permit efficient steady state flow computations for large regional domains, instead of using Richards equation to simulate variably saturated flow in the porous matrix, an option was added to solve for flow using the 3 D saturated flow equation. The height of the model domain is adjusted in accordance with the change in height of the water table at each time step, us ing an approach similar to that used in MODFLOW. Contrary to MODFLOW, however, net recharge is applied to the topmost model cells even if the water table drops below the cell. The vertical hydraulic conductivity is assumed to remain constant, at its satura ted value, to transmit recharge to the water table. Within the conduits the pipe flow equation proposed by Swamee and Swamee [ Swamee and Swamee , 2007] is implemented, assuming the conduits always remain full. This equation provides for a smooth transition between laminar and turbulent flow.
178 For laminar flow the equation approximates the Poiseuille equation. For turbulent flow, the equation approximates the Darcy Weisbach equatio n. The Swamee and Swamee equation allows conduit flow to automatically switch from laminar to turbulent conditions during c onduit evolution (Figure B 2). Figure B 2: Swamee and Swamee pipe flow equation showing smooth transition between laminar and turbu lent flow as hydraulic gradient in pipe increases (solid line). Poiseuille equation for laminar flow (dotted line) and Darcy Weisbach equation for turbulent flow (dashed line) are also shown for comparison. Numerical Solution o f Reactive Transport React ive transport is also simulated using finite differences. Accurate numerical solution of advection dispersion equations is subject to criteria for spatial as well as temporal discretization, which are typically given in terms of Courant and Peclet numbers. To avoid small space and time discretization that would result in extremely long computation times for large regional models, we use an upwind scheme that is unconditionally stable regardless of discretization. The drawback of upwinding is that it introdu ces numerical dispersion.
179 In the numerical solution of reactive transport the reaction term must be handled carefully to avoid numerical instability and time stepping restrictions. Reaction rates defined by equations ( B 10) and ( B 17) can result in numeri cal instability, with reaction rates jumping between first and higher order during non linear iterations. This behavior is likely the due to the fact that the derivative of the reaction rate is highly discontinuous at c = c * (Figure B 1). Therefore the fol lowing expression for the higher order effective rate, which has a continuous derivative at c * , was adopted [ Svensson and Dreybrodt , 1992; Andr e and Rajaram , 2005; Detwiler and Rajaram , 2007] : ( B 26) with: ( B 27) Computations with this relationship were found to be more efficient. Figure B 1 illustrates that the adapted expression com pares reasonably well to the original expression.
180 LIST OF REFERENCES Aeschbach Hertig, W., F. Peeters, U. Beyerle, and R. Kipfer (1999a), Interpretation of dissolved atmospheric noble gases in natural waters , Water Resour. Res. , 35 (9), 2779 2792, doi:10.1029/1999WR900130. Aeschbach Hertig, W., F. Peeters, U. Beyerle, and R. Kipfer (1999b), Interpretation of dissolved atmospheric noble gases in natural waters, Water Resour. Res. , 35 (9), 2779 2792. Aller, L., T Standardized Method for Evaluating Ground Water Pollution Potential Using Hydrogeologic Settings, NWWA/Epa 600/2 87 035 , 455. Amberger, A., and H. L. Schmidt (1987), The natural isotope content of nitrate as an indicator of its origin, Geochim. Cosmochim. Acta , 51 (10), 2699 2705. Andersson, J., A. M. Shapiro, and J. Bear (1984), A stochastic model of a fractured rock conditioned by measured information, Water Resour. Res. , 20 (1), 79 88. Andre, B. J., and H. Rajaram (2005), Dissolution of limestone fractures by cooling waters: Early development of hypogene karst systems, Water Resour. Res. , 41 (1). Atkinson, T. C. (1977), Diffuse flow and conduit flow in limestone terrain in the Mendip Hill s, Somerset (Great Britain), J. Hydrol. , 35 (1 2), 93 110, doi:10.1016/0022 1694(77)90079 8. Bae, H. S., F. E. Dierberg, and A. Ogram (2014), Syntrophs dominate sequences associated with the mercury methylation related gene hgcA in the water conservation ar eas of the Florida Everglades, Appl. Environ. Microbiol. , 80 (20), 6517 6526. Baedke, S. , and N. C. Krothe (2001), Derivation of effective hydraulic parameters of a karst aquifer from discharge hydrograph analysis, Water Resour. Res. , 37 (1), 13 19. Bailly C omte, V., J. B. Martin, and E. J. Screaton (2011), Time variant cross correlation to assess residence time of water and implication for hydraulics of a sink rise karst system, Water Resour. Res. , 47 (5), 1 16, doi:10.1029/2010WR009613. Batu, V. (1998), Aqui fer hydraulics: a comprehensive guide to hydrogeologic data analysis , John Wiley & Sons. Bauer, S., R. Liedl, and M. Sauter (2005), Modeling the influence of epikarst evolution on karst aquifer genesis: A time variant recharge boundary condition for joint karst epikarst development, Water Resour. Res. , 41 (9), 1 12, doi:10.1029/2004WR003321.
181 Beck, B. F. (1986), A generalized genetic framework for the development of sinkholes and karst in Florida, USA, Environ. Geol. Water Sci. , 8 (1 2), 5 18. Bengtsson, G., a nd H. Annadotter (1989), Nitrate reduction in a groundwater microcosm determined by N gas chromatography mass spectrometry., Appl. Environ. Microbiol. , 55 (11), 2861 70. Berndt, M. P., E. T. Oaksford, G. L. Mahon, and W. Schmidt (1998), Groundwater, Water R esour. Atlas Florida Florida State Univ., Inst. Public Aff. Biesterfeld, S., G. Farmer, L. Figueroa, D. Parker, and P. Russell (2003), Quantification of denitrification potential in carbonaceous trickling filters, Water Res. , 37 (16), 4011 4017, doi:10.1016 /S0043 1354(03)00302 6. Birk, S. (2003), Hydraulic boundary conditions as a controlling factor in karst genesis: A numerical modeling study on artesian conduit development in gypsum, Water Resour. Res. , 39 (1), 1 14, doi:10.1029/2002WR001308. Birk, S., R. L iedl, and M. Sauter (2004), Identification of localised recharge and conduit flow by combined analysis of hydraulic and physico chemical spring responses (Urenbrunnen, SW Germany), J. Hydrol. , 286 (1), 179 193. Blakey, N. C., and P. A. Towler (1988), The ef fect of unsaturated/saturated zone property upon the hydrogeochemical and microbiological processes involved in the migration and attenuation of landfill leachate components, Water Sci. Technol. , 20 (3), 119 128. B Ã¶ hlke, J. K., and J. M. Denver (1995), Comb ined use of groundwater dating, chemical, and isotopic analyses to resolve the history and fate of nitrate contamination in two agricultural watersheds, Atlantic coastal plain, Maryland, Water Resour. Res. , 31 (9), 2319 2339, doi:10.1029/95WR01584. BÃ¶hlke, J. K. (2002), Denitrification in the recharge area and discharge area of a transient agricultural nitrate plume in a glacial outwash sand aquifer, Minnesota, Water Resour. Res. , 38 (7), doi:10.1029/2001WR000663. Bonacci, O. (1987), Karst Hydrology , Herdelbe rg, Germany. Bonacci, O. (1993), Karst springs hydrographs as indicators of karst aquifers, Hydrol. Sci. J. , 38 (August), 51 62, doi:10.1080/02626669309492639. Bonn, M. A., and B. M. R. Group (2004), Visitor profiles, economic impacts and recreational aesth etic values associated with eight priority Florida springs located in the St. Johns River Water Management District , St. Johns River Water Management District.
182 Borghi, A., P. Renard, and S. Jenni (2010), How to Model Realistic 3D Karst Reservoirs Using a Pseudo Genetic Methodology Example of Two Case Studies, in Advances in Research in Karst Media , pp. 251 255, Springer. Borghi, A., P. Renard, and S. Jenni (2012), A pseudo genetic stochastic model to generate karstic networks, J. Hydrol. , 414 415 , 516 529, doi:10.1016/j.jhydrol.2011.11.032. Borghi, A., P. Renard, and F. Cornaton (2016), Can one identify karst conduit networks geometry and properties from hydraulic and tracer test data?, Adv. Water Resour. , 90 , 99 115. BÃ¶ttcher, J., O. Strebel, S. Voerkelius , and H. L. Schmidt (1990), Using isotope fractionation of nitrate nitrogen and nitrate oxygen for evaluation of microbial denitrification in a sandy aquifer, J. Hydrol. , 114 (3), 413 424. Bratbak, G. (1985), Bacterial biovolume and biomass estimations, App l. Environ. Microbiol. , 49 (6), 1488 1493. Bratbak, G., and I. A. N. Dundas (1984), Bacterial Dry Matter Content and Biomass Estimations, , 48 (4), 755 757. Brown, M., K. Reiss, M. Cohen, and J. Evans (2008), Summary and synthesis of the available literature on the effects of nutrients on spring organisms and systems, , (April 2008). Brunet, R. C., and L. J. Garcia Gil (1996), Sulfide induced dissimilatory nitrate reduction to ammonia in anaerobic freshwater sediments, FEMS Microbiol. Eco l. , 21 (2), 131 138, doi:10.1016/0168 6496(96)00051 7. Burgin, A. J., and S. K. Hamilton (2007), Have we overemphasized in aquatic removal Front. Ecol. Environ. , 5 (2), 89 96, doi:10.1890/1540 9295(2007)5[89:HWOTRO]2.0.CO;2. Bush, P. W., and R . H. Johnson (1988), Ground water hydraulics, regional flow, and ground water development of the Floridan aquifer system in Florida and in parts of Georgia, South Carolina, Buss, S. R., a. W. Herbert, P. Morgan, S. F. Thornton, and J. W. N. Smith (2004), A review of ammonium attenuation in soil and groundwater, Q. J. Eng. Geol. Hydrogeol. , 37 (4), 347 359, doi:10.1144/1470 9236/04 005. Cacas, M. C., E. Ledoux, G. de Marsily, A. Barbreau, P. Calmels, B. Gaillard, and R. Margritta (1990), Modeling fracture fl ow with a stochastic discrete fracture network: Calibration and validation: 2. The transport model, Water Resour. Res. , 26 (3), 491 500.
183 Campolongo, F., J. Cariboni, and A. Saltelli (2007), An effective screening design for sensitivity analysis of large mo dels, Environ. Model. Softw. , 22 (10), 1509 1518, doi:10.1016/j.envsoft.2006.10.004. Casciotti, K. L., D. M. Sigman, M. G. Hastings, J. K. Bohlke, and a Hilkert (2002), Measurement of the oxygen isotopic composition of nitrate seawater and freshwater using the dentirifier method., Anal. Chem. , 74 (19), 4905 4912. Ceryak, R., M. S. Knapp, and T. Burnson (1983), The geology and water resources of the upper Suwannee River basin, Florida , Published for the Bureau of Geology, Division of Resource Management, Flor ida Department of Natural Resources. Champion, K. M., and S. B. Upchurch (2003a), Delineation of Spring Water Source Areas in the Ichetucknee Springshed, Champion, K. M., and S. B. Upchurch (2003b), Delineation of Spring Water Source Areas in the Ichetuckn ee Springshed . Chang, Y., C. Ho, C. Chang, and S. Tseng (2006), Denitrificati on Under High Dissolved Oxygen b y a Membrane Attached Biofilm Reactor , , 29 (4), 741 745. Chang, Y., J. Wu, and L. Liu (2015), Effects of the conduit network on the spring hydrogra ph of the karst aquifer, J. Hydrol. , 527 , 517 530, doi:10.1016/j.jhydrol.2015.05.006. Chaudhuri, A., H. Rajaram, and H. Viswanathan (2013), Early stage hypogene karstification in a mountain hydrologic system: A coupled thermohydrochemical model incorporating buoyant convection, Water Resour. Res. , 49 (9), 5880 5899. Cohen, M. J., S. Lamsal, and L. V. Kohrnak (2007), Sources, Transport and Transforma tions of Nitrate N in the Florida Environment, , 1 128. Cohen, M. J., J. B. Heffernan, A. Albertin, and J. B. Martin (2012), Inference of riverine nitrogen processing from longitudinal and diel variation in dual nitrate isotopes, J. Geophys. Res. Biogeosci ences , 117 (1), doi:10.1029/2011JG001715. Confalonieri, R., G. Bellocchi, S. Bregaglio, M. Donatelli, and M. Acutis (2010), Comparison of sensitivity analysis techniques: A case study with the rice model WARM, Ecol. Modell. , 221 (16), 1897 1906, doi:10.1016/ j.ecolmodel.2010.04.021. Coplen, T. B., H. Qi, K. Revesz, K. Casciotti, and J. E. Hannon (2007), Determination of Cornaton, F. (1999), Utilisation des modeles continu discret et a double continuum pour
184 Cornaton, F., and P. Perrochet (2006), Groundwater age, life expectancy and transit time distributions in advective dispersive systems: 1. Generalized reservoir theory, Adv. Wa ter Resour. , 29 (9), 1267 1291. Covington, M. D., C. M. Wicks, and M. O. Saar (2009), A dimensionless number describing the effects of recharge and geometry on discharge from simple karstic aquifers, Water Resour. Res. , 45 (11), 1 16, doi:10.1029/2009WR00800 4. Covino, T. P., B. L. McGlynn, and R. a. McNamara (2010), Tracer Additions for Spiraling Curve Characterization (TASCC): Quantifying stream nutrient uptake kinetics from ambient to saturation, Limnol. Oceanogr. Methods , 8 , 484 498, doi:10.4319/lom.2010.8 .484. Coxon, C. (1999), Agriculturally induced impacts, Int. Contrib. to Hydrogeol. , 20 , 37 80. Dafny, E., A. Burg, and H. Gvirtzman (2010), Effects of Karst and geological structure on groundwater flow: The case of Yarqon Taninim Aquifer, Israel, J. Hydro l. , 389 (3 4), 260 275, doi:10.1016/j.jhydrol.2010.05.038. Davidson, E. A., and S. Seitzinger (2006), The enigma of progress in denitrification research, Ecol. Appl. , 16 (6), 2057 2063, doi:10.1890/1051 0761(2006)016[2057:TEOPID]2.0.CO;2. Davidson, E. A., J. Chorover, and D. B. Dail (2003), A mechanism of abiotic immobilization of nitrate in forest ecosystems: the ferrous wheel hypothesis, Glob. Chang. Biol. , 9 (2), 228 236. Denizman, C. (1998), Evolution Of Karst i n t he Lower Suwannee River Basin, Florida , Un iversity of Florida. de Rooij, R., P. Perrochet, and W. Graham (2013), From rainfall to spring discharge: Coupling conduit flow, subsurface matrix flow and surface flow in karst systems using a discrete continuum model, Adv. Water Resour. , 61 , 29 41, doi:1 0.1016/j.advwatres.2013.08.009. Detwiler, R. L., and H. Rajaram (2007), Predicting dissolution patterns in variable aperture fractures: Evaluation of an enhanced depth averaged computational model, Water Resour. Res. , 43 (4). Dewandel, B., P. Lachassagne, M. Bakalowicz, P. Weng, and A. Al Malki (2003), Evaluation of aquifer thickness by analysing recession hydrographs. Application to the Oman ophiolite hard rock aquifer, J. Hydrol. , 274 (1 4), 248 269, doi:10.1016/S0022 1694(02)00418 3. Diersch, H. J. G. (2005), FEFLOW finite element subsurface flow and transport simulation system, Inst. Water Resour. Plan. Syst. Res., Berlin .
185 Doummar, J., M. Sauter, and T. Geyer (2012), Simulation of flow processes in a large scale karst system with an integrated catchment model (Mike She) Identification of relevant parameters influencing spring discharge, J. Hydrol. , 426 427 , 112 123, doi:10.1016/j.jhydrol.2012.01.021. Dreybrodt, W. (1 990), The role of dissolution kinetics in the development of karst aquifers in limestone: a model simulation of karst evolution, J. Geol. , 639 655. Dreybrodt, W. (1996), Principles of early development of karst conduits under natural and man made conditio ns revealed by mathematical analysis of numerical models, Water Resour. Res. , 32 (9), 2923 2935, doi:10.1029/96WR01332. Eisenlohr, L. (1996), Eisenlohr, L., 1996. Variabilite des reponses naturelles des aquiferes karstiques., Universite de Neuchatel. Eisenlohr, L., L. KirÃ¡ly, M. Bouzelboudjen, and Y. Rossier (1997a), Numerical simulation as a tool for checking the interpretation of karst spring hydrographs, J. Hydrol. , 193 (1 4), 306 315, doi:10.1016/S0022 1694(96)03140 X. Eisenlohr, L., M. Bouzelboudje n, L. KirÃ¡ly, and Yvan Rossier (1997b), Numerical versus statistical modelling of natural response of a karst hydrogeological system, J. Hydrol. , 202 (1 4), 244 262, doi:10.1016/S0022 1694(97)00069 3. Ewers, R. O. (2006), Karst aquifers and the role of assu mptions and authority in science, Geol. Soc. Am. Spec. Pap. , 404 (19), 235 242, doi:10.1130/2006.2404(19). Fagerbakke, K. M., M. Heldal, and S. Norland (1996), Content of carbon, nitrogen, oxygen, sulfur and phosphorus in native aquatic and cultured bacteri a, Aquat. Microb. Ecol. , 10 (1), 15 27, doi:10.3354/ame010015. Faul, F., E. Erdfeld, A. G. Lang, and A. Buchner (2007), A flexible statistical power analysis program for the social, behavioral, and biomedical sciences., Behav. Res. Methods , 39 (2), 175 191. Faulkner, G. L. (1970), Geohydrology of the Cross Florida Barge Canal Area: With Special Reference to the Ocala Vicinity , Diane Publishing. FDEP (1990), DRASTIC Index for Floridan Aquifer, Fiorillo, F. (2014), The Recession of Spring Hydrographs, Focused o n Karst Aquifers, Water Resour. Manag. , 28 (7), 1781 1805, doi:10.1007/s11269 014 0597 z. Florea, L. J., and H. L. Vacher (2007a), Eogenetic karst hydrology: Insights from the 2004 hurricanes, peninsular Florida, Ground Water , 45 (4), 439 446, doi:10.1111/j. 1745 6584.2007.00309.x.
186 Florea, L. J., and H. L. Vacher (2007b), Eogenetic karst hydrology: Insights from the 2004 hurricanes, peninsular Florida, Ground Water , 45 (4), 439 446, doi:10.1111/j.1745 6584.2007.00309.x. Florea, L. J., H. L. Vacher, B. Donahue, and D. Naar (2007), Quaternary cave levels in peninsular Florida, Quat. Sci. Rev. , 26 (9 10), 1344 1361, doi:10.1016/j.quascirev.2007.02.011. Florida Department of Environmental Protection (2012), Critical Lands Report: Ichetucknee Trace. Florida Geologica l Survey (FGS) (2003), Sinks and Swallets Spatial Data, Florida Geological Survey (FGS) (2004), Closed Topographic Depressions, Foglia, L., M. C. Hill, S. W. Mehl, and P. Burlando (2009), Sensitivity analysis, calibration, and testing of a distributed hydr ological model using error based weighting and one objective function, Water Resour. Res. , 45 (6), 1 18, doi:10.1029/2008WR007255. Ford, D., and P. D. Williams (2013), Karst hydrogeology and geomorphology , John Wiley & Sons. Fournillon, F., S. Viseur, B. Ar fib, and J. Borgomano (2010), Insights of 3D geological modelling in distributed hydrogeological models of karstic carbonate aquifers, in Advances in Research in Karst Media , pp. 257 262, Springer. Gallegos, J. J., B. X. Hu, and H. Davis (2013), Simulating flow in karst aquifers at laboratory and sub regional scales using MODFLOW CFP, Hydrogeol. J. , 21 (8), 1749 1760, doi:10.1007/s10040 013 1046 4. Gardner, P., and D. K. Solomon (2009), An advanced passive diffusion sampler for the determination of dissolved gas concentrations, Water Resour. Res. , 45 (6). Geyer, T., S. Birk, R. Liedl, and M. Sauter (2008a), Quantification of temporal distribution of recharge in karst systems from spring hydrographs, J. Hydrol. , 348 (3), 452 463. Geyer, T., S. Birk, R. Liedl, an d M. Sauter (2008b), Quantification of temporal distribution of recharge in karst systems from spring hydrographs, J. Hydrol. , 348 (3 4), 452 463, doi:10.1016/j.jhydrol.2007.10.015. Giblin, a E., C. R. Tobias, B. Song, N. Weston, G. T. Banta, and V. H. Riv era Monroy (2013), The importance of dissimilatory nitrate reduction to ammonium (DNRA) in the nitrogen cycle of coastal ecosystems, Oceanography , 26 (3), 124 131, doi:10.5670/oceanog.2013.54.
187 Goode, D. J. (1996), Direct simulation of groundwater age, Wate r Resour. Res. , 32 (2), 289 296. Gray, C. J., and A. S. Engel (2013), Microbial diversity and impact on carbonate geochemistry across a changing geochemical gradient in a karst aquifer., ISME J. , 7 (2), 325 37, doi:10.1038/ismej.2012.105. Green, C. T., L. J. Puckett, J. K. BÃ¶hlke, B. A. Bekins, S. P. Phillips, L. J. Kauffman, J. M. Denver, and H. M. Johnson (2008), Limited occurrence of denitrification in four shallow aquifers in agricultural areas of the United States, J. Environ. Qual. , 37 (3), 994 1009. Gri mm, N. B., S. E. Gergel, W. H. McDowell, E. W. Boyer, C. L. Dent, P. Groffman, S. C. Hart, J. Harvey, C. Johnston, and E. Mayorga (2003), Merging aquatic and terrestrial perspectives of nutrient biogeochemistry, Oecologia , 137 (4), 485 501. Groffman, P. M. (2012), Terrestrial denitrification: challenges and opportunities, Ecol. Process. , 1 (1), 11, doi:10.1186/2192 1709 1 11. Groffman, P. M., E. A. Davidson, and S. Seitzinger (2009), New approaches to modeling denitrification, Biogeochemistry , 1 5, doi:10.100 7/s10533 009 9285 0. Groves, C. G., and A. D. Howard (1994), Early development of karst systems: 1. Preferential flow path enlargement under laminar flow, Water Resour. Res. , 30 (10), 2837 2846, doi:10.1029/94WR01303. Guerrero, M. G. (1985), Assimilatory ni trate reduction., Guerrero, M. G., J. M. Vega, and M. Losada (1981), The assimilatory nitrate reducing system and its regulation, Annu. Rev. Plant Physiol. , 32 (1), 169 204. Gulley, J., J. Martin, P. Spellman, P. Moore, and E. Screaton (2013a), Dissolution in a variably confined carbonate platform: Effects of allogenic runoff, hydraulic damming of groundwater inputs, and surface groundwater exchange at the basin scale, Earth Surf. Process. Landforms , 38 (14), 1700 1713, doi:10.1002/esp.3411. Gulley, J. D., J. B. Martin, P. J. Moore, and J. Murphy (2013b), Formation of phreatic caves in an eogenetic karst aquifer by CO2 enrichment at lower water tables and subsequent flooding by sea level rise, Earth Surf. Process. Landforms , 38 (11), 1210 1224, doi:10.1002/esp. 3358. Hamme, R. C., and S. R. Emerson (2004), The solubility of neon, nitrogen and argon in distilled water and seawater, Deep. Res. Part I Oceanogr. Res. Pap. , 51 (11), 1517 1528, doi:10.1016/j.dsr.2004.06.009. Hancock, P. J., A. J. Boulton, and W. F. Hump hreys (2005), Aquifers and hyporheic zones: Towards an ecological understanding of groundwater, Hydrogeol. J. , 13 (1), 98 111, doi:10.1007/s10040 004 0421 6.
188 Hanna, R. B., and H. Rajaram (1998), Influence of aperture variability on dissolutional growth of f issures in Karst Formations, Water Resour. Res. , 34 (11), 2843, doi:10.1029/98WR01528. Hanshaw, B. B., and W. Back (1979), Major geochemical processes in the evolution of carbonate aquifer systems, J. Hydrol. , 43 (1), 287 312. Heath, R. C. (1983), Basic Grou nd Water Hydrology . Heaton, T. H. E. (1986), Isotopic studies of nitrogen pollution in the hydrosphere and atmosphere: a review, Chem. Geol. Isot. Geosci. Sect. , 59 , 87 102. Heffernan, J. B., A. R. Albertin, M. L. Fork, B. G. Katz, and M. J. Cohen (2012), Denitrification and inference of nitrogen sources in the karstic Floridan Aquifer, Biogeosciences , 9 (5), 1671 1690, doi:10.5194/bg 9 1671 2012. Henry, S., E. Baudoin, J. C. LÃ³pez GutiÃ©rrez, F. Martin Laurent, A. Brauman, and L. Philippot (2004), Quantifica tion of denitrifying bacteria in soils by nirK gene targeted real time PCR, J. Microbiol. Methods , 59 (3), 327 335, doi:10.1016/j.mimet.2004.07.002. Herman, J. D., J. B. Kollat, P. M. Reed, and T. Wagener (2013), Technical Note: Method of Morris effectively reduces the computational demands of global sensitivity analysis for distributed watershed models, Hydrol. Earth Syst. Sci. , 17 (7), 2893 2903, doi:10.5194/hess 17 2893 2013. Hill, A. R. (1996), Nitrate Removal in Stream Riparian Zones, J. Environ. Qual. , 25 (4), 743, doi:10.2134/jeq1996.00472425002500040014x. Hollibaugh, J. T., S. M. Gifford, M. A. Moran, M. J. Ross, S. Sharma, and B. B. Tolar (2014), Seasonal variation in the metatranscriptomes of a Thaumarchaeota population from SE USA coastal waters., IS ME J. , 8 (3), 685 98, doi:10.1038/ismej.2013.171. Holm, P. E., P. H. Nielsen, H. J. Albrechtsen, and T. H. Christensen (1992), Importance of unattached bacteria and bacteria attached to sediment in determining potentials for degradation of xenobiotic organi c contaminants in an aerobic aquifer, Appl. Environ. Microbiol. , 58 (9), 3020 3026. Hornsby, D., and R. Mattson (1998), Surface water quality and biological monitoring network: Live Oak, Fla., Suwannee River Water Manag. Dist. Annu. Rep. WR 98 02 . Hornsby, H. D., and R. Ceryak (1998), Springs of the Suwannee River basin in Florida , Department of Water Resources, Suwannee River Water Management District. Houston, T. B. (1965), Soil Survey, Suwannee County, Florida , US Department of Agriculture, Soil Conservat ion Service.
189 Howard, A. D., and C. G. Groves (1995), Early Development of Karst Systems: 2. Turbulent Flow, Water Resour. Res. , 31 (1), 19, doi:10.1029/94WR01964. Hsu, C. H., S. T. Han, Y. H. Kao, and C. W. Liu (2010), Redox characteristics and zonation of arsenic affected multi layers aquifers in the Choushui River alluvial fan, Taiwan, J. Hydrol. , 391 (3), 351 366. Humbert, S., S. Tarnawski, N. Fromin, M. P. Mallet, M. Aragno, and J. Zopfi (2010), Molecular detection of anammox bacteria in terrestrial ecosy stems: distribution and diversity, Isme J. , 4 (3), 450 454, doi:10.1038/ismej.2009.125. Hunn, J. D., and L. J. Slack (1983), Water Resources of the Santa Fe River Basin, Florida , US Geological Survey,. Hydrogeologic (2013), Northern District Groundwater Flo w Model . Iooss, B., and P. LemaÃ®tre (2014), A Review o n Global Sensitivity Analysis Methods , arXiv:1404.2405 [math, stat] . Single ination of Microbial Activities, Ground Water , 35 (4), 619 631, doi:10.1111/j.1745 6584.1997.tb00127.x. James, A. L., and C. M. Oldenburg (1997), Linear and Monte Carlo uncertainty analysis for subsurface contaminant transport simulation, Water Resour. Res. , 33 (11), 2495, doi:10.1029/97WR01925. Jaquet, O., and P. Y. Jeannin (1994), Modelling the karstic medium: a geostatistical approach, in Geostatistical Simulations , pp. 185 195, Springer. Jaquet, O., P. Siegel, G. Klubertanz, and H. Benabderrhamane (2004), Stochastic discrete model of karstic networks, Adv. Water Resour. , 27 (7), 751 760, doi:10.1016/j.advwatres.2004.03.007. Jones, G. W., S. B. Upchurch, K. M. Champion, and E. C. DeHaven (1996), Origin of nitrate in ground water discharging from Rainbow Spri ngs, Marion County, Florida, Kandeler, E., K. Deiglmayr, D. Tscherko, D. Bru, and L. Philippot (2006), Abundance of narG, nirS, nirK, and nosZ genes of denitrifying bacteria during primary successions of a glacier foreland, Appl. Environ. Microbiol. , 72 (9) , 5957 5962, doi:10.1128/AEM.00439 06. Kartal, B. et al. (2011), Molecular mechanism of anaerobic ammonium oxidation., Nature , 479 (7371), 127 30, doi:10.1038/nature10453. Katsanou, K., N. Lambrakis, G. Tayfur, and A. Baba (2015), Describing the Karst Evolu tion by the Exploitation of Hydrologic Time Series Data, Water Resour. Manag. , 29 (9), 3131 3147, doi:10.1007/s11269 015 0987 x.
190 Katz, B. (1999), Sources and chronology of nitrate contamination in spring waters, Suwannee River Basin, Florida, U.S. Geol. Sur v. Water Resour. Investig. Rep. , #99 4252 . Katz, B., R. Copeland, T. Greenhalgh, R. Ceryak, and W. Zwanka (2005), Using multiple chemical indicators to assess sources of nitrate and age of groundwater in a karstic spring basin, Environ. Eng. Geosci. , 11 (4) , 333 346, doi:10.2113/11.4.333. Katz, B. G. (1992), Hydrochemistry of the upper Floridan aquifer, Florida , US Geological Survey. Katz, B. G. (2004), Sources of nitrate contamination and age of water in large karstic springs of Florida, Environ. Geol. , 46 ( 6 7), 689 706, doi:10.1007/s00254 004 1061 9. Katz, B. G., and D. W. Griffin (2008), Using chemical and microbiological indicators to track the possible movement of contaminants from the land application of reclaimed municipal wastewater in a karstic sprin gs basin, Environ. Geol , 55 , 801 821. Katz, B. G., J. K. BÃ¶hlke, and H. D. Hornsby (2001), Timescales for nitrate contamination of spring waters, northern Florida, USA, Chem. Geol. , 179 (1 4), 167 186, doi:10.1016/S0009 2541(01)00321 7. Katz, B. G., A. R. C helette, and T. R. Pratt (2004), Use of chemical and isotopic tracers to assess nitrate contamination and ground water age, Woodville Karst Plain, USA, J. Hydrol. , 289 (1 4), 36 61, doi:10.1016/j.jhydrol.2003.11.001. Katz, B. G., A. A. Sepulveda, and R. J. Verdi (2009), Estimating nitrogen loading to ground water and assessing vulnerability to nitrate contamination in a large karstic springs Basin, Florida, J. Am. Water Resour. Assoc. , 45 (3), 607 627, doi:10.1111/j.1752 1688.2009.00309.x. Kaufmann, G. (2003) , A model comparison of karst aquifer evolution for different matrix flow formulations, J. Hydrol. , 283 (1 4), 281 289, doi:10.1016/S0022 1694(03)00270 1. Kaufmann, G. (2009), Modelling karst geomorphology on different time scales, Geomorphology , 106 (1), 62 77. Kaufmann, G., and J. Braun (2000), Karst aquifer evolution in fractured, porous rocks, Water Resour. Res. , 36 (6), 1381 1391, doi:10.1029/1999WR900356. Kaufmann, G., D. Romanov, and T. Hiller (2010), Modeling three dimensional karst aquifer evolution u sing different matrix flow contributions, J. Hydrol. , 388 (3 4), 241 250, doi:10.1016/j.jhydrol.2010.05.001.
191 Kelso, B. H. L., R. V. Smith, R. J. Laughlin, and S. D. Lennox (1997), Dissimilatory nitrate reduction in anaerobic sediments leading to river nitr ite accumulation, Appl. Environ. Microbiol. , 63 (12), 4679 4685. Kendall, C., and R. Aravena (2000), Nitrate isotopes in groundwater systems, in Environmental tracers in subsurface hydrology , pp. 261 297, Springer. Kendall, C., and J. J. McDonnell (2012), I sotope tracers in catchment hydrology , Elsevier. Khare, Y., and R. MuÃ±oz Effects Screening Method of Morris using Sampling Uniformity ( SU ) MatLab code manual, , 0570 (352). Kim, J. H., C. Y. Ha, S. W. Oa, J. W. Lee, S. H. Park, S. Y. Kwon, S. Kim, and Y. Kim (2011), Assessing the activity and diversity of fumarate fed denitrifying bacteria by performing field single well push pull tests, J. Environ. Sci. Heal. Part A , 46 (1), 33 41. Kim, Y., J. H. Kim , B. H. Son, and S. W. Oa (2005), A single well push pull test method for in situ determination of denitrification rates in a nitrate contaminated groundwater aquifer, Water Sci. Technol. , 52 (8), 77 86. Kiraly, L. (1985), FEM 301 a three dimensional model for groundwater flow simulation , Nationale Genossenschaft fuer die Lagerung Radioaktiver Abfaelle (NAGRA). Kiraly, L. (2002), Karstification and groundwater flow, Evol. Karst Prekarst to Cessat. Postojna Ljubljana, Zaloz. ZRC , 155 190. KirÃ¡ly, L., and G. M simulÃ© par modÃ¨les mathÃ©matiques, , 1 , 37 60. Klappenbach, J. a, P. R. Saxman, J. R. Cole, and T. M. Schmidt (2001), rrndb: the Ribosomal RNA Operon Copy Number D atabase., Nucleic Acids Res. , 29 (1), 181 184, doi:10.1093/nar/29.1.181. Kolle, W., O. Strebel, and J. Bottcher (1985), Formation of sulfate by a microbial denitrification in a reducing aquifer, Water Supply , 3 (1), 35 40. Korom, S. F. (1992), Natural denitr ification in the saturated zone: A review, Water Resour. Res. , 28 (6), 1657 1668, doi:10.1029/92WR00252. KovÃ¡cs, A. (2003), Geometry and hydraulic parameters of karst aquifers: a hydrodynamic modeling approach, , 1 131. KovÃ¡c s, A., and P. Perrochet (2008), A quantitative approach to spring hydrograph decomposition, J. Hydrol. , 352 (1 2), 16 29, doi:10.1016/j.jhydrol.2007.12.009.
192 KovÃ¡cs, A., P. Perrochet, L. KirÃ¡ly, and P. Y. Jeannin (2005), A quantitative method for the charact erisation of karst aquifers based on spring hydrograph analysis, J. Hydrol. , 303 (1 4), 152 164, doi:10.1016/j.jhydrol.2004.08.023. Kuenen, J. G., B. J. Nsen, and N. P. Revsnech (1986), O xygen microprofiles of trickling filter biofilms , 20 (12), 1589 1598. Langston, A. L., E. J. Screaton, J. B. Martin, and V. Bailly Comte (2012a), Interactions of diffuse and focused allogenic recharge in an eogenetic karst aquifer (Florida, USA), Hydrogeol. J. , 20 (4), 767 781, doi:10.1007/s10040 012 0845 3. Langston, A . L., E. J. Screaton, J. B. Martin, and V. Bailly Comte (2012b), Interactions of diffuse and focused allogenic recharge in an eogenetic karst aquifer (Florida, USA), Hydrogeol. J. , 20 (4), 767 781. Lauritzen, S. E., N. Odling, and J. Petersen (1992), Modeli ng the evolution of channel networks in carbonate rocks, in ISRM Symposium: Eurock , vol. 92, pp. 57 62. Li, J., Q. Y. Duan, W. Gong, A. Ye, Y. Dai, C. Miao, Z. Di, C. Tong, and Y. Sun (2013), Assessing parameter importance of the Common Land Model based on qualitative and quantitative sensitivity analysis, Hydrol. Earth Syst. Sci. , 17 (8), 3279 3293, doi:10.5194/hess 17 3279 2013. Lichtner, P. C. (1988), The quasi stationary state approximation to coupled mass transport and fluid rock interaction in a porous medium, Geochim. Cosmochim. Acta , 52 (1), 143 165. Liedl, R. (2003), Simulation of the development of karst aquifers using a coupled continuum pipe flow model, Water Resour. Res. , 39 (3), 1 11, doi:10.1029/2001WR001206. Littlefield, J. R., M. A. Culbreth, S . B. Upchurch, and M. T. Stewart (1984), Relationship of modern sinkhole development to large scale photolinear features, in Proceedings of First Multidisciplinary Conference on Sinkholes, Orlando, Florida , pp. 189 195. Luz, B., and E. Barkan (2011), The i sotopic composition of atmospheric oxygen, Global Biogeochem. Cycles , 25 (3), 1 14, doi:10.1029/2010GB003883. Macalady, J. L., E. H. Lyon, B. Koffman, L. K. Albertson, K. Meyer, S. Galdenzi, and S. Mariani (2006), Dominant microbial populations in limestone corroding stream biofilms, Frasassi cave system, Italy, Appl. Environ. Microbiol. , 72 (8), 5596 5609, doi:10.1128/AEM.00715 06. Maramathas, A. J., and A. G. Boudouvis (2006), Manifestation and measurement of the fractal characteristics of karst hydrogeolog ical formations, Adv. Water Resour. , 29 (1), 112 116.
193 Marella, R. L. (2008), Water use in Florida, 2005 and trends 1950 2005 , Geological Survey (US). Mariotti, A. (1983), Atmospheric nitrogen is a reliable standard for natural 15N abundance measurements, Na ture , 303 , 685 687. Martin, J. B., and R. W. Dean (2001), Exchange of water between conduits and matrix in the Floridan aquifer, Chem. Geol. , 179 (1 4), 145 165, doi:10.1016/S0009 2541(01)00320 5. Martin, J. B., and S. L. Gordon (2000), Surface and ground w ater mixing, flow paths, and temporal variations in chemical compositions of karst springs, Groundw. flow Contam. Transp. carbonate aquifers. AA Balkema , 65 92. Martin, J. B., M. Kurz, and M. Khadka (n.d.), Climate Control of Decadal SCale Increases in Apa rent Ages of Eogenetic Karst Spring Water, Hydrogeol. J. Mazzilli, N., V. Guinot, and H. Jourde (2012), Sensitivity analysis of conceptual model calibration to initialisation bias. Application to karst spring discharge models, Adv. Water Resour. , 42 , 1 16, doi:10.1016/j.advwatres.2012.03.020. McClain, M. E. et al. (2003), Biogeochemical Hot Spots and Hot Moments at the Interface of Terrestrial and Aquatic Ecosystems, Ecosystems , 6 (4), 301 312, doi:10.1007/s10021 003 0161 9. Moore, P. J., J. B. Martin, and E . J. Screaton (2009), Geochemical and statistical evidence of recharge, mixing, and controls on spring discharge in an eogenetic karst aquifer, J. Hydrol. , 376 (3 4), 443 455, doi:10.1016/j.jhydrol.2009.07.052. Moore, P. J., J. B. Martin, E. J. Screaton, an d P. S. Neuhoff (2010), Conduit enlargement in an eogenetic karst aquifer, J. Hydrol. , 393 (3 4), 143 155, doi:10.1016/j.jhydrol.2010.08.008. Morris, J. T., G. J. Whiting, and F. H. Chapelle (1988), Potential denitrification rates in deep sediments from the southeastern coastal plain, Environ. Sci. Technol. , 22 (7), 832 836. Morris, M. D. (1991), Factorial Sampling Plans for Preliminary Computational Experiments, Technometrics , 33 (2), 161 174, doi:10.2307/1269043. Mulholland, A. P. J. et al. (2016), Can uptak e length in streams be determined by
194 Munch, D. A., D. J. Toth, C. Huang, J. B. Davis, C. M. Fortich, W. L. Osburn, E. J. Phlips, E. L. Quinlan, M. S. Allen, and M. J. Woods (2006), Fifty year retrospective study of the ecology of Silver Springs, Florida, Publ. Number SJ2007 SP4. St. Johns River Wat er Manag. Dist. Palatka, FL . MuÃ±oz Carpena, R., Z. Zajac, and Y. M. Kuo (2007), G lobal sensitivity and uncertainty analyses of the water quality model VFSMOD W , , 50 (5), 1719 1732. uring Nutrient Spiralling in Streams, Can. J. Fish. Aquat. Sci. , 38 (7), 860 863, doi:10.1139/f81 114. North, N. N., S. L. Dollhopf, L. Petrie, J. D. Istok, D. L. Balkwill, and J. E. Kostka (2004), Change in bacterial community structure during in situ bios timulation of subsurface sediment cocontaminated with uranium and nitrate, Appl. Environ. Microbiol. , 70 (8), 4911 4920, doi:10.1128/AEM.70.8.4911 4920.2004. Nossent, J., and W. Bauwens (2012), Multi variable sensitivity and identifiability analysis for a c omplex environmental model in view of integrated water quantity and water quality modeling, Water Sci. Technol. , 65 (3), 539 549, doi:10.2166/wst.2012.884. Odum, H. T. (1956), Primary production in flowing waters, Limnol. Ocean. , 1 (2), 102 117. Otte, S., J. G. Kuenen, L. P. Nielsen, H. W. Paerl, J. Zopfi, H. N. Schulz, A. Teske, B. Strotmann, V. A. Gallardo, and B. B. JÃ¸rgensen (1999), Nitrogen, carbon, and sulfur metabolism in natural Thioploca samples, Appl. Environ. Microbiol. , 65 (7), 3148 3157. Ottley, C . J., W. Davison, and W. M. Edmunds (1997), Chemical catalysis of nitrate reduction by iron (II), Geochim. Cosmochim. Acta , 61 (9), 1819 1828, doi:10.1016/S0016 7037(97)00058 6. Padilla, A., A. Pulido Bosch, and A. Mangin (1994), Relative Importance of Base flow and Quickflow from Hydrographs of Karst Spring, Ground Water , 32 (2), 267 277, doi:10.1111/j.1745 6584.1994.tb00641.x. Palmer, A. N. (1991), Origin and development of limestone caves, Geol. Soc. Am. Bull. , 103 (2), 1 21, doi:10.1177/030913338100500204. Palmer, A. N. (2003), Speleogenesis in carbonate rocks, Speleogenes. Evol. Karst Aquifers , 1 (1), 1 11. Palmer, A. N. (2006), Digital modeling of karst aquifers successes, failures, and promises, Geol. Soc. Am. Spec. Pap. , 404 , 243 250.
195 Pan, F., J. Zhu, M. Ye, Y. A. Pachepsky, and Y. S. Wu (2011), Sensitivity analysis of unsaturated flow and contaminant transport with correlated parameters, J. Hydrol. , 397 (3), 238 249. Pandey, S. N., A. Chaudhuri, S. Kelkar, V. R. Sandeep, and H. Rajaram (2014), Investigati on of permeability alteration of fractured limestone reservoir due to geothermal heat extraction using three dimensional thermo hydro chemical (THC) model, Geothermics , 51 , 46 62. Pappenberger, F., K. J. Beven, M. Ratto, and P. Matgen (2008), Multi method global sensitivity analysis of flood inundation models, Adv. Water Resour. , 31 (1), 1 14, doi:10.1016/j.advwatres.2007.04.009. Pardo IgÃºzquiza, E., P. A. Dowd, C. Xu, and J. J. DurÃ¡n Valsero (2012), Stochastic simulation of karst conduit networks, Adv. Wate r Resour. , 35 , 141 150, doi:10.1016/j.advwatres.2011.09.014. Peaceman, D. W. (1978), Interpretation of well block pressures in numerical reservoir simulation (includes associated paper 6988), Soc. Pet. Eng. J. , 18 (03), 183 194. Peaceman, D. W. (1983), Inte rpretation of well block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability, Soc. Pet. Eng. J. , 23 (03), 531 543. in trans ition from pressurized flow to free surface flow, Hydrol. Earth Syst. Sci. , 18 (11), 4617 4633. Perrin, J., P. Y. Jeannin, and F. Cornaton (2007), The role of tributary mixing in chemical variations at a karst spring, Milandre, Switzerland, J. Hydrol. , 332 ( 1 2), 158 173, doi:10.1016/j.jhydrol.2006.06.027. Phelps, B. G. G., and U. S. G. Survey (2004), Chemistry of Ground Water in the Silver Springs Basin , Florida , with an U . S . Department of the Interior, Ground Water . Phillips, R., J. Kondev, J. Theriot, H. Garcia, and B. Chasan (2010), Physical Biology of the Cell, Am. J. Phys. , 78 (11), 1230, doi:10.1119/1.3459039. Pinay, G., L. Roques, and A. Fabre (1993), Spatial and temporal patterns of denitrification in a riparian forest, J. Appl. Ecol. , 30 (4), 581 591, doi:10.2307/2404238. Pittman, J. R., H. H. Hatzell, and E. T. Oaksford (1997), Spring contributions to water quantity and nitrate loads in the Suwannee River during base flow in July 1995, Water Resources Investig. Rep , 97 , 4152. Rabalais, N. N. (2002 ), Nitrogen in aquatic ecosystems., Ambio , 31 (2), 102 112, doi:10.2307/4315222.
196 Raeisi, E., C. Groves, and J. Meiman (2007), Effects of partial and full pipe flow on hydrochemographs of Logsdon river, Mammoth Cave Kentucky USA, J. Hydrol. , 337 (1), 1 10. Ra ndazzo, A. F. (1997), The sedimentary platform of Florida: Mesozoic to Cenozoic, Geol. Florida. Univ. Press Florida, Gainesv. , 39 56. Randazzo, A. F., and J. I. Bloom (1985), Mineralogical changes along the freshwater/saltwater interface of a modern aquife r, Sediment. Geol. , 43 (1), 219 239. Reeves, D. M., D. A. Benson, M. M. Meerschaert, and H. P. Scheffler (2008), Transport of conservative solutes in simulated fracture networks: 2. Ensemble solute transport and the correspondence to operator stable limit d istributions, Water Resour. Res. , 44 (5), 1 20, doi:10.1029/2008WR006858. Reimann, T., T. Geyer, W. B. Shoemaker, R. Liedl, and M. Sauter (2011a), Effects of dynamically variable saturation and matrix conduit coupling of flow in karst aquifers, Water Resour . Res. , 47 (11), 1 19, doi:10.1029/2011WR010446. Reimann, T., C. Rehrl, W. B. Shoemaker, T. Geyer, and S. Birk (2011b), The significance of turbulent flow representation in single continuum models, Water Resour. Res. , 47 (9), 1 15, doi:10.1029/2010WR010133. Revsbech, N. P., P. B. Christensen, L. P. Nielsen, and J. S??rensen (1989), Denitrification in a trickling filter biofilm studied by a microsensor for oxygen and nitrous oxide, Water Res. , 23 (7), 867 871, doi:10.1016/0043 1354(89)90011 0. Rivett, M. O., J. W. N. Smith, S. R. Buss, and P. Morgan (2007), Nitrate occurrence and attenuation in the major aquifers of England and Wales, Q. J. Eng. Geol. Hydrogeol. , 40 (4), 335 352, doi:10.1144/1470 9236/07 032. Rivett, M. O., S. R. Buss, P. Morgan, J. W. N. Smith, and C. D. Bemment (2008), Nitrate attenuation in groundwater: A review of biogeochemical controlling processes, Water Res. , 42 (16), 4215 4232, doi:10.1016/j.watres.2008.07.020. Robertson, W. D., B. M. Russell, and J. A. Cherry (1996), Attenuation of nitrat e in aquitard sediments of southern Ontario, J. Hydrol. , 180 (1 4), 267 281, doi:10.1016/0022 1694(95)02885 4. Rosenau, J. C., G. L. Faulkner, C. W. Hendry Jr, and R. W. Hull (1977), Springs of Florida , Florida Department of Natural Resources: Bureau of Geo logy. Saller, S. P., M. J. Ronayne, and A. J. Long (2013), Comparison of a karst groundwater model with and without discrete conduit flow, Hydrogeol. J. , 21 (7), 1555 1566, doi:10.1007/s10040 013 1036 6.
197 Saltelli, A., S. Tarantola, and K. S. Chan (1999), A quantitative model independent method for global sensitivity analysis of model output, Technometrics , 41 (1), 39 56. Saltelli, A., S. Tarantola, F. Campolongo, and M. Ratto (2004), Sensitivity analysis in practice: a guide to assessing scientific models , J ohn Wiley & Sons. Satoh, H., H. Ono, B. Rulin, J. Kamo, S. Okabe, and K. I. Fukushi (2004), Macroscale and microscale analyses of nitrification and denitrification in biofilms attached on membrane aerated biofilm reactors, Water Res. , 38 (6), 1633 1641, doi :10.1016/j.watres.2003.12.020. Scanlon, B. R., R. E. Mace, M. E. Barrett, and B. Smith (2003), Can we simulate regional groundwater flow in a karst system using equivalent porous media models? Case study, Barton Springs Edwards aquifer, USA, J. Hydrol. , 27 6 (1 4), 137 158, doi:10.1016/S0022 1694(03)00064 7. Schippers, A., L. N. Neretin, J. Kallmeyer, T. G. Ferdelman, B. a Cragg, R. J. Parkes, and B. B. JÃ¸rgensen (2005), Prokaryotic cells of the deep sub seafloor biosphere identified as living bacteria., Natu re , 433 , 861 864, doi:10.1038/nature03302. (1998), Spatial Variability in In Situ Aerobic Respiration and Denitrification Rates in a Petroleum Contaminated Aquifer, Gro und Water , 36 (6), 924 937. Schroth, M. H., J. D. Istok, and R. Haggerty (2000), In situ evaluation of solute retardation using single well push pull tests, Adv. Water Resour. , 24 (1), 105 117. Scott, T. M., G. H. Means, R. P. Meegan, R. C. Means, S. Upchurc h, R. E. Copeland, J. Jones, T. Roberts, and A. Willet (2004), Springs of Florida, Seitzinger, S., J. A. Harrison, J. K. B??hlke, A. F. Bouwman, R. Lowrance, B. Peterson, C. Tobias, and G. Van Drecht (2006), Denitrification across landscapes and waterscape s: A synthesis, Ecol. Appl. , 16 (6), 2064 2090, doi:10.1890/1051 0761(2006)016[2064:DALAWA]2.0.CO;2. Sepulveda, A. A., B. G. Katz, and G. L. Mahon (2006), Potentiometric surface of the Upper Floridan aquifer in the Ichetucknee springshed and vicinity, north ern Florida, September 2003 . Sexstone, A. J., N. P. Revsbech, T. B. Parkin, and J. M. Tiedje (1985), Direct Measurement of Oxygen Profiles and Denitrification Rates in Soil Aggregates1, Soil Sci. Soc. Am. J. , 49 , 645, doi:10.2136/sssaj1985.0361599500490003 0024x. Shin, M. J., J. H. A. Guillaume, B. F. W. Croke, and A. J. Jakeman (2013), Addressing ten questions about conceptual rainfall runoff models with global sensitivity analyses in R, J. Hydrol. , 503 , 135 152, doi:10.1016/j.jhydrol.2013.08.047.
198 Shoemaker , W. B., E. L. Kuniansky, S. Birk, S. Bauer, and E. D. Swain (2005), Documentation of a Conduit Flow Process (CFP) for MODFLOW 2005, U.S. Geol. Surv. Tech. Methods, B. 6, chapter A24 , 50. Shoemaker, W. B., K. J. Cunningham, E. L. Kuniansky, and J. Dixon (2 008), Effects of turbulence on hydraulic heads and parameter sensitivities in preferential groundwater flow layers, Water Resour. Res. , 44 (3), 1 11, doi:10.1029/2007WR006601. Siemers, J., and W. Dreybrodt (1998), Early development of karst aquifers on pers olation networks of fractures in limestone, Water Resour. Res. , 34 (3), 409 419. SJRWMD (2014), Intermediate Confining Unit Thickness, Smith, R. L., and J. H. Duff (1988), Denitrification in a sand and gravel aquifer., Appl. Environ. Microbiol. , 54 (5), 1071 1078. Smith, R. L., J. K. Bohlke, B. Song, and C. R. Tobias (2015), Role of anaerobic ammonium oxidation (anammox) in nitrogen removal from a freshwater aquifer, Environ. Sci. Technol. , 49 (20), 12169 12177. Snider, D. M., J. Spoelstra, S. L. Schiff, and J . J. Venkiteswaran (2010), Stable oxygen isotope ratios of nitrate produced from nitrification: 18O labeled water incubations of agricultural and temperate forest soils, Environ. Sci. Technol. , 44 (14), 5358 5364, doi:10.1021/es1002567. Sobol, I. M. (2001), Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Math. Comput. Simul. , 55 (1), 271 280. Song, B., J. A. Lisa, and C. R. Tobias (2014), Linking DNRA community structure and activity in a shallow lagoonal estuarin e system, Front. Microbiol. , 5 (SEP), doi:10.3389/fmicb.2014.00460. Srivastava, V., W. Graham, R. MuÃ±oz Carpena, and R. M. Maxwell (2014), Insights on geologic and vegetative controls over hydrologic behavior of a large complex basin Global Sensitivity An alysis of an integrated parallel hydrologic model, J. Hydrol. , 519 (PB), 2238 2257, doi:10.1016/j.jhydrol.2014.10.020. Sun, W., C. Xia, M. Xu, J. Guo, A. Wang, and G. Sun (2014), Diversity and distribution of planktonic anaerobic ammonium oxidizing bacteria in the Dongjiang River, China, Microbiol. Res. , 169 (12), 897 906, doi:10.1016/j.micres.2014.05.003. Svensson, U., and W. Dreybrodt (1992), Dissolution kinetics of natural calcite minerals in CO 2 water systems approaching calcite equilibrium, Chem. Geol. , 100 (1), 129 145.
199 Swamee, P. K., and N. Swamee (2007), Full range pipe flow equations, J. Hydraul. Res. , 45 (6), 841 843. Szymczak, P., and A. J. C. Ladd (2009), Wormhole formation in dissolving fractures, J. Geophys. Res. Solid Earth , 114 (B6). Szymczak, P ., and A. J. C. Ladd (2011), The initial stages of cave formation: Beyond the one dimensional paradigm, Earth Planet. Sci. Lett. , 301 (3), 424 432. Taylor, C. J., and E. a Greene (2005), Hydrogeologic Characterization and Methods Used in the Investigation o f Karst Hydrology, F. Tech. Estim. Water Fluxes Between Surf. Water Gr. Water , 75 114, doi:PNR61. Tesoriero, A. J., H. Liebscher, and S. E. Cox (2000), Mechanism and rate of denitrification in an agricultural watershed: Electron and mass balance along grou ndwater flow paths, Water Resour. Res. , 36 (6), 1545 1559, doi:10.1029/2000WR900035. Teutsch, G., and M. Sauter (1991), Groundwater modeling in karst terranes: Scale effects, data acquisition and field validation, Proc. Third Conf. Hydrogeol. Ecol. Monit. M anag. Gr. Water Karst Terranes, Nashville, TN , 17 35, doi:PNR61. Thayalakumaran, T., K. L. Bristow, P. B. Charlesworth, and T. Fass (2008), Geochemical conditions in groundwater systems: Implications for the attenuation of agricultural nitrate, Agric. Wate r Manag. , 95 (2), 103 115, doi:10.1016/j.agwat.2007.09.003. Tiedje, J. M. (1988), Ecology of denitrification and dissimilatory nitrate reduction to ammonium, Environ. Microbiol. Anaerobes , 179 244. Toth, D. J., and B. G. Katz (2006), Mixing of shallow and d eep groundwater as indicated by the chemistry and age of karstic springs (Hydrogeology Journal (2006) (827 847) 10.1007/s10040 005 0478 x)), Hydrogeol. J. , 14 (6), 1060 1080, doi:10.1007/s10040 006 0099 z. Ugolini, F., R. Henneberger, H. BÃ¼rgmann, J. Zeyer, and M. H. Schroth (2014), In Situ Sonication for Enhanced Recovery of Aquifer Microbial Communities, Groundwater , 52 (5), 737 747. Upchurch, S. B., and J. R. Littlefield Jr (1988), Evaluation of data for sinkhole development risk models, Environ. Geol. Water Sci. , 12 (2), 135 140. USFW S (2014), Florida Wetlands , US Fish and Wildlife Service. USGS (1998), Floridan Aquifers , U.S. Geological Survey, Reston, Va.
200 van Griensven, A., T. Meixner, S. Grunwald, T. Bishop, M. Diluzio, and R. Srinivasan (2006), A global sensitivity analysis tool f or the parameters of multi variable catchment models, J. Hydrol. , 324 (1 4), 10 23, doi:10.1016/j.jhydrol.2005.09.008. Vernon, R. O. (1951), Geology of Citrus and Levy Counties, Florida , Tallahassee, Fl. Vidon, P., and A. R. Hill (2004), Denitrification and patterns of electron donors and acceptors in eight riparian zones with contrasting hydrogeology, Biogeochemistry , 71 (2), 259 283, doi:10.1007/s10533 004 9684 1. Vogel, J. C., A. S. Talma, and T. H. E. Heaton (1981), Gaseous nitrogen as evidence for denitr ification in groundwater, J. Hydrol. , 50 , 191 200. Wang, L., and H. Liu (2006), An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling, Int. J. Geogr. Inf. Sci. , 20 (2), 193 213. Weber, K. A., M. M. Urrutia, P. F. Churchill, R. K. Kukkadapu, and E. E. Roden (2006), Anaerobic redox cycling of iron by freshwater sediment microorganisms, Environ. Microbiol. , 8 (1), 100 113, doi:10.1111/j.1462 2920.2005.00873.x. White, W. B. (1977), Ro le of solution kinetics in the development of karst aquifers, Karst Hydrogeol. Int. Assoc. Hydrogeol. 12th Mem. , 503 517. White, W. B. (1993), Analysis of karst aquifers , Van Nostrand Reinhold. White, W. B., D. C. Culver, J. S. Herman, T. C. Kane, and J. E . Mylroie (1995), Karst lands, Am. Sci. , 83 (5), 450 459. Worthington, S. R. H. (1999), A comprehensive strategy for understanding flow in carbonate aquifers, Karst Model. Spec. Publ. , 5 , 30 37. Worthington, S. R. H. (2009), Diagnostic hydrogeologic charact eristics of a karst aquifer (Kentucky, USA), Hydrogeol. J. , 17 (7), 1665 1678, doi:10.1007/s10040 009 0489 0. Worthington, S. R. H., and D. C. Ford (2009), Self organized permeability in carbonate aquifers, Ground Water , 47 (3), 326 336, doi:10.1111/j.1745 6584.2009.00551.x. Xu, Z., B. X. Hu, H. Davis, and S. Kish (2015), Numerical study of groundwater flow cycling controlled by seawater/freshwater interaction in a coastal karst aquifer through conduit network using CFPv2, J. Contam. Hydrol. , 182 , 131 145, doi:10.1016/j.jconhyd.2015.09.003.
201 BIOGRAPHICAL SKETCH Wesley Henson grew up tumbling around the Nevada Desert; where water is understood to be a precious resource. His first graduate deg ree was a Master of Science in h y drogeology in the fall of 2008. Wanting to explore other side of the water cycle, he continued his graduate education with a Doctor of Philosophy in agricultural and biological e ngineering in Florida. His research explores the interface between hydrologic and biogeochemical processes, examining process and model complexity in human engineered and natural environments. While simultaneously completing his graduate work, Wesley has been a Hydrologist for the U.S. Geological Survey where he develops tools to be tter manage water resources and develop parsimonious integrated hydrologic models. Wesley obtained his doctorate in summer 2016.