UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 INVESTIGATION OF GROUNDWATER FLOW IN KARST AQUIFER AND APPLICATION OF GEOMORPHOLOGIC KARST MODEL By LILI YU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORID A IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2010 PAGE 2 2 2010 Lili Yu PAGE 3 3 To all who nurtured my intellectual curios ity, academic interests, and sense of scholarship throughout my study, making this research possible PAGE 4 4 ACKNOWLEDGMENTS I would like to thank my advisor, Dr. Kirk Hatfield, for providing me the opportunity to pursue my Ph.D. degree and his continu ed support throughout my graduate career. I appreciate his invaluable guidanc e, encour agement, patience, generosity, concern and sense of humor. I would like to thank Dr. Ronald Green for gi ving me the opportunity to work on this project and his advice and patienc e throughout my research. I would like to thank Dr. Louis Motz His advice and encouragement were invaluable to me. I would especially like to thank Dr. Ma rk Newman. He gave me many important suggestions for my research and was a great source of encouragement and camaraderie. I would like to thank Dr. Alexandru Sher emet and Dr. Jonathan Martin for their technical guidance. I would like to thank my family, my gr andparents, my parents and my husband, for their love and support. And of course, I cannot forget to thank Nancy Been, Carolyn Carpenter, Doretha Ray, and all of the civil engineering staff for their help over the years. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDG MENTS .................................................................................................. 4LIST OF TABLES............................................................................................................ 7LIST OF ABBR EVIATIONS........................................................................................... 13ABSTRACT................................................................................................................... 14 CHA PTER 1 INTRODUCTION.................................................................................................... 162 LITERATURE REVIEW .......................................................................................... 212.1 Application of Time Series Techniques in Karst An alysis................................. 222.2 Concepts of Geomorphologic Stru cture Model................................................. 272.3 Numerical Models of Groundwater Flow in a Ka rst Aquifer............................... 363 DERIVATION OF GEOMORPHOL OGIC KARST FL OW M ODELS....................... 403.1 Derivation of a Third Order Geomorphologic Karst Model................................ 403.2 Derivation of a Fourth Or der Geomorphologi c Karst Model.............................. 443.3 Derivation of The Fifth Or der Geomorphologic Karst Model............................. 534 APPLICATION OF GEOMORPH OLOGIC KARST MODEL ................................... 624.1 Application of DCM to Construct a Hypothetical Karst A quifer.......................... 624.2 Cross Correlation Analysis of t he Hypothetical Kars t Flow System.................. 654.3 Geomorphologic Karst Flow Modeling wi th Matrix Flow as the First Flow Segm ent......................................................................................................... 704.3.1 Calibration of the Fourth Order Geomorphologic Karst Flow Model...... 704.3.2 Analysis of RL and V in the Fourth Order Geomorphologic Karst Model................................................................................................... 784.4 Treatment of the Por ous Medium as the First Two Flow Segments of a Geomorpholog ic Model.................................................................................. 824.4.1 Calibration of the Fifth Order Geomorphologic Model of the Hypothetical Aquifer............................................................................ 834.5 Analysis of the Travel Time of Rainfall in the Geomorphologic Models............ 854.5.1 Travel Time Analysis in the Fourth Order Geomorphologic Karst Model................................................................................................... 864.5.2 Travel Time Analysis for t he Fifth Order Geomorphologic Model.......... 884.5.3 Application of Laplace Transfo rms to Estimate the Travel Time........... 91 PAGE 6 6 5 GEOMORPHOLOGIC KARST MODEL WITH COMMON CONDUIT FLOW FOR DIFFERENT FL OW COMP ONENTS ............................................................. 955.1 The Conceptual Model of Common Conduit Flow for Different Flow Components................................................................................................... 955.2 Calibration of Common C onduit Geomorphol ogic Model................................ 1025.3 Application of the Common Condui t Geomorphologic Model to a Revised Conduit Ne twork........................................................................................... 1065.3.1 The Revised Conduit Ne twork............................................................. 1065.3.2 Cross Correlation Analysis of Rainfall with Calculated Spring Dischar ge.......................................................................................... 1085.3.3 Calibrate Common Conduit G eomorphologic Karst Model for the Revised Condui t Network.................................................................. 1146 ANALYSES OF OBSERVED DATA IN THE BLUE SPRING AREA AND GEOMOR PHOLOGIC MODEL CALI BRATION.................................................... 1196.1 Observation Data Anal ysis for Blue Spring Area............................................ 1196.2 Cross Correlation Analysis of Data from the Bl ue Spring Area....................... 1216.2.1 Autocorrela tion Anal ysis...................................................................... 1216.2.2 Correlation Analysis of Ra infall with Blue Sp ring Disc harge................ 1226.2.3 Correlation Analysis of Rain fall Versus Groundwater Levels.............. 1336.2.4 Correlation Analysis of Groundw ater Levels versus Blue Spring Dischar ge.......................................................................................... 1396.2.5 Di scussion........................................................................................... 1436.3 Calibration and Validation of the Geomorphologic Karst Model of Blue Spri ng........................................................................................................... 1496.4 Analysis of the Travel Time of Rainfall in the Geomorphologic Models.......... 1536.5 Application of Laplace Transforms to Estimate t he Travel Time..................... 1547 CONCLUSIONS AND RE COMMENDAT IONS ..................................................... 1577.1 Conc lusions.................................................................................................... 1577.2 Recomm endations.......................................................................................... 159APPENDIX: PARAMETERS FOR EQ UATION 340 .................................................... 160LIST OF RE FERENCES............................................................................................. 171BIOGRAPHICAL SKETCH.......................................................................................... 174 PAGE 7 7 LIST OF TABLES Table page 41 Hydraulic conductivities of conduits in the hypothetic al karst model ................... 64 42 The calibrated results of 2V for different 2 R L when 1,1,12 VrtRLrt are constant.............................................................................................................. 80 43 The calibrated results of 1V for different 1 R L when 2,1,22 VrtRLrt are constant.............................................................................................................. 81 44 Values of ij p in fourth order geomorphologic model............................................ 88 45 Probabilities of the four flow paths for the four th order geomorphologic model.. 88 46 Estimation of travel time of rain fall in fourth order geomorphologic model......... 89 47 Values of ij p in the fifth order geomorphologic model......................................... 91 48 Probabilities of the flow paths for the fift h order geomorphol ogic model............. 91 49 Estimation of travel time of rain fall in fifth order geomorphologic model............. 92 410 Central moments of I UH for the two flow components calculated from the fourth order geom orphologic model.................................................................... 94 51 Calibrated results of Vc for different R L for fifth order geomorphologic karst model................................................................................................................ 104 52 Calibrated results of Vc for different R L for fifth order geomorphologic karst model................................................................................................................ 115 61 Results of cross correlation analysis of rainfall at Deland with Blue Spring dischar ge.......................................................................................................... 145 62 Results of cross correlation analysis of rainfall at Deland and groundwater levels................................................................................................................ 148 63 Calibrated results of Vc for different R L............................................................ 153 64 Analysis mean travel time of different fl ow com ponents................................... 157 65 Central moments of I UH for the two flow components calculated from the common conduit model for Blue Spring............................................................ 158 PAGE 8 8 LIST OF FIGURES Figure page 11 Schematic drawing of the commo n types of aquifer m edia................................ 19 12 Epikar st Zone..................................................................................................... 20 21 Two boxes models for carbonate aqui fers, with fraction of storage (in percentage) and residence time in each box and flow (m3s1) between boxes (Worthington 2007)............................................................................................. 25 22 Third order surface basin with orderi ng system (I. Rodriguez, et al. (1979))...... 30 23 Common patterns of c onduits (N.Pal mer 1991).................................................. 32 24 Representation of a third order basin as a continuous Markov process (I. Rodriguez, et al. (1979))..................................................................................... 35 31 Considering a karst aquifer as a third order geomorphologic flow system.......... 43 32 Conceptual model of a third order karst aquifer using five model flow segments identified in par entheses or as standalone numbers (e.g., 1,2, 3,4,5).................................................................................................................. 44 33 Conceptual model of a fourth order karst aquifer using six model flow segments identified in parentheses or as standalone numbers (i.e., 1,2,3,4, 5,6)......................................................................................................... 46 34 Conceptual model of a fifth order geomorphologic Karst system using seven model flow segments identified in par entheses or as standalone numbers (i.e., 1,2,3, 4,5,6,7).............................................................................................. 56 41 Conduits network in the hypothetical ka rst aquifer............................................. 63 42 Daily rainfall data used in the hypothetical karst model...................................... 64 43 Spring discharge calculated by the DCM............................................................ 65 44 Cross correlation of rainfall and sp ring discharge of the hypothetical Karst aquife r................................................................................................................ 66 45 Smoothed spring discharge (calculated spring discharge deduct fluctuations).. 67 46 Fluctuations of spring di scharge......................................................................... 67 47 Cross correlation of rainfall wit h fluctuations of spring di scharge....................... 68 48 Cross correlation of rainfall versus smoothed spring di scharge......................... 68 PAGE 9 9 49 Scatter plot of rainfall and fl uctuations of spring dischar ge................................. 69 410 Order numbers of each conduit in the hypothetical Karst aquifer....................... 71 411 Comparison of smoothed spring discharge data from DCM and discharges simulated from the fourth order geomorphologic Karst model assuming one flow co mponent.................................................................................................. 72 412 Comparison of results with different RB values.................................................. 73 413 Comparison of spring discharge from the DCM and fourth order two flow components geomorphologic Karst model in the calib ration period.................... 75 414 IUH of the two flow components calculated from the fourth order two flow components geomorpholog ic Kars t model......................................................... 75 415 Spring discharge from the two flow components of the geom orphologic Karst model.................................................................................................................. 76 416 Comparison of DCM simulated spri ng discharges to discharges obtained by combining simulations from t he fourthordertwocomponents geomorphologic karst model and shortterm fluctuations from equation 41...... 76 417 Comparison of spring discharge from the DCM and results from GEO model for the validat ion per iod...................................................................................... 78 418 Cross correlation of rainfall wit h spring discharge from the first flow component .......................................................................................................... 79 419 Cross correlation of rainfall wit h spring discharge from the second flow component .......................................................................................................... 79 420 Plot of V2 versus RL2 when the values of V1, rt1, RL1, rt2 are held constant.... 81 421 Plot of V1 versus RL1 when the val ues of V2, rt1, RL2, rt2 are constant........... 82 422 Comparison of results from different values of RL and V of fourth order geomorphologic model....................................................................................... 83 423 Order numbers of each conduit in the hypothetical aquifer with considering the porous medium as the first two flow segments............................................. 84 424 Comparison of smoothed spring disc harge resulted from DCM and fifth order geomorphologic Karst model with one flow component...................................... 85 425 Four possible paths in fourth order geomor phologic model................................ 87 426 Possible paths of the fi fth order geomorphologic m odel..................................... 88 PAGE 10 10 51 Two different flow components t hat share a common c onduit network.............. 97 52 Conceptual models for the slow (a) and fast (b) components of a geomorphologic karst flow model usi ng a common karst conduit network......... 98 53 Plot of Vc versus RL......................................................................................... 104 54 Comparison of results from the new geomorphologic model and DCM............ 105 55 Instantaneous unit hydrographs of the two fl ow com ponent............................. 106 56 Simulated spring discharge fr om the two flow components.............................. 107 57 Validation of the Geom orphologic model using common conduits for different flow co mponents............................................................................................... 108 58 Revised c onduit net work.................................................................................. 109 59 Daily rainfall data used in the hypothetical karst model.................................... 109 510 Spring discharge calculated by the DCM.......................................................... 110 511 Cross correlation between rainfall and spring discharge of the hypothetical Karst groundwat er system................................................................................ 110 512 Smoothed spring discharge in the calibrati on period........................................ 111 513 Fluctuations of spring disc harge in the calib ration period................................. 112 514 Cross correlation of rainfall with Fluctuations of spring di scharge.................... 112 515 Cross correlation of rainfall with smoothed spring disc harge............................ 113 516 Scatter plot of rainfall and fl uctuations of spring discharge............................... 113 517 Plot of Vc versus RL......................................................................................... 116 518 Comparison of results fr om geomorphologic model wit h DCM......................... 116 519 Instantaneous unit hydrographs of the two fl ow com ponent............................. 117 520 Spring discharge from t he two flow components.............................................. 118 521 Validation of the Geom orphologic model using common conduits for different flow co mponents............................................................................................... 119 61 Study Area........................................................................................................ 122 62 Available observed data of rainfall and spring discharge in Blue Spring area.. 123 PAGE 11 11 63 Available data of groundw ater levels at four observation wells and rainfall ...... 123 64 Autocorrelation functions of rainfall, water level and spring di scharge............. 124 65 Cross correlation analysis of rainfa ll with spring discharge from 12/8/2001 to 9/7/ 2004........................................................................................................... 126 66 Spectral dens ity of spring................................................................................. 127 67 Coherency and P hase diagr ams...................................................................... 128 68 Cross correlation analysis of rainfa ll with spring discharge from 3/21/2005 to 8/21/ 2006......................................................................................................... 129 69 Monthly averaged rainfall dat a from January 2001 to Ma y 2006...................... 130 610 Cross correlation analysis of ra infall with spring discharge in 2002.................. 130 611 Cross correlation analysis of ra infall with spring discharge in 2003.................. 131 612 Cross correlation analysis of ra infall with spring discharge in 2004.................. 131 613 Cross correlation analysis of rain fall with spring discharge in 2005.................. 132 614 Cross correlation analysis of ra infall with spring discharge in 2006.................. 132 615 Cross correlation analysis of rainfa ll with spring discharge from 6/23/2002 to 9/26/ 2002......................................................................................................... 133 616 Cross correlation analysis of rainfa ll with spring discharge from 7/26/2004 to 9/7/ 2004. .......................................................................................................... 134 617 Cross correlation analysis of rainfa ll with spring discharge from 7/6/2005 to 9/10/ 2005. ........................................................................................................ 134 618 Cross correlation of groundwater leve ls at the three observation wells with the well V0083................................................................................................. 136 619 CCF of rainfall versus GW levels for the period 1/1/2004 to 6/30/2006............ 137 620 Cross correlation of rainfall with GW levels fo r the year 2004.......................... 138 621 CCF of rainfall with GW levels for t he year 2005.............................................. 139 622 Cross correlation analysis of rainfall with GW levels from 6/3/2004 to 10/15/2004 (wet season in 2004)..................................................................... 139 623 Cross correlation analysis of rainfa ll with GW levels from 5/30/2005 to 11/2/2005 (wet s eason in 2005)....................................................................... 140 PAGE 12 12 624 CCF or rainfall with GW at four ob servation wells (5/24/2005 to 6/30/2006) .... 140 625 Available data of spring discharge and groundw ater levels.............................. 141 626 Cross correlation analysi s of GW levels with Blue Spring discharge from 1/1/2004 to 7/25/2004....................................................................................... 142 627 Cross correlation analysi s of GW levels with Blue Spring discharge from 7/26/2004 to 9/7/2004....................................................................................... 143 628 Cross correlation analysi s of GW levels with Bl ue Spring discharge from 5/24/2005 to 12/31/ 2005................................................................................... 143 629 Cross correlation analysi s of GW levels with Blue Spring discharge from 1/1/2006 to 8/21/2006....................................................................................... 144 630 Plot of Vc versus RL......................................................................................... 154 631 Comparison of observed spring discharge with results from the calibrated geomorphologic karst model in the calibrati on period....................................... 154 632 Comparison of observed spring discharge with results from the calibrated geomorphologic karst model in the validat ion period........................................ 155 PAGE 13 13 LIST OF ABBREVIATIONS ACF Autocorrelation Function CCF Cross Correlation Function IUH Instantaneous Unit Hydrograph CFS Cubic feet per second DCM Dual conductivity model BSEACD Barton Springs Edwards Aquifer Conservation District PAGE 14 14 Abstract of Dissertation Pr esented to the Graduate School of the University of Florida in Partial Fulf illment of the Requirements for t he Degree of Doctor of Philosophy INVESTIGATION OF GROUNDWATER FLOW IN KARST AQUIFER AND APPLICATION OF GEOMORPHOLOGIC KARST GROUNDWATER FLOW MODEL By Lili Yu May 2010 Chair: Kirk Hatfield Cochair: Alexandru Sheremet Major: Civil Engineering Karst aquifers are extremely heterogeneous systems, in which conduits are indeterminate, and recharge is spatially variable. Developing a physicallybased numerical model of karst systems is difficult at best. As an alternative several analytical geomorphologic karst models were devel oped based on the concept of a geomorphologic structure model for surface basins. These new models express the hydrologic response of a karst basin in terms of calibrated conduit Horton numbers and mean flow velocities. The geomorphologic ka rst flow models were validated against a numerical model MODFLOWDCM Geomorphologic karst models of different orders and complexity were compared. In addition, analytical relationships were derived from governing equations that explicitly express karst basin hydraulic response times to rainfall inputs as a function of conduit Horton number RL and a mean basin velocity. One of the several geomorphologic kars t models developed in this dissertation was calibrated and validated to simulate t he spring discharges at Blue Spring near Deland, Florida. Time series analysis te chniques were used to analyze observed data for rainfall, groundwater levels and spring di scharge in the Blue Spring area in Volusia PAGE 15 15 County, Florida. Results indicate that multiple flow components might exist in the karst aquifer, and recharge conditions appear to influence the hydrodynamic characteristics of the karst aquifer. Temporal lag interv al results generated from the time series analysis are consistent with response times estimated from the calibrated geomorphologic model. PAGE 16 16 CHAPTER 1 INTRODUCTION It is estimated that the Un ited States derives 40% of its potable wat er needs from karst aquifers (Quinlan and Ewers, 1989) and for the total population of the earth, it is 25% (Ford and Williams, 1989). Clearly, the importance of karst aquifers cannot be overstated; and, understanding their evolution, hydraulic behavior, critical functions in a watershed, and vulnerability to contami nation are topics of high importance. Karst is a terrain with distinctive l andforms and hydrolog y created from the dissolution of soluble carbonate rocks such as limestone and dolomite. Over millions of years, dissolution transforms the underlying car bonate rock aquifers from a laminar or diffusive dominated flow regime to a system where diffusive flow feeds a welldeveloped network of solution conduits wit h primarily turbulent flow. This conduit network flow is typically manifested at ground surface as spri ng discharge. When the water table falls beneath the level of surface streams, theses streams lose water to underlying developed cave systems. If th is condition persists, more and more surface drainage will be diverted underground, stream valleys will virtually disappear only to be replaced by closed basins or sinkholes. Thus, karst terrain is characterized by springs, sinkholes, and subsurface caves and conduits that produce highly productive aqui fers which are extremely vulnerable to contamination from surface sources (i.e., agriculture, urban storm water, etc.). The issue of vulnerability is extremely pertinent to the karst environment. Groundwater flow through unconsolidated aquifer media tends to move slowly and the interand intragranular surface area of the porous media acts to a limited extent attenuate contaminants (i.e. through sorption) [see Figur e 11]. For Karst a quifers, flow can be PAGE 17 17 quite rapid once it reaches a conduit; consequent ly, the short hydraulic residence times and limited opportunity for su rfacedependent attenuation leav e karst aquifers highly vulnerable to contamination. Figure 11. Schematic drawing of the co mmon types of aquifer media. Note the reduction in inter granular surface area from intergranular media to fr actured or solution enhanced. (1Kentucky Geological Survey 19972010.). Karst hydrogeology is defined by a network of interconnected fissures, fractures and condu its emplaced in a relatively lowpermeability rock matrix. Most of the groundwater flow and transport occurs through the network of openings, while most of the groundwater storage occurs in the matrix. As a result, most karst aquifers are highly heterogeneous and anisotropic. For many karst aquifers, a large perc entage of the water stored underground is perched, or suspended, above the main part of the aquifer in t he "epikarst" (Figure 12) (Alexander, 2003). The epikarst is situat ed beneath the topsoil and above the unaltered bedrock. Water in the epikarst is stored in enlarged joints and bedding planes, spaces around pieces of float (rocks that have been detached from the bedrock), porosity within residual chert rubble, and the sma ller conduits in the bedrock (1Kentucky Geological Survey 19972010.). During floods, the epikarst may recharge the saturated zone by direct and fast infiltration through a transient conduit connection that is established PAGE 18 18 between infiltration and phreatic zones. During other periods, the epikarst contributes delayed recharge by slow infiltration. Figure 12. Epikarst Zone The karst aquifer is commonly modeled as a branching of tributary conduits, which connect together to drain a groundwater basin and discharge to a perennial spring (Worthington, 2002). The system is roughly analogous to a roofedover creek. In comparison, granular or fractured bedro ck aquifers have no equivalent "underground river" or channel. For the most part the conduit locations are indeterminate, and as such various types of surface features, su ch dolines, fracture traces, depressed water tables, and sinkholes are used to guide t he construction of a conduit network for numerical flow transport modeling. Due to the complexity and heterogeneity of karst aquifers, the information gathered from classical borehole techniques of hydrologic testing and analysis (e.g., pump tests) are difficult to interpret. In addition karst aqui fers typically exhibit dual PAGE 19 19 groundwater flow regimes, that is, a slow (dif fuse) flow in porous medium that feeds fast flow in the conduit network. In particula r, groundwater modeling tools based on Darcy law alone and developed for unconsolidat ed porous media type aquifers cannot adequately accommodate both the rapid flow of groundwater through conduits and the slow flow and storage of groundwater in the matrix of karst aquifers. The objective of this study is to develop simple analytical models for investigating and simulating groundwater flow in karst aquifers that can use data that is readily available. The purpose of these models is to 1) provide re liable predictions of spring flow, and 2) generate additional hydrogeologic information on the conduit network that could facilitate the development of comple x numerical models. Chapter 2 presents a literature review on modeling techniques to be used to construct the analytical karst flow models. In Chapter 3, the development of three differ ent analytical, geomorphologic, karst flow models is presented. In Chapt er 4 MODFLOWDCM is used to create a spring discharge data set for a hypothetical karst aquifer within known hydrogeologic parameters. Then, two analyt ical geomorphologic models are calibrated and validated using the hypothetical spring discharge data. Results from both models are compared. Based on the findings from Chapter 4, a new analytical geomorphic karst flow model is developed and compared with previous models in Chapter 5. A new hypothetical karst aquifer is established using MODFLOWDCM and the new geomorphic karst flow model is applied to predict simulated spring dischar ges from a new hypothetical karst system. In Chapter 6, the analytical karst flow models are tested using monitoring data from Blue Springs in Volusia County, Florida. Time series analysis was performed, to analyze observations of rainfall and groundwa ter levels near Blue Springs and the PAGE 20 20 spring discharge itself, in an effort to invest igate the hydrodynamic c haracteristics of the aquifer near Blue Spring. Chapter 7 presents conclusions and recommendations. PAGE 21 21 CHAPTER 2 LITERATURE REVIEW A karst aquifer is a complex hydr ogeologic system composed of a complex network of solution channels and conduits that is for the most part indeterminate. Four types conceptual models are typically used to characterize and simulate groundwater flow in karst aquifers: lum pedparameter models, contin uum models, discrete models, and dual porosity models (M. Ra hnemaei et al. 2005; S. L. Painter, A. Sun, and R. T. Green, 2004; Ford and Williams, 1989 Continuum models assume groundwater flow in karst aquifers is similar to flow in unconsolidated porous media (B. R. Scanlon, et al). For these model s, the grid cells representing a conduit are assigned high transmi ssivity and low storage. This approach is how MODFLOW is used to represent ka rst conduit flow. Discrete models represent karst conduits as pipes embedded in a diffuse flow system. Discrete models are intrinsically linked to finiteelement numeric al models (e.g. FEFLOW). In dual porosity model, two interacting flow systems representing conduits and diffuse flow system are modeled. Painter, Sun, and Green, 2004 conducted a series of numerical experiments to identify the limitati ons and advantages of the conti nuum models and dual porosity models to simulate groundwater flow in karst aquifer. Results suggested continuum models were flexible for a wide range of conditions; however, continuum models become inaccurate if the dens ity of conduits becomes large. In comparison, the dual porosity model provides more modeling flexibility. Lumped parameter models include simple water balance and timeseries models; these models are not physically based or emulate the physics of flow (e.g., N. Massei et PAGE 22 22 al., 2006; Samani, 2001; Panagopoulos et al ., 2006; M. Larocque et al., 1997; D.Labat et al., 2000; Rahnemaei et al., 2005; L. R. Ve stena, M. Kobiyama., 2007 .). Time series modeling usually treat a rainfall as an i nput and spring discharge and/or piezometric head as outputs to generate an analysis of the hydrodynamics of karst aquifers. These models yield estimates of lag times or aqui fer response times to rainfall events from which analyses typically show more than one flow component exists in a spring discharge series. This finding is consistent with several conceptual karst flow models proposed by Worthington (2007) [see Figure 21]. Worthingtons ( 2007) investigations of the hydraulic residence time distributi on support the general assumption that conduit and matrix flow coexist; but that the majority of gr oundwater is through the porous matrix (Figure 21). In general, the response function is not easy to interpret in the context of an unknown conduit network; however, there is some ev idence that previous research on the hydrologic response of surface drainage basins (Ignacio, et al. 1979) may provide a simple geomorphic analog for interpreting spri ng flows from karst aquifers. This concept is explored further in this dissertation to develop simple analytical, lumpedparameter models of karst aquifers. 2.1 Application of Time Series Techniques in Karst Anal ysis Several investigators (e.g., N. Massei et al., 2006; Samani, 2001; Panagopoulos et al., 2006; M. Larocque et al., 1997; D.Labat et al., 2000; Rahnemaei et al., 2005) have applied time series analysis to interpret trans ient flows through karst aquifers. Rainfall is treated as an input and sp ring discharge and/or piezomet ric head as outputs. The elements of the methods of time series analysis including correlation and spectral PAGE 23 23 analyses as applied to interpret karst aquife r flow are described in the following paragraphs. Figure 21. Twobox models for carbonate aquifers, with fraction of storage (in percentage) and residence time in each box and flow (m3s1) between boxes (Worthington 2007) Each time series is correlated to itself for increasing lag times. The autocorrelation function of a time series represents the linear dependence of an event on the remainder of the time series, the function characteri zes the duration that an event then influence subsequent data (M.Larocque, A.Mangin, 19 97). Thus, the autocorrelation function maps the effective memory in a time seri es (N. Massei, J.P. Dupont, 2006) and is defined as (Chatfield 1982): () () (0) ck rk c (21) PAGE 24 24 in which 0,1,...1 kN represents an integer number of ti me increments in the lag; (0) cis the variance or autocovariance with zero lag;() ckis autocovariance as a function of lag k and is calculated using following function (Chatfield 1982) 1 1()()()nk tkt tckNxxxx (22) where, N = number of observations, t x is the observation at timet, and x is the mean of the observations. Larocque, et al. (1997) anal yzed the autocorrelation of rainfall and spring discharge and showed the absence of mem ory in the rainfall series; whereas, autocorrelations for spring discharge tended to diminished very slowly when lag time increases. This suggested the aquifer, dr ained by the spring, had a large storage capacity. The cross correlation function measures the correlation between t x and tk y and is calculated from equation (23), (Chatfield 1982) : 00xy xy xxyyck rk cc (23) where,xyckis cross covariance function: 1 11 ,0,1,...,1 cov(,) 1 ,1,2,...,(1)Nk ttk t xyttk N ttk tkxxyykN N ckxy xxyykN N (24) Similarly, ()cov(,) x xt t kckxx; ()cov(,)yy ttkckyy in which x and y represent two different time series. That is, t x represents observations of series x at timet; x is the mean of the observations of series x ; t y represents the observations of series y at timet; and yis the mean of the obs ervations of series y. PAGE 25 25 The correlogram is a graphical repres entation of correlation coefficient ()xyrkas a function of lag k, where values of ()xyrkare plotted as ordinates against their respective values of k (Nozar 2001). A perfectly symmetrical crosscorrelation function centered at 0 k indicates both input and output signals react at the same time to a third independent signal (M.Laroc que, A.Mangin. 1997). The maximum values of ()xyrkobserved at a lag d indicates that one series is re lated to the other when delayed by time periodd. The lag time between rainfall and spring discharges or groundwater levels gives an estimation of the particle travel time s through the aquifer and also indicates the degree of karsification and type of flow regime. A short time response to rain of less than 10 15 days corresponds to high degree of karstification, and a long time response of more than 2 months, corresponds to low degree of karstifiction (Obarti et al. 1988). The power spectral density function is t he Fourier transform of the autocovariance function and is used to investigate the periodi c structure of a time series. Spectral density function() Sfis calculated as (Shumway and Stoffer 2006): 2()ifk kSfcke 1/21/2f (25) Where f is the frequency, expressed in cycles per point;() ckis the auto covariance function (22); and 1 i The spectral density is an even function of period one. Due to this evenness, the function of spectral dens ity can be expressed as: 102cos2kSfcckfk (26) PAGE 26 26 The cross spectrum function is the Fourie r transform of the cross covariance. The coherency function expresses the linearity of a system. And a system is linear when a change in input creates a proportional change of the output (Larocque et al.1998). This linearity is a characteristic of a highly karstified aquifer, where a heavy rainfall event produces a strong discharge of the karst spring over a short time. The coherence function ()xyChfcan be derived from the cross spectrum function and is defined as (Nozar 2001): 22()() () ()()xy xy xy xyCfqf Chf SfSf (27) Here ()xSfand ()ySfare spectral densities of time series of t x andt y x yCis the real part of the crossspectrum called the cospectrum, is given by (Nozar 2001): 11 ()(0)cosxy xy xy yx kCfcckckfk (28) x yqis the complex part of the cross spectr um, called the quadrature spectrum and is given by (Nozar 2001): 11 ()sinxy xy yx kqfckckfk (29) where x ycand yxcare the cross covariance functions (24). The value of coherence spectrum varies between 0 and 1, where 1 indicates that the two series are fully dependent and for zero means that the two series are independent (N. Samani 2001). The phase spectrum ()xyPf determines the time lag between the two series and defined as (Nozar 2001) PAGE 27 27 () ()arctan ()xy xy xyqf Pf Cf (210) Where x yqis quadrature spectrum (29), x yCis the cospectrum (28). The plot of ()xyPfagainst frequency f is the phase diagram. The phase between the two series is obtained from the phase diagram at the frequency at which coherence is maximum. In turn, division of phase by angular frequency gives the time lag between the tw o series. 2 T k (211) Where, kis the time lag, is the angular frequency, and Tis period. Samani (2001) estimated the period of spring discharge from a karst aquifer us ing a plot of the spectral density, and lag time from equation (211), based on plots of coherency and phase function. 2.2 Concepts of Geomor phologic Structure Model Time series analysis can assist with esti mation of aquifer response times (lag times) to rainfall events; but, the physical ba sis of the transient response function of karst basin requires further analysis. I. R odriguez, et al. (1979) conducted research on surface watersheds that suggests stream flow response at basin outlet can be expressed in terms of geomor phologic characteristics of the basin depicted in Figure 22. They developed a geomorphologic structur e model to link the instantaneous unit hydrograph (IUH) of a surfac e basin, which represents the hydrologic response of the surface basin, to Hortons numbers. The geomorphic modeling approach may have utility in the simulation of spring flows from karst aquifers. PAGE 28 28 Figure 22. Third order surf ace basin with ordering system (I. Rodriguez, et al. (1979)) Karst conduits, as suggested by Worthi ngton (2007), are analogous to roofedover creeks. In general, there are four type s of karst conduits pa tterns (Figure 23) (N.Palmer 1991). Branchwork conduits consist of passages that join as tributaries (Figure 23a). Network conduits are intersecting fissures formed by the widening of nearly all major fractures within favorable areas of soluble rock (Figure 23b). Anastomotic conduits consist of curvilinear t ubes that intersect in a braided pattern with many closed loops (Figure 23c). Ramiform conduits consist of irregular rooms and galleries wandering three dimensionally with branches extending outwa rd from the main areas of development (Figure 23d). The conduit pattern of the first type is very similar to a surface basin shown in Figure 22. The remaining three conduit network configurations appear different the first; however, the flows in these netwo rks can be recreated using smaller and less PAGE 29 29 bifurcated networks of the first type. In this way, the tributary conduits in a karst aquifer are analogous to a surface basin system. In this dissertation, several analytical geomorphologic structure models are develo ped as new simulation tools to predict spring discharge from karst aquifers. The theor etical basis of these geomorphic karst models is rooted in the geomorphologic stru cture model for surface basins developed by Ignacio, et al. (1979). The next secti on presents the fundament als of analytical geomorphologic modeling. Figure 23 shows a hypothetical waters hed and ordered streams such that: (a) channels that originate at a source are defined as firstorder st reams; (b) when two streams of order join, a stream of order 1 is created; and (c) when two streams of different orders join, the channel segment i mmediately downstream has the higher order of the two combining streams. The expressions of Hortons laws are (Schumm, 1956): Hortons law of stream numbers is expressed as: 1/BNNR B R is the bifurcation ratio, whose value is normal by between 3 and 5. Hortons law of stream lengths is defined as: 1/LLLR L R corresponds to the length ratio, whose value is no rmally between 1.5 and 3.5. Hortons law of stream areas equates to: 1/ A AAR A R is the area ratio, whose value is normally between 3 and 6. PAGE 30 30 Figure 23. Common patterns of conduits (N.Palmer 1991) To evaluate the instantaneous unit hydrogr aph (IUH) for a watershed such as the third order drainage basin depicted in Figure 22, it is assu med the rainfall input is a unit volume of effective precipitation uniformly distributed over the basin and instantaneously delivered to the land surface. The cumulat ed water volume yielded at the outlet as output up to a certain time t can be interpreted in terms of probabilities (i.e., what is the probability that a drop of rainfa ll chosen at random arrives at the basin outlet at time t ). The process which describes the fate of a rainfall drop flowing through the basin depends on the residence time in each stream reach (i.e., ordered stream segment) as well as the number of reaches traversed to arrive at a specific reach at time t In a basin, the residence time in each stream reac h depends on the location of the rain drop, because different reaches in the same basin possess different hydraulic characteristics. PAGE 31 31 Thus, the fate of a rain drop, as recogniz ed by Ignacio, et al. (1979), is considered a semi Markovian process in which a drop of rainfall occupies successive stream segments or reaches (or st ream orders) described by transition probabilities of a Markov process and the residence time in each segment is described by a random variable that depends on the specific reach cu rrently occupied and on the next reach to which the next trans ition will be made. The probability that a drop of rainfall c hosen at random has reached a given stream order or segment at (or before) timet, was defined by Ignacio et al. (1979) using the follo wing equation: ()(0)() tt (212) where () t is the flow segment probability matrix, whose elements ()it give the probabilities that drops occupy flow segment iat timet; (0) represents the initial probability vector, which depends on the spatial character of the rainfall and whose elements0igive the probabilities that the process starts at flow segment i, or, in other words, that the drop starts its travel in streams of orderi; and ()ijt is the interval transition probabilities of a continuous ti me semiMarkov process Howard (1971). Creating an expression of the interval transition probabilities ij for a general model of continuous time semiMarkov process is quite complicated (Howard 1971). Ignacio, et al. (1979) simplified the problem assu ming first, the hydraulic residence time distribution of a flow segment was independent of destina tion segment, and second, the residence time distribution iw was exponential for a rain drop in stream of order i; thus, PAGE 32 32 ()iiiwe (213) The second assumption is equivalent to assuming each stream reach behaves like a linear reservoir. The transition probability matrix defines the probabilities that the drop in one stream reach transition to another. The matrix is generally defined as: 1213 1 111213 1 23 2 212223 2 3 1230 ...0 ... 00 0 0000 ... ...... ... 000...01N N NNNNNppp pppp pp pppp P p pppp (214) where, is the order of the basin and N is the total number of model segments or 1 N. The surface basin shown in Figure 22 can be represented by the flow system shown in Figure 24. Different orders str eams are represented by different flow model segments (ellipses and the circle shown in Figure 24). Each elem ent of the transition matrix P expresses the probability a rain drop ma kes a transition from flow segment of orderito flow segment of order j Thus, ij p represents the proportion of drops in flow segment order i, that move to flow segment order j The highest order stream of the geomorphic model is represented using two model segments to insure the basin produces a discharge hydrograph at the outlet which increases continuously from a zero flow rate at t=0 (Ignacio et al 1979). For the third order basin depicted in Figure 22, secondorder streams can on ly go to the third order stream, but now the thirdor der stream is represented by model segments 3 and 4. Segment 3 receives drops from all secondorder streams and some firstorder streams. PAGE 33 33 Drops received by segment 3 are passed to s egment 4, which feeds the outlet of the basin (the fifth segment). Together segments 3 and 4 have a mean hydraulic residence time of1 3 Assuming the hydraulic residence time follows an exponential distribution, 3w the following expression is relevant: 2 3we (215) where the mean value of 11 32 3 4 5 Figure 24. Representation of a third order basin as a c ontinuous Markov process (I. Rodriguez, et al. (1979)) In the process of using two model segment s to simulate the third order stream reach and another model segment to repr esent the basin outlet, the transition probability matrix is now given by 1213000 00100 00010 00001 00001 pp P PAGE 34 34 To obtain the basin dischar ge at the outlet (model segment N=5), an expression is needed for the probability that a drop chosen at random has reached the outlet at (or before) time t. Thus, t t t t35 3 25 2 15 1 50 0 0 (216) And for uniform precipitation, *1 10TA A *1 20TA A *1 30TA A (217) where *iA1,2,3 i represents the area draini ng directly into all of the streams of order iand TAis the total area of the basin. From equation 216, t he IUH is defined: 51 52 53 5 123000 dtdtdtdt IUHt dt dt dt dt (218) Ignacio, et al. (1979) devel oped equations which express terms in the equation (218) as functions of geomor phologic Hortons numbers: 2 12 222 2BB BBRR p RR (219) 2 13 232 2BB BBRR p RR (220) 22 10BA R R (221) 32 2 222 0 21BBBB AAB R RRR RRR (222) 2 3 232 1 01 21BBB B AABRRR R RRR (223) The derivative terms in equation (222) are given by (Ignacio, et al. 1979): PAGE 35 35 ** 33 1215 1 2 3 41tt tttAeAeAteAe (224) 33 22 32 3 2 32 25 22 32 32 322 ()1tt ttee te (225) 3335 3()1tttete (226) where, 2 31132 1 2 2131p A 2 3112 2 2 2132p A ** 3121313 3 ** 1332 p A 322 2 ******** 43113312332313121231313 222 *** 31323322 App in which 332 andi is the inverse of the mean wait ing time in segment of order i. Assuming v is the average stream flow ve locity in the basin, then: /iivL (227) Where, iLrepresents the mean length of the flow segment of orderi. Equation (227) implies 11/vL 1 21L R 2 31L R (228) Substituting equations from (219) to (228) into equation (218), IUH can be expressed as a function of Hortons numbers and the mean flow velocity in the basin. PAGE 36 36 With the equation (218), the structure of the hy drologic response is linked to the geomorphologic parameters of the surface water basin. Due to dual porosity of karst aquifers, more than one lag time is typically observed, which is related to flow of water th rough karst conduits and matrix of the rock (Rahnemaei et al., 2005). When using the geomorphologic model to analyze the response of a karst aquifer to rainfall, it may be necessary to simulate flow subdivided into two or more components. 2.3 Numerical Models of Groundwat er Flow in a Karst Aquifer DCM, a dual conductivity model, developed by Southwest Research Institute (SwRI), is a new module in MODFLOW. Painte r, et al. (2004) succe ssfully applied DCM to the Barton Springs Segment of the Edwa rds Aquifer. Results showed DCM was able to simulate spring flows and it offered greater flexibility in simulati ng flows as compared to using a continuous matrix fl ow model such as MODFLOW. DCM simulates two coupled flow system s, one representing conduits and the other diffusive flow. Both diffuse and conduit systems are modeled as coupled continuums; whereby, fluid exchange bet ween systems is assumed linear. Conduit flows are assumed turbulent in DCM; thus, Da rcys law is not valid. Instead, the DarcyWeisbach equation is used to expr ess frictional energy losses: 22HLv hff Dg (229) Where h represents the energy head lost. L is the length of the conduit, H Dis the mean hydraulic diameter, vis the mean velocity in the conduit, gis the acceleration due PAGE 37 37 to gravity, and f fis the friction factor. For straight conduits, f fdepends on the relative roughness of the conduit and on the Reynolds numbere R Natural conduits of karst aquifers are not straight and they do not have a constant diameter; therefore, they suffer head losses due to bends, ex pansions, and contractions in addition to frictional losses. Thus, a mo re comprehensive expression for energy losses is equation 230: 222222bends ec HLvvv hffCC Dggg (230) where bendsCis an empirically coefficient acc ounting for head loss due to bends, and ecCis a coefficient for crosssectional exp ansions/contractions. Each term in equation (230) uses same power expression of ve locity. Hence, the e ffects of the bends and crosssectional variations can be grouped into an effective friction factor e f f. Hence, equation (230) can be rewritten as follo ws: 22e HLv hff Dg (231) For conduit flow, the following analogous equation for specific discharge is used in DCM: ck qh h (232) where, qis the Darcy velocity, and ckis an effective conductivity for the conduit. The governing groundwater flow equat ions in DCM are therefore: ,()mmm mmxmmymmcmcm mhhh ShThThQhhhh txxyy (233a) PAGE 38 38 ,()ccc ccxccycccmcm chhh ShThThQhhhh txxyy (233b) where the subscript c denotes the conduit system and m denotes the background porous matrix flow system; his hydraulic head; x yTTis transmissivity in the x andy directions; Sis a storage term; and Qis the volumetric source term per unit area of the aquifer. The final term in equation (233) r epresents the exchange of fluid between two systems. is the exchange coefficient. The e xchange is linear and depends on the hydraulic heads of conduit syst em and background system. Thus, 0min,max,,top botbot cmccc topbot ccZhhZZ ZZ (234) where, top Z denotes the elevation of the top of the aquifer and bot Z the bottom elevation. To model turbulent flow in conduits, equation (233) can be replaced with the following equation: ,()mmm mmxmmymmcmcm mhhh ShThThQhhhh txxyy (235a) ,()c ccxcyccmcm ch ShqqQhhhh txy (235b) where, 1/2 1 cc cxcxchh qTh x x (236a) In the situation of turbulent flow, and x h hTqc ccx cx (236b) PAGE 39 39 and the parameter 1cxT is a proportionality constant inco rporating geometrical properties of the conduit and friction factor. The DCM package is currently limited to simulating flow in singlelayer aquifers; however, DCM package uses two MODFLOWlike layers to create a singlelayer karst flow model. Conduits in one MODFLOW layer si mulate diffusive matr ix flow and in the other conduit flow. Within the conduit layer, ac tive cells are those cells intersected by the conduit. The successful application of DCM to the Barton Springs Segment of the Edwards Aquifer served to demonstrate DCM as a credible karst aquifer simulation tool (Scott Painter, Sun, and Green, 2004). In this research, the MODFLOWDCM is used to create spring discharge data from hypothetical karst flow aquifer wit h known aquifer parameters and conduit networks. The virtues and ut ility of different analytical geomorphic models are then explored using this hypothet ical spring discharge data. PAGE 40 40 CHAPTER 3 DERIVATION OF GEOMORPHOL OGIC KARST FLOW M ODELS The technical definition of a karst aquifer is a body of soluble rock that conducts water principally by conduits formed by t he dissolution of the rock. The aquifer is commonly structured as a branc hed network of conduits, which together convey flow to a spring (Worthington, 2002). Here, we will assume these conduits are analogous to a network of stream segments t hat converge to discharge basin flows through a surface spring. The basic idea is to follow Ignacio, et al. (1979) and develop simple karst flow models based on Hortons geomorphic numbers for streams but now interpreted as Hortons numbers for a karst network of conduits. Thus,B R is the law of conduit numbers (the bifurcation ratio); L R is the conduit length ratio; and ARrepresents the law of conduit recharge areas. Properly developed, these simple kars t flow models could be used estimate effective values ofB R L R andAR for an indeterminate karst conduit network using monitored rainfall inputs and spring discharges. The optimal values ofB R L R andAR could then be used to design equivalent cond uit networks in more complex numerical flow models that would then describe transi ent spring discharges. Following the surface watershed work of Ignacio, et al. (1979), th ree new geomorphologic structure models for karst aquifers are developed in this chapter. 3.1 Derivation of a Third Or der Geomorphologic Karst Model The karst conduit network is essentially an indeterminate system, because it is not visible for inspection. Thus, Hortons numbe r for conduits cannot be determined in the same way as Hortons numbers for a surface basin (by direct inspection). Instead, PAGE 41 41 preliminary assumptions and specificat ions will be necessary before model development can begin. In general karst systems show evidence of both laminar flow in the porous matrix and turbulent flow through conduits. The porous medium is a ssumed here to be the first order flow segment, and it receives all ra infall. The conduit network is composed of second and third order flow segments. The outlet of a kars t aquifer, which is normally a karst spring, is referred to as the fourth order flow segment (Fi gure 31). Using two linear reservoirs to represent the third order flow segment (this ensures the IUH begins at zero at t=0) and yet another model segment for the sp ring, the initial third order flow system is essentially modeled as the fifth order system illustrated in figure 32. Figure 31. Considering a karst aquifer as a third order geomor phologic flow system PAGE 42 42 Figure 32. Conceptual model of a third order karst aquifer using five model flow segments identified in par entheses or as standalone numbers (e.g., 1,2, 3,4,5). We assume the karst aquifer receives ra infall through the porous medium at the ground surface or the first order flow segment in Figure 32. Thus, the elements of the initial probability row vector iiin equation (217) are: 101 23000 (31) With these initial probability values, equation (218) reduces to: 15dt IUH dt (32) Substituting equation (224) into equation (3 2), yields an expression for the IUH of the karst system (Figure 31) equation 33. *** 333 12** 112233343 ttt tt I UHAeAeAeAteAte (33) where, 2 31132 1 2 2131p A PAGE 43 43 2 3112 2 2 2132p A ** 3121313 3 ** 1332 p A 322 2 ******** 43113312332313121231313 222 *** 31323322Ap p 2 12 222 2BB BBRR p RR 2 13 232 2BB BBRR p RR 11/vL 1 21L R 2 31 L R 332 Assuming rainfall as the primary source of aquifer recharge, karst spring discharge (the output of a karst aquife r) is calculated using equation 33 and rainfall data in the following formulation: 1()()(1)t iQtPiUti (34) where, Q is spring discharge, p is rainfall, Uis the instantaneous unit hydrograph of the karst basin, tis the incremental time step, 1,2,...i PAGE 44 44 3.2 Derivation of a Fourth Or der Geomorphologic Karst Model It is likely, depending on the level of karsif ication, that an aquifer could be best represented as a fourth order flo w system as shown in Figure 33. Again, the intact porous matrix could be the first segment followed by a third order co nduit network. In its entirety, the aquifer would be a fourth or der system. The conduit network would be composed of second, third and f ourth order flow segments. Using two linear reservoirs to represent the fourth order conduit flow segment (this ensures the IUH begins at zero at t=0) and yet another model segment for the sp ring, the initial f ourth order flow system is essentially modeled as the sixth or der system illustrated in Figure 33. Figure 33. Conceptual model of a fourth order karst aquifer using six model flow segments identified in parentheses or as standalone numbers (i.e., 1,2,3,4,5,6). For the flow shown in Figure 33, the probability that a drop of rain chosen at random has reached the outlet, which is referred as the 6th order flow segment, at (or before) time t is given by: PAGE 45 45 61162263364460()0()0()0()ttttt (35) Where, 6t represents the probability that the drop of rain is found in segment i= 6 at transition stept. 6(),(1,2,3,4)iti is the probability that the process goes from segment ito 6 at transition stept. 0,(1,2,3,4)ii is the probability that the rain drop is initially located in segment i. The initial conditions for a flow system shown in Figure 33 is 101, 2340000 (36) Hence, the probability that a drop chosen at random has reached the outlet at (or before) time t can be expressed as: 6116160()()ttt (37) And, the instantaneous unit hydrograph (IUH) of the a quifer can be calculated from: 61 6dtdt IUH dtdt (38) The next step in the derivation of the fl ow model is to express equation (38) in terms of conduit Hortons numbers. For a flow system shown in Figure 33, the transition probability matrix is now given by: P=121314 232400 0 0000 000100 000010 000001 000001ppp pp (39) The mean waiting time matrix is 1 where PAGE 46 46 1 2 3 4 400000 00000 00000 000200 000020 000000 (310) where i is the inverse of the mean re sidence time in systems of orderi. By defining a transition rate matrix as: A PI= A te11 2 1 3 1 4 22 3 2 4 3 4 4000001 00 000000100 00000001100 000200000110 000020000011 000000000000ppp pp 1112113114 1112113114 2223224 2223224 33 33 ** 44 ** 4400 00 00 0 00 0 00 0000 00 000220000 0 0000220000 000000000000ppp ppp pp pp (311) where, 42 ; the interval transition probability can be written as [Howard, 1971; and Ignacio et al. 1979]: A tte (312) in which A teis defined as 22(/2!)...IAtAt Howard [1971] showed t hat the exponential transform of (312) is given by: PAGE 47 47 1etsIA (313) where, 1112113114 2223224 33 ** **00 000 0000 0000 0000 00000sppp spp s sIA s s s Since, we are only interested in the 6th column of()t only the 6th column of 1sIA is needed: 2 *2 1431421431321224 1 23132312232314231224 *2 2124323324 *2 312 2 *2 123 123 ** 123 2 1231spspspspspp pppppp ssppp ss sssss sss ssss ssss (314) The term 1sIA can be expressed in a partial fraction expansion: 1 *2* 123111111 ()()ij ij ij ij ij ijsIAabcd e f ssssss (315) Equation (315) is the expression of()et, and the interval transition probability matrix is obtained by inverse ex ponential transformation. Hence, PAGE 48 48 ** 3 12t tt tt ijijijijijijtaebecedteeef (316) and the sixth column of ()t can be expressed as: ** 3 12161616 16 16 16 16()t tt tttabecedeetefe (317) From allowing equation (3 14) to equal the 6th column of equation (315), the parametersija ijb ijc ijd ije ij f can be calculated as: 1612*241312*23141apppppp (318a) 2 23232323 2 16 *1312*1 2 121213131122413122314 /// 122414141314pppppp b pppppp = 22 23121213131* 2 1312*1122414141314 pppppp (318b) 21 2 2*32 2 *1 2 3 3 1612242324 pppp c (318c) 2 1*223 16 2 *3231312231313pppp d (318d) 2 **2*2*3*3 161* *3*2*1 32323232141412241314 // 122413122314 pppppp e pppppp (318e) PAGE 49 49 32 232 *2 3*13*2 22 2 *23 123* 22 2 132 123* 42 2 2*1*23 1*23 22 1*23 1*23 2 123* 2 1621421421223 31223214 1223213 41441224413 31223414 2122361224 pppp ppp ppp pppp ppp pppp f 3 *3 32222 2*3132231 2243 1233*1*2 22322 13*13*13* 3234 3*13**1 42222 *3*32*23 222 *23*814 1314 1224313214 1321314 21421314 414413414 4143 pp p pppp ppp ppp ppp p 22 2 23*21 22 222 2 *32 *23*2 3 22 22 *32 *2 1 22 22 *2 3* 23 22 *231*132 2222 *1 32*32 2 *121314 21321431224 312231224 4122421224 214213 2122421223 212 pp pppp pppp pppp pp pppp pp 222 1*3*2*1 53 3*12* 233 2*23* 34 23*2*/// 2421421224 21224613 4122331224 ppp ppp pppp (318f) 622423 app 620 b 2 233* 62 2 *232242324 ppp c 2 2* 62 2 *33223 p d *332* 62 *2*3242324 ppp e PAGE 50 50 2222 3*3*2323 22 22 623*3*23*23*2*3*2 23 2**2242232324 323424224223// 24224 pppp fpppp pp 631 a 630 b 630 c 2 63 2 *3d 3* 63 *3e 3*3 63 2 *32 f 641 a ; 640 b ; 640 c ; 640 d ; 64*e ; 641 f 651 a ; 650 b ; 650 c ; 650 d ; 650 e ; 651 f 661 a ; 660 b ; 660 c ; 660 d ; 660 e ; 660 f The transition probability ij p expresses the probability that a drop goes from a segment of order i to a segment of order j. ij p can be calculated using the following equation [IGNCIO, JUAN. 1979]. i ijorderofconduitsofnumber total jordertodrainingiorderofconduitsofnumber p (319) 1,2;2,3,4ij There are N1 conduits of order 1 of whic h 2N2 combine to make conduits of order 2. The remaining (N12N2) chan nels of order 1 drain into karst pipes of orders 2 and 3. PAGE 51 51 Following Smart [1968], it is assumed that the lengths of interior conduits in a given network are independent random variables drawn from a common population. This assumption implies that the distribution of interior conduit lengths is independent of order, magnitude, or any other topologic char acteristic, from which we may write that the (N12N2) karst pipes of order 1 join conduits of order 2, 3 and 4 according to 1 1 )221( orderofconduitsofnumber total jordertodraining orderofconduitsofnumber NN (320) 2,3,4 j The mean number of karst conduits of order in a finite network of order is [Smart, 1971] is given by: 1 21 ,,2,3,... 21a a aN EvN N (321) Thus, the number of conduits of order 2 equals 21 21 21 NN N ; the number of karst pipes of order 3 equals 12 3 2311 2121 NN N NN ; and the number of condui ts of order 4 equals to 12 43 2311 1 2121 NN NN NN (where, 41 N ). On the average the number of conduits of order 1 that drain into secondorder conduits is 21 2 212 21 1212 343 223231 21 22 1 1111 1 2121212121 NN N NNN NN NNNN NNN NNNNN (322) Then, PAGE 52 52 21 2 212 21 1212 343 22323 12 11 21 22 1 1111 1 2121212121 NN N NNN NN NNNN NNN NNNNN p N 2 3222 2 4221 RBRBRB R BRBRBRB (323a) Similarly, 212 3 23 21 21 1212 343 22323 13 111 2121 22 1 1111 1 2121212121 NN N NN NNN NN NNNN NNN NNNNN p N 3 322 1 4221 RBRB RBRBRBRB (323b) 212 43 23 21 21 1212 343 22323 14 111 1 2121 22 1 1111 1 2121212121 NN NN NN NNN NN NNNN NNN NNNNN p N 32 3221 1 3221 RBRBRB RBRBRBRB (323c) Finally, there are N2 conduits of order 2 of which 2N3 make up for the conduits of order 3. The remaining (N22N3) karst pipes of order 2 drain into conduits of orders 3 and the (N22N3) pipes of order 2 join ka rst pipes of order 3 and 4 according to 2 2 )322( orderofconduitsofnumber total jordertodraining orderofconduitsofnumber NN (324) 3,4 j PAGE 53 53 Hence, 23 p and 24 p can be calculated as 12 3 23 323 1212 343 2323 23 211 2121 22 1111 1 21212121 NN N NN NNN NNNN NNN NNNN p N 22 21 RB R BRB (325a) 12 43 23 23 1212 343 2323 24 211 1 2121 2 1111 1 21212121 NN NN NN NN NNNN NNN NNNN p N = 2 232 2 RBRB R BRB (325b) From equations (38) and (316), the in stantaneous unit hydrogr aph (IUH) for the fourth order flow system shown in Figure 33 can be calculated from: 3 12***16 116 216 316 *16 16 *t tttttIUHbecedeeteeefe (326) After substituting the expressions of ij p and expressions of ,,,,,ijijijijijijabcdef to equation (326), the IUH (326) for the kars t aquifer is expressed in terms of hortons conduit numbers and residence time i 3.3 Derivation of The Fifth Or de r Geomorphologic Karst Model If a karst aquifer is believed to be best conceptualized as a fifth order flow system shown is Figure 34, equations for a fi fth order geomorphologic karst model can be PAGE 54 54 developed using same techniques previous ly applied to derive the above third and fourth order geomorphologic karst models. Figure 34. Conceptual model of a fifth or der geomorphologic Kars t system using seven model flow segments identified in par entheses or as standalone numbers (i.e., 1,2,3,4,5,6,7). The probability that a rain drop chosen at random has reached the outlet (the 7th order system) at (or before) time t is given by: 71172273374475570()0()0()0()0() tttttt (327) Where, 7t gives the probability that the drop is found in segment i=7 at transition step t. 7(),(1,...,5)iti represents probabilities that the process goes from segment ito 7 at transition stept. 0,(1,...,5) ii represent probabilities that the process starts at segmenti, or, in other words, that the drop starts its travel in flow segments of orderi. PAGE 55 55 Only the first order segment, which repr esents the upper part of a karst aquifer, receives rainfall; therefore, the initial cond ition for fifth order geom orphologic karst flow system is 101, 234500000 (328) The probability that a drop chosen at r andom has reached the outlet (the seventh order segment in the system) can be expressed as: 7117170()() ttt (329) And, then the instantaneous unit hydrograph (IUH) for the flow system shown in Figure 34 can be calculated as: 71 7dtdt IUH dtdt (330) The next step of this derivation is to express 17tas a function of geomorphologic Hortons numbers. For the fl ow system shown in Figure 34, the transition probability matrix is: P=12131415 232425 343500 0 00 00 00000 0000100 0000010 0000001 0000001 pppp ppp pp (331) The mean waiting time matrix is 1 in which PAGE 56 56 1 2 3 4 5 5000000 000000 000000 000000 0000200 0000020 0000000 (332) where i 1,2,3,4,5 i are the inverse of the mean wait ing time in flow segments of orderi. The transition rate matrix is then defined as A PI=Ate11 2 1 3 1 4 1 5 22 3 2 4 2 5 33 4 3 5 4 5 50000001 00 0000000100 00000000100 0000000001100 00002000000110 00000200000011 00000000000000 pppp ppp pp 1112113114115 1112113114115 2223224225 2223224225 3334335 3334335 44 44 55 5500 00 00 0 00 0 00 0000 00 000 00000 0 0000220 0000022 0000000 pppp pppp ppp ppp pp pp ** **0 00000 00000 0000000 (333) where, 52 The interval transition probability can be written after [Howard, 1971; IGNACIO, et al 1979]: PAGE 57 57 Atte (334) where Ateis defined as 22(/2!)... IAtAt The exponential trans form of (334) is given by: 1etsIA (335) where, 1112113114115 2223224225 3334335 44 ** **00 000 0000 00000 00000 00000 0000000 spppp sppp spp sIA s s s (336) The term, 1 s IAcan be expressed in a partial fraction expansion form: 1 *2* 12341111111 ()()ij ij ij ij ij ij ijsIAabcde f g sssssss (337) Hence, 2 12341111111e ijijijijijijijtabcdefg ssssss s (338) Using the method of inverse exponential transformation, interval transition probability matrix can be expressed as: ** 3 124t ttt tt ijijijijijijijtaebecedeetefeg (339) The seventh column of () t is: PAGE 58 58 3 124**17171717 1717 17 17()t ttttttabecedeeeftege (340) To calculate 17() t using equation (340), parameters ,,,,,ijijijijijijabcdef must be determined. By equating equations 337 and 338, the parameters ,,,,,ijijijijijijabcdef can be calculated (See the appendix). The final step of this derivation is to express the transition probability ij p as functions of geomorphologic Hortons numbers Similar to the development for the fourth order geomorphologic karst flow system, the transit ion probability ij p for the fifth order geomorphologic karst system is: 6543 12 65432422 22 (1) 8442221 RBRBRBRB p R BRBRBRBRBRBRBRB (341a) 6532 13 6543222 2 (1) 8442221 RBRBRBRB p R BRBRBRBRBRBRB (341b) 643 14 654322 (1) 8442221 RBRBRBRB p R BRBRBRBRBRBRB (341c) 6542 15 6543212 (1) 8442221 RBRBRBRBRB p R BRBRBRBRBRBRB (341d) 32 23 3222 2 (1) 4221 RBRB p R BRBRBRBRB (341e) 3 24 322 (1) 4221 RBRB p R BRBRBRB (341f) 32 25 3212 (1) 4221 RBRBRB p R BRBRBRB (341g) 3422 (1) 21 RB p R BRBRB (341h) PAGE 59 59 3512 (1) 21 RB p R BRB (341i) Substituting the expressions of ij p (338) and expressions of ,,,,,ijijijijijijabcdef to (340), 17 can be calculated. From equation (330), the following formulation for IUH is derived: 3 124***17 117 217 317 417 17 *17 *()()()() ()()t ttttttIUHbecedeeefeftege (342) The spring discharge for a fifth order kars t aquifer can be simulated using equation (342). PAGE 60 62 CHAPTER 4 APPLICATION OF GEOMORPHOLOGIC KARST MODEL In chapter 3, governing equations for thir d, fourth, and fifth order geomorphologic karst flow models wer e developed. The focus of this chapter is to apply these models to simulate the spring discharge from a hypot hetical karst aquifer and then assess their performance. As noted above, the conduit net work is in general indeterminate for most karst aquifers, i.e., there is limited informa tion on conduit locations and geometries. This becomes an issue with the application of geomorphologic models to karst aquifers; because, such models are expressed in terms of Hortons numbersB R and L R which cannot be assessed by direct measurement as is possible for surface basins. To overcome this issue, field observations (e.g., spring discharge and heads) must be used to calibrate these parameters. In this chapter, calibration/validation studies are conducted with the geomorphologic models. Synthetic rainfall and spring discharge data are created and used to conduct these studi es. The synthetic spring discharge data are generated from simulating a 2dim ensional hypothetical karst aquifer with known characteristics (e.g., conduit network geometry and matrix permeability) using MODFLOW2000 DCM. 4.1 Application of DCM to Construct a Hypothetical Karst Aquifer This section describes the used of the DCM package is to create synthetic spring discharge data set for a hypothetical single layered karst aquifer. The karst aquifer possesses a known conduit network that discharges to a spring. A synthetic rainfall data set is the only source of aquifer re charge and the aquifer parameters and conduit network are loosely based on the data extracted from the BSEACD (Barton Springs/Edwards Aquifer Conservation Distric t) model (Painter, Sun, and Green, 2007). PAGE 61 63 The conduit network of the hypothetical ka rst aquifer is illustrated in Figure 41. Specifics on the DCM model of the hypothetical aquifer include: grid cell dimensions of 1,000 ft long and 500 ft wide; noflow boundaries that encompass the entire aquifer; a spring that is modeled as a drain at t he end of the conduit network; and a spatially uniform, time varying, synthetic rainfall time series that is the only source of recharge over a total area of 910*141.2square feet. Figure 41. Conduits network in the hypothetical karst aquifer Hydraulic conductivities for each conduit and storage parameters in both conduits and diffuse layer are the same those in t he BSEACD model. The hydraulic conductivity values are shown in Table 41. The exchange parameter 0 is 0.001 1day for all conduits. The hydraulic conductivity for diffus e layer ranges from 0.5 ft/day to 100ft/day. Storage parameters varied over two orders of magnitude. Values of specific storage ranged from 510in conduits to 5310 in the diffuse system. Specific yield varied from 310in conduits to a value 3310 in the diffuse systems. PAGE 62 64 Table 41. Hydraulic conductivities of conduits in the hypothetical karst model Conduit Number Hydraulic Conductivity (ft/day) A 160000 B 15000 C 10000 D 10000 E 5000 f 5000 G 4200 H 5000 I 10000 J 5000 Daily rainfall data for 5,286 days are shown in Figure 42. The hypothetical karst aquifer is simulated as transient model with 5,286 stress periods. The resulting spring discharge generated from the MODFLOW2000 DCM simulation is shown in Figure 43. Figure 42. Daily rainfall data used in the hypothetical karst model PAGE 63 65 Figure 43. Spring discharge calculated by the DCM 4.2 Cross Correlation Analysis of the Hypothetical Karst Flow Sy stem Cross correlation analysis was conducted to elucidate hydraulic response times from the DCM simulated hypothe tical karst aquifer. Figure 44 shows cross correlations of rainfall (Figure 42) versus simulated sp ring discharge (Figure 43) From Figure 44, it is evident the cross correlation function (CCF) peaks (nearly 0.7) for a lag time equal to zero. And for lag times greater than zero, values are relatively low. The notably high cross correlation at a lag time equal to zero suggests the existence of a flow component wh ich is linearly correlated to rainfall. Comparing rainfall data (Figure 42) and DCM gener ated spring discharges (Figure 43), major discharge fluctuations occur on the very day of most ra infall events. It is important to emphasize this observation is generated from a model and not from direct observations of spring flow. PAGE 64 66 Figure 44. Cross correlation of rainfall and spring discharge of the hypothetical Karst aquifer. To distinguish the nature of shortterm sp ring discharge fluctuations due to rainfall from spring discharge more representative of longterm recharge, it is necessary transform (smooth) the dischar ge series, or in this case use linear interpolation to estimate the running baseflow discharge fo r those days when rainfall was recorded. The transformed discharge series is refe rred, hereafter as the smoothed spring discharge. Shortterm dischar ge fluctuations were obtained by subtracting values of smoothed spring discharge series from spring discharge data originally calculated by DCM. The smoothed spring discharge series is plotted in Figure 45, while Figure 46 illustrates fluctuations. Figure 47 and Figure 48 present s respective plots of cross correlation of rainfall versus spring discharge fluctuations and ra infall versus the smoothed spring discharge series. In Figure 47, the cross correlati on coefficient between rainfall and discharge fluctuations is 0.99 for a lag time of ze ro. This strongly sugge sts shortterm spring PAGE 65 67 discharge fluctuations from DCM are linearly correlated to rainfall. A scatter plot of rainfall versus discharge fluctuations is shown in Figure 49. Figure 45. Smoothed spring discharge (calcula ted spring discharge deduct fluctuations) Figure 46. Fluctuations of spring discharge PAGE 66 68 Figure 47. Cross correlation of rainfa ll with fluctuations of spring discharge Figure 48. Cross correlation of rainfall versus smoothed spring discharge PAGE 67 69 Figure 49. Scatter plot of rainfa ll and fluctuations of spring discharge The resultant linear regression equitation is: 556* f luctuationsrainfall (41) where, rainfall has units of ft/d and spring discharge units of ft3/s. From the aquifer recharge area of 2,141,000,000 ft2, the contributing fraction of rainfall to discharge fluctuations is: 02244.0 2141000000 /3600*24*556 (42) Thus, approximately 2.24% of rainfall is simulated in DCM simulated as elevated spring discharge within one day. The remain ing 97.756% is mani fested within the smoothed spring discharge series. Compared to the original cross correlation analysis of rainfall and spring discharge (Figure 44), the current CCF plot of ra infall and smoothed discharge (Figure 48) does not have a peak at lag time zero. The peak is absent because shortterm discharge fluctuations were removed. Similar to the original CCF analysis, there is a cyclical pattern in the plot. The first peak value is occurred when lag time is equal to 74 days. PAGE 68 70 4.3 Geomorphologic Karst Flow Modeling w ith Matrix Flow as the First Flow Segment The previous analysis revealed two com ponents in the DCM simulated spring discharge from the hypothetical karst aquifer. The first represents about 2% of the total spring flow and can be represented as shortterm fluctuations (Figure 46) estimated from a simple linear equation (41). T he second component, the smoothed discharge series (Figure 45), appears contains multiple flow components manifested as multiple hydraulic responses in the discharge series. T he objective of this section is to calibrate a geomorphologic karst model to simula te the smooth disch arge component. To achieve this objective, flow through t he porous matrix of the DCM simulated hypothetical karst aquifer is treated as the first order flow segment in each geomorphologic karst flow model tested. C onduits, however, are r epresented in second, third, and fourth order flow segments (Fi gure 410). Finally the spring corresponds to the fifth order flow segment. T he whole aquifer is in turn treat ed as a fourth order flow system, as shown in Figure 33. The IUH for simulating th is system is equation (326). 4.3.1 Calibration of the Fourth Orde r Geomorphologic Karst Flow Model A wellcalibrated geomor phologic model reproduc es the smoothed spring discharge series from the same rainfall input used in DCM less the volume in the raininduced shortterm discharge fluctuations. Calib ration typically involves the estimation of four parameters including Hortons numbers R Band R L; the mean flow velocity through the basin v; and the characteristic length for the karst aquifer. Here, the length of the highest order conduit, or the on e immediately feeding the spring, is taken to represent the characteristic length of the basin. And from this length, the lengths of all other segments are estimated from R L. In the field, this length of the highest order PAGE 69 71 conduit is often estimated from visible surface features and cave dives. In this analysis, the characteristic length for the hypothetical aquifer is 5.68 miles or the length of the fourth order conduit illustr ated in Figure 410. The DCM simulation produced 5,286 days of spring discharge from a record of an equal number of days of rainfall, the last 903 days of rainfall and simulation discharge are used to calibrate the geomorphic model, while the first 4,383 days of rainfall input are needed to bring the geomorphologic karst aquifer and spring discharge to create baseflow conditions (in effect provide sufficient rainfall history t hat IUH convolution fill the aquifer and spring discharge approaches base flow). Since 2. 24 % of rainfall contributes to the fluctuations of spring discharge, rainfall data used in the geomorphologic models were reduced to 97.76% of values used as input to the DCM model. Figure 410. Order numbers of each c onduit in the hypothetical Karst aquifer 1. Modeling a single flow component. It may be possible to simulate the smoothed DCM spring discharge series a ssuming a single flow component; however, modeling more than one component may be necessary to accommodate the existence PAGE 70 72 of multiple embedded system response times (e.g., the CCF betwe en rainfall and spring discharge suggests multiple response times exis t). Here, the system is simulated as a single flow component, where three paramet ers are calibrated using PEST. The smooth discharge series is shown with the geom orphologic karst model results in Figure 411 for optimum parameter values: 0.009/ vmileday ;1.5 RL ; 5.0 RB (43) Figure 411. Comparison of smoothed spring discharge data from DCM and discharges simulated from the fourth order geomorphologic Karst model assuming one flow component. Figure 411 compares t he simulated spring discharge data from the calibrated geomorphologic karst model wit h the smoothed spring disc harge calculated from the DCM package for a calibration period of 903 days. Clearly, the geomorphologic model cannot capture smoothed variations in sp ring discharge produced from the DCM simulation. PAGE 71 73 Figure 412 does show the simulated spring discharge from the geomorphologic karst model is not sensitive to RB. As resu lt, in all subsequenc e geomorphologic model calibrations, RB is equated 5.0; which simplifi es calibration to estimating two classes of parameters for each flow component: R L and the mean flow velocity through the basin, v. Figure 412. Comparison of results with different RB values 2. Modeling two flow components. Since calibrated results assuming a single flow component geomorphologic karst model were not satisfactory, two groundwater flow components were assumed. Six parameters mu st be specified or calibrated to model two flow components. T he parameters are: t he ratios of rainfall contributing to each flow components: rt1, rt2; Hort ons numbers RL1 and RL2 for each flow components; basinwide mean flow velociti es for both flow components V1 and V2. Because the summation of rt1 and, rt2 equals 0.9776, the number of calibration parameters decreases to five. They are rt1, RL1, RL2, V1 and V2. And rt2=0.9776rt1. PAGE 72 74 The values of parameter calibrated using pest: 10.0051/,10.8702,13.5 20.086/,23.5 VmiledayrtRL VmiledayRL (44) Note that RL is the same for both flow com ponents. Hence, only four parameters need be calibrated if RL can be assumed the sa me for both flow components. Figure 413 shows a comparison of spring discharges gener ated from the calibrated fourth order, two flow components geomorphologic karst flow model and smoothed discharges from DCM for the calibration period. The compar ison is quite good and the average relative error between discharge series is 0.0228. Figure 414 illustrates the instantaneous unit hydrographs (IUH) for the two groundwater flow components of the geomorphologic model. Figure 415 shows the simulated spring discharges from each flow component. The first component has the slowest velocity of .005 miles/d and repres ents 87% of the rainfall volume, whereas the second component is much faster 0.086 mi les/d and dynamic and constitutes 10.7 percent of the rainfall. Combining discharge fluctuations calcul ated by equation (41) to the simulated spring discharges from the geomorphologic ka rst model, will enable direct comparisons to DCM generated discharges (see Figure 416). Figure 416 illustrates that the calibra ted fourth order, two flow components geomorphologic karst model with fluctuations can fit the spring discharge calculated by the DCM package. The average relative error is 0.0242. The final calibrated two flow components geomorphologic karst model with fluctuations is expressed as: 11()1*()*1(1)2*()*2(1)3*()tt iiQtrtpiUtirtpiUtirtpt (45) PAGE 73 75 Figure 413. Comparison of spring discharge from the DCM and fourth order two flow components geomorphologic Karst m odel in the calibration period Figure 414. IUH of the two flow components ca lculated from the four th order two flow components geomorphologic Karst model PAGE 74 76 Figure 415. Spring discharge from the two flow compon ents of the geomorphologic Karst model Figure 416. Comparison of DCM simulated spring discharges to discharges obtained by combining simulations from the fourthordertwocomponents geomorphologic karst model and shortterm fluctuations fr om equation 41. PAGE 75 77 Where Qis spring discharge, p is rainfall, 1Uis the instantaneous unit hydrograph of the first groundwater flow component, 2 Uis the instantaneous unit hydrograph of the second groundwater flow component. 1 rt,2 rt, and 3 rt are the ratios of rainfall contributeed to the first groundwater flow component, the second groundwater flow component, and the fluctuations, respectively. t is the time step. 1,2... i When the dimension rainfall intensity is ft /day, the unit of spring discharge is CFS, and IUH is unit less, the equati on (45) can be written as: 1 1()1*()*1(1)*/24/3600 2*()*2(1)*/24/36003*()/24/3600t i t iQtrtpiUtiarea rtpiUtiareartpt (46) Where the recharge area is 2,141,000,000 square feet fo r the hypothetical karst model (Figure 41), the calibrated values of the ratios of rainfall are rt1=0.8702, rt2=0.97760.8702=0. 1074, rt3=0.0224. 1 Ucan be calculated from equation (326), with the velocity is equal to 0.0051 mile/day and RL1 =3.5; 2U can be calculated from equation (326), with the velocity is equal to 0.0794 mile/day and 23.5 RL To validate the composite geomorphologic karst model (equation 46), DCM and equation 46 were used to simulate dischar ges for an additional 1000 days beyond the calibration period and results compared: to conduct these simulations and additional 1000 days of synthetic rainfall were used. Figure 417 shows the comparison: the first 903 days is the calibration period, after whic h is the validation period. It can be seen from Figure 417 results from the calibra ted model (46) produc e a good approximation of discharges generated from DCM giving an average relative error of 3.1%. PAGE 76 78 Figure 417. Comparison of spring dischar ge from the DCM and results from GEO model for the validation period A cross correlation analyses of rainfall with respect to discharges from each groundwater flow components are show n in Figures 418 and 419. From figure 418, the cross correlation of rainfall and spring discharge from the first groundwater flow component is very low. Howeer, the cross correlation of rainfall and the second flow component is similar to the cross correlation of rainfall with smoothed spring discharge (Figure 48), which is approximately equal to the summation of the spring discharges from the first and second flow components. This suggests that the cross correlation of rainfall with the smoo thed spring discharge is due to the second or quick flow component alone, which cons titutes ~10% of total spring flow. 4.3.2 Analysis of RL and V in the Fourth Order Geomorphologic Karst Model The succes sful calibration of the fourth order two flow components geomorphologic karst model ( equation 46) to the DCM generated discharges, suggests that given rainfall and spring discharge fr om true karst aquife r, one would need the PAGE 77 79 characteristic aquifer length, a single value for Horton number R Land the basinwide average flow velocities V1 and V2 of each flow components to calibrate the geomorphologic model. However, calibra tion does not require that parameters R Land Vbe unique. For example, if parameter values of the first flow component (i.e., 1,1,12 VrtRLrt) are retained, Table 42 shows how 2Vmust vary with respect to2 R L, to maintain a calibrated model. A plot of resu lts from Table 42 is shown in Figure 420 with a fitted regression equation. In a similar analysis, parameter values of the second flow component (i.e., values of 2,1,22 VrtRLrt) can be held constant to examine allowable variations in the velocity of the first flow component velocity 1V with respect to 1 R L (see Table 43). A plot of results from Table 43 is shown in Figure 421, again with a fitt ed regression equation. The exponents of both regression equations are very similar. Values of parameter combinations listed in both Tables 42 and 43 produce the same spring discharge from geomorphologic model (see Figure 422). Figure 418. Cross correlation of rainfall with spring discharge from the first flow component PAGE 78 80 Figure 419. Cross correlation of rainfall with spring discharge from the second flow component Table 42. The calibrated results of 2V for different 2 R L when 1,1,12 VrtRLrt are constant 2 R L 2V 2 R L 2V 2 R L 2V 3.5 0.086067 2.8 0.093319 2.1 0.108 3.4 0.086907 2.7 0.094763 2.05 0.109626 3.3 0.087738 2.6 0.096361 1.9 0.115434 3.2 0.08869 2.5 0.098136 1.8 0.120062 3.1 0.089687 2.4 0.100227 1.7 0.125548 3.0 0.090777 2.3 0.102497 1.6 0.132137 2.9 0.091964 2.2 0.105065 1.5 0.140346 PAGE 79 81 Figure 420. Plot of V2 versus RL2 when the values of V1, rt1, RL1, rt2 are held constant Table 43. The calibrated results of 1V for different 1 R L when 2,1,22 VrtRLrt are constant 1 R L 1V 1 R L 1V 1 R L 1V 3.5 0.005079 2.8 .005465 2.1 .006284 3.4 0.005121 2.7 .005545 2.05 .006376 3.3 0.005166 2.6 .005633 1.9 .006702 3.2 0.005216 2.5 .005732 1.8 .006969 3.1 0.005270 2.4 .005844 1.7 .007288 3.0 0.005329 2.3 .005971 1.6 .007674 2.9 0.005393 2.2 .006115 1.5 .008146 PAGE 80 82 Figure 421. Plot of V1 versus RL1 when the values of V2, rt1, RL2, rt2 are constant 4.4 Treatment of the Porous Medium as the First T wo Flow Segments of a Geomorphologic Model In the previous analysis of geomorphologic karst model the entire porous medium is included in the first order flow segment, which is the only segment receiving recharge. Based on this assumption, the hypothet ical karst aquifer is a fourth order flow system. In this Section, the upper part of the aquifer is assumed to be the first order flow segment, which including the infiltration zone and the epikarst. This first order system is the only flow segment receiving recharge (rainfall). The remaining portion of porous medium is represented as second order flow segments. The conduits network is composed 3rd, 4th, and 5th order flow segments as shown in Figure 423. The spring is a 6th order flow segment. The flow process of the fifth order flow system is shown in Figure 34. PAGE 81 83 Figure 422. Comparison of results from diffe rent values of RL and V of fourth order geomorphologic model 4.4.1 Calibration of the Fi fth Order Geomorphologic Model of the Hy pothetical Aquifer As with the previous analysis, the sm oothed spring discharge series will be used to calibrate the fifth order geomorphologic karst model (equation 342). This discharge represents 97.76% of the rainfall. In addition, the characteristic aquifer length is again the length of the highest order conduit (the fi fth order conduit in Figure 423) or 5.68 miles. 1. Modeling a Single Flow Component: To calibrate the fifth order geomorphologic model for a single flow component, three parameters are calibrated using PEST. Figure 424 is presentation of both the smooth discharge series and simulation results from the fifth order geomorphologic karst model for optimum parameter values: PAGE 82 840.0085/,1.5,5.0VmiledayRLRB (47) Results from the fifth order geomorphologic model are ess entially the same as the fourth order model (Figure 411). Introduci ng an extra order of complexity does not improve results, and parameter estimates are in fact the same. In addition, model simulations continue to be relatively insensitive to the parameter RB. In the following analysis RB are equated to 5.0. Figure 423. Order numbers of each conduit in the hypothetic al aquifer with considering the porous medium as the first two flow segments 2. Modeling Two Flow Components: In this section, the fifth order geomorphologic karst model and two flow com ponents are modeled to simulate spring discharge from the hypothetical karst aqui fer. Using the same 903 day calibration period, optimal parameter estimates are: 10.005053/,10.869808,13.5 20.08558/,23.5V miledayrt RL V miledayRL (48) PAGE 83 85 Figure 424. Comparison of smoothed spring discharge re sulted from DCM and fifth order geomorphologic Karst m odel with one flow component Parameter estimates are essentially the same as th ose determined from calibrating the above fourth order model. In turn, model performance is the same, as are previously observed cross correlation rela tionships between rainfall and discharge, and above described functional relationships between RL1 and V1 and RL2 and V2. Hence, in study a fifth order geomorphologic model does not appear to offer any advantages over a fourth order model with re spect to simulating karst flow. 4.5 Analysis of the Travel Time of Rainfall in the Geomorphologic Models Cross correlation analysis of rainfall ve rsus spring discharge for the hypothetical karst model shows that the lag time between rainfall and smoothed part of spring discharge is 74 days. In addition, cross corre lations of rainfall with the second flow component (quick flow component) (Figure 419) are similar to cross correlations of rainfall to smoothed spring discharges from t he DCM (Figure 48). An analysis of the PAGE 84 86 flow process in geomorphologic model can gi ve estimates of mean travel time for different flow components. 4.5.1 Travel Time Analysis in the Fourth Order Geomorphologic Karst Model For the fourth order geomorphologic model there are four possible paths for rainfall to reach the spring. They are depicted in Figure 433. Knowing the length of the fourth order conduit (Figure 410), which is 5. 68 miles, and the value of RL, is feasible to estimate the mean length of each flow segm ent. From the value of RB, the transition probability ij p can be calculated from equation 323. T he probabilities of rainfall to taking one of four paths can be calc ulated and are shown in Table 45. The average travel distance through the whole system can be estimated by the leng th of every path and its probability. The mean travel time is then estimated using this mean travel length and the mean flow velocity. Table 46 shows this calculation in which: iLis the length of the ithorder flow segment (columns (2 ), (3), (4), (5)); column (6 ) is the length of Path 1, 1LP, which is equal to 1234LLLL ; Column (7) is the length of Path 2, 2LP, which is equal to 124LLL ; Column (8) is the length of Path 3, 3LP, which is equal to 134LLL ; Column (9) is the length of Path 4, 4LP, which is equal to 14LL ; Column (10) is the estimated total travel length (LT) of a rain drop to reach th e spring, which is estimated as 4 1 ii iLTLPp; Column (11) is the flow velociti es of the first flow component (slow flow) 1V calibrated d for different values of 1 R L listed in table 43; Co lumn (12) is the flow velocities of the second flow component (quick flow) 2V calibrated for different values of PAGE 85 872 R Llisted in Table 42; Column (13) is the estimated travel time of the first flow component, which is calculated by 1/1TimeLTV ; and Column (14) is the estimated travel time of the second flow component, which is calculated by 2/2TimeLTV. Travel time for the firs t flow component ranges from 1426 to1439 days, and for the second flow component 83.5 to 84.2days. Figure 425. Four possible paths in fourth order geomorphologic model Table 44. Values of ij p in fourth order geomorphologic model 12 p 13 p 14 p 23 p 24 p 0.71 0.1633 0.1306 0.7333 0.2667 Table 45. Probabilitie s of the four flow paths for the fourth order geomorphologic model Path No. Equations Values Path 1 11223 p pp 0.51778313 Path 2 21224 p pp 0.18831687 Path 3 313 p p 0.1633 Path 4 414 p p 0.1306 PAGE 86 88 4.5.2 Travel Time Analysis for th e Fifth Order Geomorphologic Model The analysis for the fifth order geomorphologi c model is similar to the analysis for the fourth as shown below. Figure 426. Possible paths of th e fifth order geomorphologic model For a fifth order flow process shown in Figure 34, there are eight possible paths that a rain drop can take to reach the spring, shown in Figure 426. The probabilities for the rainfall to take any one of the ei ght paths are shown in Table 48. The average travel distance for a rain drop can be estimated by the length of every path and its probability. The travel time can be estimated using the values of travel length and flow velocity. PAGE 87 89 Table 46. Estimation of travel time of rainfall in fourth order geomorphologic model (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) RL L1 L2 L3 L4 1LP 2LP 3LP 4LP LT V1 V2 Time 1 Time 2 3.5 0.1325 0.4637 1.6229 5.68 7.899 6.2762 7. 4353 5.8125 7.245 0.00508 0.0867 1426.471 84.179 3.4 0.1445 0.4913 1.6706 5.68 7.9865 6.3159 7. 4951 5.8245 7.3091 0.00512 0.0869 1427.287 84.102 3.3 0.1581 0.5216 1.7212 5.68 8.0808 6.3569 7. 5593 5.8381 7.3785 0.00517 0.0877 1428.281 84.066 3.2 0.1733 0.5547 1.775 5.68 8.183 6.408 7. 6283 5.8533 7.4538 0.00522 0.0887 1429.025 84.043 3.1 0.1907 0.5911 1.8323 5.68 8.294 6.4617 7. 7029 5.8707 7.5358 0.00527 0.0897 1429.941 84.023 3.0 0.2104 0.6311 1.8933 5.68 8.4148 6.5215 7. 7837 5.8904 7.6254 0.00533 0.0908 1430.921 84.001 2.9 0.2329 0.6754 1.9586 5.68 8.5469 6.5883 7. 8715 5.9129 7.7236 0.00539 0.0920 1432.158 83.985 2.8 0.2587 0.7245 2.0286 5.68 8.6918 6.6632 7. 9673 5.9387 7.8318 0.00547 0.0933 1433.082 83.924 2.7 0.2886 0.7791 2.1037 5.68 8.8514 .7477 8. 0723 5.9686 7.9514 0.00555 0.0948 1433.974 83.908 2.6 0.3232 0.8402 2.1846 5.68 9.028 6.8434 8. 1878 6.0032 8.0842 0.00563 0.0964 1435.154 83.895 2.5 0.3635 0.9088 2.272 5.68 9.2243 6.9523 8. 3155 6.0435 8.2325 0.00573 0.0981 1436.235 83.888 2.4 0.4109 0.9861 2.3667 5.68 9.4437 7.077 8. 4575 6.0909 8.3989 0.00584 0.1002 1437.187 83.798 2.3 0.4668 1.0737 2.4696 5.68 9.6901 7.2206 8. 6164 6.1468 8.5868 0.00597 0.1025 1438.087 83.776 2.2 0.5334 1.1736 2.5818 5.68 9.9688 7.387 8. 7953 6.2134 8.8004 0.00612 0.1051 1439.143 83.761 2.1 0.6133 1.288 2.7048 5.68 10.2861 7.5813 8. 9981 6.2933 9.0448 0.00628 0.108 1439.334 83.747 2.05 0.6593 1.3516 2.7707 5.68 10.4616 7.6909 9. 11 6.3393 9.1806 0.00638 0.1096 1438.867 83.744 1.9 0.8281 1.5734 2.9895 5.68 11.071 8.0815 9. 4976 6.5081 9.655 0.00670 0.1154 1440.615 83.641 1.8 0.9739 1.7531 3.1556 5.68 11.5626 8.407 9. 8095 6.6539 10.041 0.00697 0.1201 1440.782 83.630 1.7 1.1561 1.9654 3.3412 5.68 12.1427 8.8015 10. 1773 6.8361 10.499 0.00729 0.1255 1440.631 83.628 1.6 1.3867 2.2188 3.55 5.68 12.8355 9.2855 10. 6167 7.0667 11.051 0.00767 0.1321 1440.061 83.633 1.5 1.6829 2.5244 3.7867 5.68 13.6741 9.8874 11. 1496 7.363 11.724 0.00815 0.1403 1439.271 83.539 PAGE 88 90 Table 49 shows this calculation in which: iLis the length of the ithorder flow segment (column (2), (3), (4 ), (5), (6)); column (7) is the length of Path 1, 1LP, which is equal to 12345LLLLL ; Column (8) is the length of Path 2 2LP, which is equal to 1235LLLL ; Column (9) is the length of Path 3, 3LP, which is equal to 1245LLLL ; Column (10) is the length of Path 4, 4LP, which is equal to 125LLL ; Column (11) is the length of Path 5, 5LP, which is equal to 1345LLLL ; Column (12) is the length of Path 6, 6LP, which is equal to 135LLL ; Column (13) is the length of Path 7, 7LP, which is equal to 145LLL ; Column (14) is the length of Path 8, 8LP, which is equal to 15LL Column (15) is the estimated total length (LT) of rainfall to take to reach the spring, which is estimated as 8 1 ii iLTLPp ; Column (16) is the flow velocities of the first flow component (slo w flow) calibrated for different values of 1 R L; Column (17) is the flow velocities of t he second flow component (quick flow) calibrated for different values of 2 R L; Column (18) is the estimated travel time of the first flow component, which is calculated by 1/1 TimeLTV ; and Column (19) is the estimated travel time of the second flow component, which is calculated by 2/2 TimeLTV Travel time for the first flow component ranges from 1426 to1442 days, and for the second flow component 83.3 to 84.1days. Table 47. Values of ij p in the fifth order geomorphologic model p12 p13 p14 p15 p23 p24 p25 p34 p35 0.701205 0.15 0.081305 0.065 0.7061 0.1633 0.13061 0.73333 0.2667 PAGE 89 91 Table 48. Probabilities of the flow paths for the fifth order geomorphologic model Path No. Equations Values Path No. Equations Values Path 1 1122334** p ppp 0.3631 Path 5 11334* p pp 0.1118 Path 2 1122335** p ppp 0.1320 Path 6 11335* p pp 0.0407 Path 3 11224* p pp 0.1145 Path 7 114 p p 0.0813 Path 4 11225* p pp 0.0916 Path 8 115 p p 0.0650 4.5.3 Application of La place Transforms to Estimate the Travel Time The travel times of a given flow com ponent can be calculated directly from the analytical equation for the normalized central moment of the IUH. From the Laplace transform of the IUH function (326) for fourth order geomorphol ogic model, temporal moment equations can be deriv ed. The Laplace transform ( sc) to the IUH function (326) is: 316 16 2 *16 3 316 2 216 1 1161 1 1 1 1 1 s f s e s e s d s c s bIUHLsc (49) The moment equations tnm, can be determined from the Laplace domain solution by utilizing the following relationship. n n s n tnds scd mlim0 ,1 (410) Then, scms t0 ,1lim = 1 *16 2 *16 1 316 1 216 1 1161 fedcb (411) ds scd ms tlim0 1 ,01=16161616fdcb (412) PAGE 90 92 Table 49. Estimation of travel time of rainfall in fifth order geomorphologic model (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17) (18) (19) RL L1 L2 L3 L4 L5 1LP 2LP 3LP 4LP 5LP 6LP 7LP 8LP LT V1 V2 Time 1 Time 2 3.5 0.0379 0.1325 0.4637 1.6229 5.68 7. 9369 6.314 7.4732 5.8503 7.8044 6.1815 7.3407 5.7179 7.1994 0.00505 0.08558 1425.63 84.12 5 3.4 0.0425 0.1445 0.4913 1.6706 5.68 8.029 6..358 4 7.5376 5.867 7.8844 6.2139 7.3931 5. 7225 7.2625 0.00509 0.08655 1426.81 83.91 1 3.3 0.0479 0.1581 0.5216 1.7212 5.68 8. 1287 6.4075 7.6072 5.8859 7.9707 6.2495 7.4491 5.7279 7.3309 0.00514 0.08739 1426.24 83.88 7 3.2 0.0542 0.1733 0.5547 1.775 5.68 8. 2372 6.4622 7.6825 5.9075 8.0639 6.2889 7.5092 5.7342 7.4054 0.00519 0.0883 1426.85 83.86 6 3.1 0.0615 0.1907 0.5911 1.8323 5.68 8. 3555 6.5232 7.7644 5.9322 8.1648 6.3326 7.5738 5.7415 7.4868 0.00524 0.08934 1428.78 83.80 1 3.0 0.0701 0.2104 0.6311 1.8933 5.68 8. 4849 6.5916 7.8538 5.9605 8.2746 6.3812 7.6435 5.7501 7.5761 0.0053 0.09044 1429.46 83.77 0 2.9 0.0803 0.2329 0.6754 1.9586 5.68 8. 6272 6.6686 7.9518 5.9932 8.3943 6.4357 7.7189 5.7603 7.6746 0.00537 0.09164 1429.16 83.74 7 2.8 0.0924 0.2587 0.7245 2.0286 5.68 8. 7842 6.7556 8.0597 6.0312 8.5255 6.4969 7.801 5.7724 7.7835 0.00544 0.09297 1430.80 83.72 1 2.7 0.1069 0.2886 0.7791 2.1037 5.68 8. 9583 6.8546 8.1792 6.0755 8.6697 6.566 7.8906 5.7869 7.9047 0.00552 0.09445 1432.01 83.69 2 2.6 0.1243 0.3232 0.8402 2.1846 5.68 9. 1523 6.9677 8.3121 6.1275 8.8291 6.6445 7.9889 5.8043 8.0402 0.00561 0.0961 1433.19 83.66 5 2.5 0.1454 0.3635 0.9088 2.272 5.68 9. 3697 7.0977 8.4609 6.1889 9.0062 6.7342 8.0974 5.8254 8.1926 0.00571 0.09795 1434.78 83.64 1 2.4 0.1712 0.4109 0.9861 2.3667 5.68 9. 6149 7.2482 8.6287 6.2621 9.204 6.8373 8.2179 5.8512 8.3652 0.00582 0.10003 1437.32 83.62 6 2.3 0.203 0.4668 1.0737 2.4696 5.68 9. 8931 7.4235 8.8194 6.3498 9.4263 6.9567 8.3525 5.883 8.5619 0.00596 0.1025 1436.57 83.53 1 2.2 0.2425 0.5334 1.1736 2.5818 5.68 10. 211 7.6295 9.0377 6.4559 9.6778 7.096 8.5043 5.9225 8.7881 0.00611 0.10524 1438.31 83.50 5 2.1 0.2921 0.6133 1.288 2.7048 5.68 10. 578 7.8734 9.2901 6.5854 9.9648 7.26 8.6768 5.9721 9.0502 0.00629 0.10841 1438.83 83.48 1 2.05 0.3216 0.6593 1.3516 2.7707 5.68 10. 783 8.0125 9.4316 6.6609 10.124 7.3532 8.7723 6.0016 9.1974 0.00638 0.11019 1441.61 83.46 9 1.9 0.4358 0.8281 1.5734 2.9895 5.68 11. 507 8.5174 9.9334 6.944 10.679 7.6893 9.1053 6.1158 9.7204 0.00674 0.11649 1442.20 83.44 4 1.8 0.5411 0.9739 1.7531 3.1556 5.68 12. 104 8.9481 10.351 7.195 11.13 7.9742 9.3766 6.2211 10.156 0.00704 0.12171 1442.56 83.44 1 1.7 0.6801 1.1561 1.9654 3.3412 5.68 12. 823 9.4816 10.857 7.5162 11.667 8.3255 9.7012 6.3601 10.684 0.0074 0.12822 1443.83 83.32 8 1.6 0.8667 1.3867 2.2188 3.55 5.68 13. 702 10.152 11.483 7.9334 12.315 8.7654 10.097 6.5467 11.337 0.00786 0.13607 1442.34 83.31 6 1.5 1.122 1.6829 2.5244 3.7867 5.68 14. 796 11.009 12.272 8.4849 13.113 9.3264 10.589 6.802 12.156 0.00843 0.14589 1442.05 83.32 6 PAGE 91 93 The central moment tu,1can be calculated by t t tm m u,0 ,1 ,1 16161616 1 116 1 116 1 116 1 116 1 116fdcb fedcb (413) The central normalize moment tu,1represents the mean re sponse time of the system. Substituting the expressions of 16,16161616,,, fedcb(equation 318), and the values of parameters R Land vin Table 46 into equation (413), values of the central moment tu,1of the IUH of different flow com ponents for different combinations of R Land vcan be calculated and shown in Table 410. It can be seen that values of tu,1 for different combinations of R Land vare nearly same. The central moment of IUH for the first flow component is around 1430 days, and the central moment of IUH for the second flow component is around 84 days. These values are consistent with the response times calculated from the geomor phologic model (see Table 46). Similar analysis can be performed for the fifth or der geomorphologic model. PAGE 92 94 Table 410 Central moments of IUH for t he two flow components calculated from the fourth order geomorphologic model RL V1 V2 Central moment of IUH2 (days) Central moment of IUH1 (days) 3.5 0.00508 0.0867 83.593 1426.677 3.4 0.00512 0.0869 84.138 1428.052 3.3 0.00517 0.0877 84.162 1427.662 3.2 0.00522 0.0887 84.062 1428.417 3.1 0.00527 0.0897 84.040 1430.428 3 0.00533 0.0908 84.009 1431.141 2.9 0.00539 0.092 83.981 1433.444 2.8 0.00547 0.0933 83.971 1432.261 2.7 0.00555 0.0948 83.904 1433.172 2.6 0.00563 0.0964 83.890 1436.409 2.5 0.00573 0.0981 83.948 1437.227 2.4 0.00584 0.1002 83.850 1438.663 2.3 0.00597 0.1025 83.803 1438.82 2.2 0.00612 0.1051 83.762 1438.459 2.1 0.00628 0.108 83.777 1440.744 2.05 0.00638 0.1096 83.793 1439.457 1.9 0.0067 0.1154 83.694 1441.539 1.8 0.00697 0.1201 83.632 1441.070 1.7 0.00729 0.1255 83.689 1440.731 1.6 0.00767 0.1321 83.685 1441.308 1.5 0.00815 0.01403 83.595 1439.060 PAGE 93 95 CHAPTER 5 GEOMORPHOLOGIC KARST MODEL WI TH COMMON CONDUIT FLOW FOR DIFFERENT FLOW C OMPONENTS In chapter 4, it is s een that the best geomorpholog ic modeling results are obtained when separate slow and fast flow components are considered in parallel. Furthermore, simply increasing the order of t he model (i.e. from fo urth to fifth order) did not improve calibration; in fact, m odel permanence, validat ion performance, and parameter estimates are ess entially the same between fourth and fifth order models. In this chapter, a novel geomorphologic m odel is developed in which fast and slow flow components share the same conduit network. The conceptual design of this new model is believed to share greater similarity to natural karst systems. 5.1 The Conceptual Model of Comm on Conduit Flow for Different Flow Components Groundwater can travel rapidly in karst c onduits but relatively slow in porous medium. Initially rainfall enters a karst aquifer as slow matrix flow, and then the travel velocity increases rapidly once the flow reaches a conduit and enters the conduit network. If recharge takes a long tr avel path as laminar flow through porous medium, the travel time to a conduit will ty pically exceed travel time in the conduit alone and total travel time can be quite long. If recharge takes a relatively short path through the porous matrix to reach a conduit or recharge occurs directly into a conduit (as with sinking stream s from autogenic or allogenic sources), travel time in the conduit may dominate and in turn produce flow with relatively brief hydraulic residence times in the aquifer. Hence, all fl ows (fast and slow) spend travel time as conduit flow and the variance on conduit travel time is probably much smaller than the variance on travel time through porous medium. PAGE 94 96 As shown in Figure 51, the slow groundwater flow component can be considered a combination of a longer trav el times or longer tr avel paths as flow through porous media and then conduit flow. While, the fast groundwater flow component can be considered as a combinati on of a brief travel times or short travel distances as flow through porous m edia and then conduit flow. For the slow flow component, longer travel times or l onger travel paths as flow through porous media are defined by the flow velocity in porous medium Vm and the mean length of a long travel path through the porous mediumLml For the fast flow component, the brief travel times or short travel distances as fl ow through porous media are defined by the flow velocity in porous medium Vmand the mean length of a short travel path through the porous mediumLms. For both flow components, the conduit flow is defined as the flow velocity in the conduit networkVc, and the geomorphologic parameters of karst conduits ( L R L, R B). As to the fast flow component, consider fi rst an extreme condition, in which, the travel path through the porous medium tak en by the rainfall to reach the conduit network is very short and can be neglected, as shown in the dash line box in Figure 51. Then a model of the fast flow co mponent needs to include travel time in conduits alone. From the conceptual model of flows thr ough a karst aquifer illustrated in Figure 51, a new geomorphologic m odel is needed to simulate these flows in a karst aquifer such as the hypothetical aquifer shown in Figure 41. The slow flow component is essentially a fourth order fl ow system illustrated in Figure 52(a); PAGE 95 97 whereas, the fast flow com ponent is a third order system conceptualized in Figure 52(b). Figure 51. Two different flow component s that share a common conduit network To create a geomorphologic karst flow model that will simulate the distribution of flows depicted in Figure 52, several adaptati ons must be made to the system of model equations developed in c hapter 3. These adaptations include: 1. Modification of the fourth order geomorphologic model to simulate the slow flow: a. Two flow velocities will be used to describe the slow flow component. One for the mean velocity of c onduit flow and another for the mean velocity of flow through porous media. b. The length of the trav el path in porous medium Lml is independent of the length of length of the highes t order conduit and the Hortons numbers RLand RB of the network; and, c. The probabilities of a rain drop flo wing from the first flow segment (the porous medium) to others flow segments (conduits of differing order) are determined by the mean length of each order of karst conduit The original equation for the IUH of the fourth order geomorphologic karst model is equation (326). The above adaptat ions are manifested in equation 326 as several specific changes. First, the hydraulic residence time in the first flow segment (the porous medium) 1 is now defined as: PAGE 96 98 11/ VmL (51) (a) Flow process of t he slow flow component (b) Flow process of t he fast flow component Figure 52. Conceptual models for the slow (a) and fast (b) components of a geomorphologic karst flow model usi ng a common karst conduit network PAGE 97 99 where, Vmis groundwater flow in porous medium; 1Lis the mean travel length in the porous medium. The characteristic length of the aquifer is continues to be the length of the highest order conduit (in this case the 4th order) 4L, and the length ratio of the conduit networks is defined by R L. From this the m ean length of the 3rd flow segment 3L (the 2nd order karst conduits) is: 34/ LLRL In addition, the mean length of the 2nd flow segment 2L (the 1st order karst conduits) is: 2 24/ L LRL Thus, the probabilities of a rain drop flo wing from the first flow segment (the porous medium) to the others flow segm ents (conduits of different order) are: 2 2 12 2 234* ** LRB p LRBLRBL (52a) 3 13 2 234* ** LRB p L RBLRBL (52b) 4 14 2 234** L p L RBLRBL (52c) Substitute equations (51) and (52) into equation (326) and retaining all other terms produces the IUH for the slow flow component shown in Figure 52 (a). 2. Modification of the th ird order geomorphologic model to simulate the fast flow: The fast flow component (Figuer 52(b) ) can be modeled as a third order flow system. The flow path in the porous medium is assumed very short and neglected. Rainfall is assumed to enter the karst c onduit network directly. The resulting third PAGE 98 100 order model differs from t he one developed in chapter 3. First, all of the flow segments receive rainfall and not just the first. Second, the distribution of rainfall received among conduits of different order, depends on the total length of conduit in each order. Given the length of the highes t order (3rd) conduit (the 4th flow segment in Figure 52) is 4L, the conduit length ratio R L, and the conduit number ratio is R B, The total length of conduit in each order is as follows: The total length of the 3rd conduit (4th order flow segment) is 4L; The total length of the 2nd conduit (3rd order flow segment) is 4 L RB R L; The total length of the 1st conduit (2nd order flow segment) is 2 4 2LRB R L. Hence, the initial conditions of th e different flow segments are: 2 2 4 22 2 22 44 4 2 20 1 LRBRB RLRL LRBLRBRBRB L R LRL R LRL (53a) 4 3 22 44 4 2 20 1 LRBRB RLRL L RBLRBRBRB L R LRL R LRL (53b) 4 4 22 44 4 2 21 0 1 L LRBLRBRBRB L R LRL R LRL (53c) From the equation definin g third order geomorphologic model (equation 33), the IUH of the fast flow component is: PAGE 99 101 52 53 54 5 234000fdtdtdtdt IUHt dtdtdtdt (54) From the equations (221), (221) and (223), it is known that: ** 3 24425 1 2 3 41t ttttAeAeAteAe (55) 3 442 4343 43 35 22 43 43432 ()1t tttee te (56) 4445 4()1tttete (57) where, 2 42243 1 2 3242p A 2 4223 2 2 3243p A ** 4232424 3 ** 2443p A 2 43 2 2 4 2 4 24 3 42 2 43232 42 43 2 4 432 4242 3 44223 p p A 2 23 222 2RBRB p R BRB PAGE 100 102 2 24 232 2RBRB p R BRB in which, 22/VcL ,1 32 R L ,2 42 R L iL represents the mean length of the flow segment of orderi,Vcis the flow velocity in conduit network. With the IUH defined for both flow com ponents, the smoothed spring discharge of the hypothetical karst aquifer (the fluct uations extracted from the original DCM simulation) can be calculated as: 11()1*()*1(1)2*()*2(1)tt iiQtrtpiUtirtpiUti (58) Where Qis spring discharge, pis rainfall, 1 Uis the instantaneous unit hydrograph of the slow groundwater flow component, 2 Uis the instantaneous unit hydrograph of the fast flow component. 1 rt and 2 rtare the respective ratios of rainfall contributed to the slow flow component and the fast flow component. tis time step. 1,2...i 5.2 Calibration of Common C onduit Geomorphologic Model Given, the previously hypothetical ka rst aquifer depicted in Figure 41, the DCM simulated spring discharges the smoothed discharge series (fluctuations removed) and the ancillary dat a on the recharge area (2141000000 ft2), the length of the highest order conduit (5.68 miles), and the sum contribution of rainfall to the fast and slow flow components (97.756%), the new geomorph ologic model (58) was calibrated using PEST and the smoothed spring discharge series. The optimum calibrated values of the parameters were: PAGE 101 103 Mean length of the porous medium: mile Lm 0023.7 ; Velocity of porous medium flow: daymile Vm / 006002.0 ; Ratio of rainfall to the slow component: mile rt 9341.01 which gave The ratio of rainfall to the conduit flow: mile rt 04346.01 The 50% less flow is associated with the fast flow components in the new geomorphologic model the found with the fourth and fifth order models developed in Chapter 4. In addition, as with the other models, calibrated values for R L and Vc are not unique. Table 51 lists calibrated values of cV and R L. Figure 53 shows the plot of velocity Vcversus the Hortons number R L. Knowing the mean length taken by slow flow component to reach the conduit network and the mean velocity of the slow fl ow component, the mean travel time to reach the conduit network equates to: days VmLm 66.1166/ (59) During the calibration, it is found that t he calculated values of Lm and Vm are not unique. Mathematically, every combinat ion of Lm and Vm satisfying equation (59) produce same results. Simulations results from the 903 day calibration period are shown and compared with smoothed spring discharge calculated from DCM in Figure 54. The relative error is around 0.0103, which is half that of the fourth order and fifth order models in chapter 4 (the previous error is around 0.025). Figure 55 illustrates instantaneous unit hydrographs for both flow components. Figure 56 shows the separate discharge seri es for both flow components for entire calibration. The both slow and fast flow components of the new geomorphologic PAGE 102 104 model demonstrate greater flow variability than respective components in fourth and fifth order models devel oped in Chapter 4. Table 51 Calibrated results of Vc for different R L for fifth order geomorphologic karst model R L Vc R L Vc R L Vc 3.5 0.269 2.8 0.291 2.1 0.333 3.4 0.271 2.7 0.295 2.05 0.337 3.3 0.274 2.6 0.3 1.9 0.352 3.2 0.277 2.5 0.305 1.8 0.363 3.1 0.28 2.4 0.311 1.7 0.377 3.0 0.283 2.3 0.317 1.6 0.392 2.9 0.287 2.2 0.324 1.5 0.41 Figure 53. Plot of Vc versus RL PAGE 103 105 To validate the combined model, i.e. one incorporating bot h geomorphologic model (58) and the shortterm discharge fl uctuations model (41), additional simulations were conducted and results compared to an additional 1000 days of DCM simulated spring discharges. Figure 57 shows the comparison. In Figure 57, the first 903 da ys in Figure 57 is the calibration period, after this, it is the validation period. It can be seen that the results from the calibrated karst model fit simulation results generated from the DCM in both calibration and validation periods. For the calibration period, the average relative error is 1.03%; for the validation period, the aver age relative error is 2.13%. Again this error is smaller than results generated from fourth and fifth order geomorphologic models calibrated to the same hypothetical spring discharge data in chapter 4. From these results, it can be concluded that the new geomorphologic karst model, in which, different flow components have same conduit network, per form better than the geomorphologic karst models applied in chapter 4. Figure 54. Comparison of results from the new geomorphologic model and DCM PAGE 104 106 5.3 Application of the Common Condui t Geomorphologic Model to a Re vised Conduit Network 5.3.1 The Revised Conduit Network Karst flows through the revised conduit net work is shown in Figure 58 were modeled using MODFLOWDCM. The rechar ge area of the revised hypothetical aquifer system was 2,083,000,000 ft2. Daily rainfall data used as recharge are shown in Figure 59. There are 5,286 days of data in this simulation. The resulting simulated spring discharge is shown in Figure 510. (a) the IUH of the sl ow flow component (b) the IUH of the fast flow component Figure 55. Instantaneous unit hydrographs of the two flow component PAGE 105 107 (a) the spring discharge from the slow flow component (b) the spring discharge from the fast flow component Figure 56. Simulated spring dischar ge from the two flow components PAGE 106 108 Figure 57. Validation of the Geomorphologic model us ing common conduits for different flow components Figure 58. Revised conduit network 5.3.2 Cross Correlation Analysis of Rainfall with Calculated Spring Discharge Figure 511 displayss the cross correlation function of ra infall (Figure 59) versus spring discharge (Figure 510). In Figure 511, the function peaks at lag time is equal to zero. PAGE 107 109 Figure 59. Daily rainfall data used in the hypothetical karst model Figure 510. Spring dischar ge calculated by the DCM PAGE 108 110 Figure 511. Cross correlation between rainfall and spring discharge of the hypothetical Karst groundwater system As with previous analyses, a high cro ss correlation coefficient at lag times equal to zero suggests the presence of a flow component which is a linearly correlated with rainfall (see Figure 511). Th e discharge fluctuations were removed from the original spring discharge as de scribed previously gener ating the smoothed series and fluctuations (see Figure 512 and Figure 513). Figure 514 shows the plot of cross correla tions of rainfall versus fluctuations of spring discharge, and in Figure 515 cross correlations of rainfall and the smoothed spring discharge series. In Figure 514, the cross correlation coefficient of rainfall with the fluctuations of rainfall is 0.99 for a lag time equal to zero. A scatter plot of rainfall and fluctuations data is shown in Figure 516. The linear equation that approximates values of the discharge fluctuation series is: fallrain nsfluctuatio *5.424 (510) PAGE 109 111 Figure 512. Smoothed spring dischar ge in the calibration period Figure 513. Fluctuations of spring discharge in the calibration period PAGE 110 112 Figure 514. Cross correlation of rainfall with Fluctuations of spring discharge Figure 515. Cross correlation of rainfall with smoothed spring discharge PAGE 111 113 Figure 516. Scatter plot of rainfa ll and fluctuations of spring discharge The unit for rainfall is ft/day and the unit for spring discharge is ft3/second. Knowing the recharge area of the revised hypothetical karst aquifer is 2,083,000,000 square feet, the ratio of rainfall that contri butes to the fluctuations of spring discharge is: 01761.0 2083000000 /3600*24*5.424 (511) This means that ~1.76% of rainfall flows through the karst aquifer and becoming spring discharge within one day. T he other 98.239% of rainfall produces the smoothed spring discharge component. Figure 515 shows the CCF plot of rainfall with smoothed spring discharge.. Similar to the CCF plot of rainfall with or iginal spring discharge, the CCF plot of rainfall with smoothed spring di scharge seems to be cyclical. PAGE 112 114 5.3.3 Calibrate Common Conduit Geomor phologic Karst Model for the Revised Conduit Netw ork Although the revised conduit network is simpler than the first hypothetical conduit network (Figure 4.1), it possesse s a third order conduit flow system. The order number of each conduit is shown in Fi gure 58. The length of the third order conduit, which is the highest order conduit, is 10.908 miles. In applying the common conduit geomorphologi c model to simulate the spring discharge, the fast flow component (conduit flow) is a third order flow system shown as figure 2(b). The slow flow component is a fourth order flow system shown as figure 2(a). The calibrated model parameters are: The mean l ength of the porous medium: mile Lm 0.7 The velocity of the porous medium flow: daymile Vm /0046.0 The ratio of ra infall to the slow component: mile rt 9401.01 Then, the ratio of rainfall to the conduit flow: mile rt 0423.02 The mean travel time of the slow flow component to reach the conduit network equates to: Lm / Vm=1521.74 days (512) Mathematically, Lm and Vm are not unique, but as long as they satisfy the equation (59), they produce same spring discharge results. Table 52 shows the calibrated results of cVfor different values of R L. Figure 517 shows the plot of velocity Vcversus the Hortons number R L. The calculated results from the calibrated geomorpholog ic karst model and the smoothed spring PAGE 113 115 discharge calculated from DCM are compared in Figure 518. The relative error is around 0.0119, which is small and acceptable. Table 52 Calibrated results of Vc for different R L for fifth order geomorphologic karst model R L Vc R L Vc R L Vc 3.5 0.3888 2.8 0.4162 2.1 0.4704 3.4 0.3925 2.7 0.4208 2.05 0.4753 3.3 0.3946 2.6 0.4273 1.9 0.4962 3.2 0.3991 2.5 0.433 1.8 0.5118 3.1 0.403 2.4 0.4408 1.7 0.5291 3.0 0.4061 2.3 0.4492 1.6 0.5499 2.9 0.4115 2.2 0.4598 1.5 0.5752 Figure 517. Plot of Vc versus RL Figure 519 shows the instantaneous unit hydrographs of the two flow component. Figure 520 shows the spring di scharge from the two flow components in the whole calibration period. PAGE 114 116 Figure 518. Comparison of result s from geomorphologic model with DCM In order to validate the geomorphologic model, results of the calibrated model are combined with the fluctuations calc ulated by equation (5 10) for direct comparison to DCM results. An additio nal 1000 days of simulation results were generated using DCM and the geomorphologic model following the calibration period. The simulations generated from both models appear in Figure 521. The results from the calibrated karst model co mpare well with results from the DCM for validation period. For the calibration period, the relative error is 1.119%; for the validation period, the relative error is 2.26%. PAGE 115 117 (a) the IUH of the slow flow component (b) the IUH of the fast flow component Figure 519. Instantaneous unit hydrogr aphs of the two flow component PAGE 116 118 Figure 520. Spring discharge fr om the two flow components Figure 521. Validation of the Geomorphologic model using common conduits for different flow components PAGE 117 119 CHAPTER 6 ANALYSES OF OBSERVED DATA IN THE BLUE SPRING AREA AND GEOMOR PHOLOGIC MODEL CALIBRATION In this chapter, time se ries analysis is applied to analyze the observed data of rainfall and groundwater levels near Blue Springs and the spring discharge itself. Following this analysis, a geomorphologic kars t model is calibrated to simulate Blue Spring discharge. The series analysis is critical as it lays the groundwork for calibrating a geomorphologic model and fo r interpreting the model results. Blue Spring lies in Volusia County in Florida. The geology of the study area consists of carbonate formations overlain by surfical classic s ediments. There are three recognized significant aquifers in the study region: the surficial aquifer, the upper Floridan aquifer, and the lower Floridan aquifer. there are two confining layers, the intermediate confining layer overlyi ng the Floridan aquifer and the middle confining unit separating the upper floridan aquifer from the lower floridan aquifer. Blue Spring is a 1st magnitude spring (discharge greater than or equal to 100 Cfs). The longterm average discharge of Blue Spring is 157 cfs calculated from 654 instantaneous manual flow m easurements over a 75year period of record (St, Johns River Water Management District, 2008; Green, 2008). 6.1 Observation Data Analysis for Blue Spring Area The available time series data near Blue Spring area include rainfall data from rainfall station DELAND, discharge rates at Blue Spring, and groundwater levels at four observation wells: V1091, V0867, V0083, and L0059. Figure 61 shows the map of the study area, incl uding wells and Blue Spring locations. Figure 62 shows observed rainfall data at DELAND and Bl ue Spring discharge. Figure 63 presents observed water levels at t he four observation wells. PAGE 118 120 Figure 61. Study Area PAGE 119 121 Figure 62. Available observ ed data of rainfall and spring discharge in Blue Spring area Figure 63. Available data of groundwater levels at four observation wells and rainfall 6.2 Cross Correlation Analysis of Data from the Blue Spring Area 6.2.1 Autocorrelation Analysis The autocorrelation functions for rainfalls, spring discharge, and groundwater levels are shown in Figure 64. The autoco rrelation functions of the rainfall series decrease sharply for lag values greater than one day. This means that little or no PAGE 120 122 memory effect is detected in rainfall after a few days, after which the rainfall signal appears fairly random. Longer memory effect s, however, are observed for spring discharge and groundwater levels. Figure 64. Autocorrelation functions of rainfall, water level and spring discharge 6.2.2 Correlation Analysis of Rain fall w ith Blue Spring Discharge Figure 62 shows available rainfall dat a at Deland and Blue Spring discharge. Based on available continuous records, t he observed spring discharge series was divided into two periods. The first extends from 12/8/2001 to 9/7/2004 and the other from 3/21/2005 to 8/21/2006. Cross correla tions between rainfall at Deland and Blue Spring discharge were subsequently analyzed based on these two periods. Period 1: (12/8/20019/7/2004): There are 1,005 records in this period. The maximum lag tested is 250 days. Before cross correlation analysis was performed, a PAGE 121 123 sevendays centered moving average trans formation was performed to smooth observed data (Figure 65). The peak value of the cross correlation f unction (CCF) of rainfall versus spring discharge calculated from the sevendays moving averaged series was 0.198, for a lag equal to 153 days. From the equation (211), lag time can be calculated based on the plots of phase and coherence functions. Figure 66 depicts the spectral density of spring flow. The spectrum density function exhi bits a peak value when the frequency is equal to 0.00312 cycles per day. The peak val ue indicates a significant amount of the variance, having a 320 days periodicity, wh ich indicates that the spring discharge may have an annual periodic component. The coherence spectrum and the phase spec trum are plotted in Figure 67. The coherence function expresses the in tensity of the input output link and its linearity. The average value of coherenc e between rainfall and spring discharge is around 0.3. Such a small coherence va lue suggests a nonlinear relationship between spring discharge and rainfall. From Figure 67, the phase values co rresponding to the two frequencies, at which the values of coherence are the maximum and second maximum, are 2.8 and 1.1, corresponding to 1.1 56 0.00312*2*3.14 k days, 2.8 143 0.00312*2*3.14 k days, respectively. These results are consistent with the aforementioned cross correlation analysis between rainfall and spring discharge. PAGE 122 124 (a) Daily data of the analyz ed period (total: 1,005 days) (b) CCF of the analyzed period Figure 65. Cross correlation analysis of rainfall with spring discharge from 12/8/2001 to 9/7/2004 PAGE 123 125 0.50 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 Frequency 4.424E5 1.628E5 5.987E4 2.203E4 8.103E3 2.981E3 1.097E3 4.034E2 1.484E2 5.46E1 2.009E1 7.389E0 1.0E0 Density Figure 66. Spectral density of spring Period 2: (3/21/2005 to 8/21/2006) There are 467 records in this period (F igure 68), and the maximum lag tested is 250 days. There are two peak values in the plot of cross correlation coeffi cients of rainfall versus spring discharge in this period (Figur e 68(b)). The first peak value of cross correlation for rainfall with spring discharge is 0.348, for a lag equal to 58 days; the second value is 0.218, for a lag equal to 159 days. The lag time (153 days) calculated from the first period is very similar to the lag time (159 days) from the se cond period. However, besides this lag time, there is another lag time (58 days) in the second period analysis. Furthermore, the cross correlation functions of rainfall and spring discharge have higher peak values for the second period than for the first period. PAGE 124 126 0.52 0.50 0.48 0.46 0.440.42 0.40 0.38 0.36 0.340.320.300.28 0.26 0.24 0.22 0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.060.040.020.00 0. 02 Frequency 1.04 0.98 0.92 0.86 0.80 0.74 0.68 0.62 0.56 0.50 0.44 0.38 0.32 0.26 0.20 0.14 0.08 0.02 0.04 Coherency (a) Coherence of rainfall and spring (b) Phase spectrum of rainfall and spring Figure 67. Coheren cy and Phase diagrams PAGE 125 127 (a) Daily data of the analyzed perio d (total: 467 days) (b) CCF of the analyzed period (c) Figure 68. Cross correlation analysis of rainfall with spring discharge from 3/21/2005 to 8/21/2006 Figure 69 shows the monthly average ra infall data for the year from 2001 to 2006. Normally, for each year, the rainfall ra te is lowest during January and the most rainfall occurs during July, Aug, and Septem ber. Cross correlations between rainfall PAGE 126 128 and Blue Spring discharges for each year are analyzed separately and are shown in Figure 610, 611, 612, 613 and 614. Figure 69. Monthly averaged rainfa ll data from January 2001 to May 2006 (a) (b) Figure 610. Cross correlati on analysis of rainfall with spring discharge in 2002, (a) Daily data of 2002 (total: 365), (b) CCF of rainfall with spring discharge in 2002. PAGE 127 129 (a) (b) Figure 611. Cross correlati on analysis of rainfall with sp ring discharge in 2003. (a) Daily data of 2003 (total: 365), (b) CCF of rainfall with spring discharge in 2003. (a) (b) Figure 612. Cross correlati on analysis of rainfall with sp ring discharge in 2004. (a) Daily data of 2004 (total: 251), (b) CCF of rainfall with spring discharge in 2004. PAGE 128 130 (a) (b) Figure 613. Cross correlati on analysis of rainfall with sp ring discharge in 2005. (a) Daily data of 2005 (total: 286), (b) CCF of rainfall with spring discharge in 2005. (a) (b) Figure 614. Cross correlati on analysis of rainfall with spring discharge in 2006. (a) Daily data of 2006 (total: 181), (b) CCF of rainfall with spring discharge in 2006 PAGE 129 131 For the year 2002 (Figure 610), the peak value of CCF for rainfall with spring discharge is 0.314, for a lag equal to 104 da ys. For the year 2003 (Figure 611), all of CCF values are very low; however, it app ears that there are tw o peaks, the first is 0.015, for a lag equal to 35 days; the second is 0.127, for a lag time equal to 154 days. For the year 2004 (Figure 612), ther e is no obvious peak, and CCF values are around 0.2, which equates to lags bet ween 130 and 160 days. For the year 2005 (Figure 613), the first CCF peak is 0.43, for a lag equal to 56 days; the second 0.103, for a lag equal to 158 days. For t he year 2006 (Figure 614), the peak is around 0.38, for a lag time of 5060 days. (a) (b) Figure 615. Cross correlation analysis of rainfall with spring discharge from 6/23/2002 to 9/26/2002. (a) Daily data of the analyzed per iod (total: 96 days), (b) CCF of the analyzed period. The results of cross correlation analysi s of rainfall with spring discharge for different years are quite di fferent. The karst aquifer mi ght have different hydraulic responses for different recharge conditions. To evaluate this concern, cross correlations were evaluated between rain fall and spring discharge during wet PAGE 130 132 seasons for the years 2002, 2004 and 2005. Because the rainfall amounts in 2003 were very low, and data for 2006 is only available through 6/30, data from 2003 and 2006 are not included in t he wet season analysis). (a) (b) Figure 616. Cross correlation analysis of rainfall with spring discharge from 7/26/2004 to 9/7/2004. (a) Daily data of the analyzed pe riod (total: 44 days), (b) CCF of the analyzed period. (a) (b) Figure 617. Cross correlation analysis of rainfall with spring discharge from 7/6/2005 to 9/10/2005. (a) Daily data of the analyzed pe riod (total: 67 days), (b) CCF of the analyzed period PAGE 131 133 Results of cross correlation analysis of rainfall with spring discharge for wet seasons in 2002, 2004 and 2005 are very diff erent. For the analysis of wet seasons in 2002 and 2005, there are no obvious peak values in the CCF plots of rainfall versus Blue Spring discharge, when the test lag time is only 50 or 60 days. This might be because it took more than 60 days for most rainfall events, during the wet seasons for these two years, to influenc e the Blue Spring discharge. However, for the analysis of the wet season in 2004, there are two peak values in the CCF plot. The first is 0.465, for a lag equal to 1 day ; the second CCF peak of value 0.168, for a lag equal to 10 days. The high cross co rrelation between rainfall and spring discharges for the wet season in the year 2004 might be because of extremely high rainfall volumes (flood condition). 6.2.3 Correlation Analysis of Rainfall Versus Groundwater Levels Figure 63 depicts available data of groundwater levels at the four observation wells and rainfall. It can be seen that the data of groundwater le vels at the four observation wells hav e similar trends. First, the cross correlation between groundwater levels at different wells is analyzed and shown in Figure 618. The peak values of cross correlation of groundwater levels at the wells V1091,V0867 and L0059 with groundwater levels at the well V0083 are very high (near 0.99), and the lag time corresponding to this peak value is 0 day. Also, these cross correlation plots are almost symmetr ical and centered at lag time=0. This indicates that the groundwater levels at these four observation wells are highly correlated and they react at the same time to a factor, which might be rainfall in our analysis. PAGE 132 134 (a) V1091~V0083 (1/1/2004 to 8/16/2004) (b) V0867 ~ V0083 (1/1/2004 to 8/16/2004) (c) L0059 ~ V0083 (5/ 24/2005 to12/31/2006) Figure 618. Cross correlati on of groundwater levels at the three observation wells with the well V0083 PAGE 133 135 The cross correlation analysis suggests groundwater levels at the four observation wells react at the same time. If the response is to rainfall, then cross correlations of rainfall with gr oundwater levels at different wells should be consistent. The longest continuous records of wate r level data exist for wells V1091 and V0867 for the period extending from 1/1/2004 to 6/30//2006 (Figur e 63). Groundwater levels data at these two wells were used to analyze the cross correlation of rainfall and groundwater levels. Figure 619 shows cross correlation of rainfall and groundwater levels at V1091 and V0867 fo r the period from 1/1/2004 through 6/30/2006. Figure 619. CCF of rainfall versus GW leve ls for the period 1/1/2004 to 6/30/2006 For well V1091, the peak value of the CCF is 0.379, for a 72 day lag. And for V0867, the CFF peaks at 0.443, for a lag equal to 74 days. Other than the obvious peak value at 72 days for well V1091, there is another small increase to 0.2338, for a lag equal to 146 days. PAGE 134 136 Cross correlations between rainfall and groundwater levels were also examined for each year. Results are s hown in Figure 620 and Figure 621. In the year 2004, the peak value of t he cross correlation function at V1091 is 0.524, for a lag equal to 73 days. And at 75 days there is a CCF peak of 0.557at V0867. In the analysis of the year 2005, the di fference between cross correlation of rainfall with groundwater levels at the wells V1091 and V0867 is bigger than other periods of analysis. The highest value of t he CCF for rainfall versus groundwater levels at V1019 is 0.423, for a lag equal to 108 days. The CCF at V0867 peaks at 0.469 after 69 days. Cross correlation analysis of rainfall with groundwater levels for the wet seasons in 2004 and 2005 was also performed. The analysis is shown in Figure 622 and Figure 623. Two CCF peaks are evident for the 2004 wet season (Figure 622). The first occurs at one day for a peak value of 0.479 at V1091 and 0. 447 at V0867. The second peak CCF value is 0.418 at V1091 and 0.455 at V0867; both at day 8. For the wet season in 2005 (Figure 623), the cross correlation of rainfall with groundwater levels at V1091 and V0867 ar e very different. No obvious peak appears within the maximum lag time tested of 20 days. For the period exending from 5/24/2005 th rough 6/30/2006, there is continuous data for all of the f our observation wells (Figure 63). Figure 624 shows the cross correlations of rainfall with groundwater leve ls data at the four observation wells from 5/24/2005 to 6/30/2006. The maximum lag tested is 200 days. PAGE 135 137 (a) (b) Figure 620. Cross correlation of rainfall wi th GW levels for t he year 2004. (a) Daily data of 2004 (total: 365), (b) CCF of rainfall and GW levels in 2004. (a) (b) Figure 621. CCF of rainfall with GW levels for the year 2005. (a) Daily data of 2005 (total: 365), (b) CCF of rainfall and GW levels in 2005 PAGE 136 138 (a) (b) Figure 622. Cross correlati on analysis of rainfall with GW levels from 6/3/2004 to 10/15/2004 (wet season in 2004). (a) Da ily data of the analyzed period (total: 135 days), (b) CCF of the analyzed period (a) (b) Figure 623. Cross correlati on analysis of rainfall with GW levels from 5/30/2005 to 11/2/2005 (wet season in 2005). (a) Daily data of the analyzed period (total: 135 days), (b) CCF of the analyzed period PAGE 137 139 Figure 624. CCF or rainfall with GW at four observation wells (5/24/2005 to 6/30/2006) It is apparent that the cross correlations of rainfall versus groundwater levels at different wells are very similar. From the CCF plots, it is apparent cross correlation values rise to a high value and plateau for 5075 days. Consequently, it is difficult to designate a peak CCF value. However, there appears to be a peak value after 147 days. 6.2.4 Correlation Analysis of Gr oundw ater Levels versus Blue Spring Discharge Figure 625 plots all available data on gr oundwater levels at four observation wells and Blue Spring discharge in a same figure. An examination of the cross correlation of groundwater levels at f our observation wells with Blue Spring discharge for different periods is presented. PAGE 138 140 Figure 625. Available data of sp ring discharge and groundwater levels Period 1: 1/1/2004 ~ 7/25/2004: There are 207 data in Period 1(Figure 626(a)), and the maximum lag time tested is 150 days. For this period, there are no data of groundwater levels at the we ll L0059. Figure 626(b) shows cross correlation of groundwater levels at th ree wells V1091, V0867, and V0083, with Blue Spring discharge. (a) (b) Figure 626. Cross correlation analysis of GW levels with Blue Spring discharge from 1/1/2004 to 7/25/2004. (a) Daily data of the analyzed pe riod (total: 207 days), (b) CCF of the analyzed period PAGE 139 141 Period 2: 7/26/2004 ~ 9/7/2004: Period 2 corresponds to the 2004 wet season. There are 44 data in this period (Figure 627a) and no available gr oundwater data at the wells V0083 and L0059. Figure 627(b) shows the cross correlation of groundwater levels at V1091 and V0867 with Blue Spring discharge. Period 3: 5/24/ 2005 ~ 12/31/2005 : Period 3 includes 233 data (Figure 628(a)), and the maximum lag time tested is 100 days. Figure 628(b) shows cross correlations of groundwater le vels at the four observati on wells with Blue Spring discharge. Period 4: 1/1/2006 ~ 8/21/2006: There are 233 data in this Period 4 (Figure 629(a)), and the maximum lag time tested is 100 days. Figure 629(b) shows cross correlation of groundwater leve ls at the four observation wells with Blue Spring discharge. (a) (b) Figure 627. Cross correlation analysis of GW levels with Blue Spring discharge from 7/26/2004 to 9/7/2004. (a) Daily data of the analyzed pe riod (total: 44 days), (b) CCF of the analyzed period. PAGE 140 142 (a) (b) Figure 628. Cross correlation analysis of GW levels with Blue Spring discharge from 5/24/2005 to 12/31/2005. (a) Daily data of the analyzed period (total: 233 days), (b) CCF of the analyzed period. (a) (b) Figure 629. Cross correlati on analysis of GW levels with Blue Spring discharge from 1/1/2006 to 8/21/2006. (a) Daily data of the analyzed pe riod (total: 233 days) (b) CCF of the analyzed period. PAGE 141 143 From the plots of cross correlation of groundwater levels at four observation wells with Blue spring discharge (Figures 626, 627, 628, 629), one can see that the peak values of cross correlations of gro undwater levels at four observation wells versus Blue Spring discharge for all analyz ed periods are very high. Most peak values are near to 0.8. The common lag ti me corresponding to these peak values is 0 days. 6.2.5 Discussion (1) Cross Correlation of Rainfall w ith Blue Spring Discharge: Table 61 summarizes cross correlation results for rainfall at Deland versus Blue Spring discharge. Table 61. Results of cross correlation analy sis of rainfall at Deland with Blue Spring discharge Analyzed periods Period length (days) Average d rainfall (in/day) 1st lag (days) 1st peak 2nd lag (days) 2nd peak 12/8/2001 ~ 9/7/2004 1005 153 0.198 Two periods 3/21/2005 ~ 8/21/2006 467 58 0.348 159 0.218 1/1/2002 ~ 12/31/2002 365 0.172 104 0.314 1/1/2003 ~ 12/31/2003 365 0.137 35 0.147 154 0.127 1/1/2004 ~ 9/7/2004 251 0.248 No obvious peaks. The CCF values for lags between 130 days and 160 days are around 0.2. 1/1/2005 ~ 12/31/2005 286 0.221 56 0.43 158 0.103 Yearly Analysis 1/1/2006 ~ 6/30/2006 181 0.079 4960 0.375 6/23/2002 ~ 9/26/2002 96 0.326 All the CCF values for the positive test lags are less than 0.1 7/26/2004 ~ 9/7/2004 44 0.813 1 0.465 10 0.168 Wet Seasons Analysis 7/6/2005 ~ 9/10/2005 67 0.308 All the CCF values for the positive test lags are less than 0.1 PAGE 142 144 Results of lag times can be divided into two groups. The first group includes lag times that are a little more than 150 days (153, 159, 154, and 158 days. The second group includes lag times that are a little more than 50 days, which are 58 days, 56 days, 4960 days. These two groups suggest the existence of multiple flow components. Other than the two dominant groups of estimated lag times, there are three other values. The first is 104 days, extr acted from 2002 data a nd corresponding to a CCF equal to 0.314. The second value is 35 days, analyzed from 2003 data and corresponding to a CCF peak of 0.147. The th ird and last estimate of lag was for 1 day, from the wet season data of 2004, and corresponding to a CCF equal to 0.465. For the year 2004, the data from 9/8/ 2004 to 12/31/2004 is missing, so the yearly analysis for 2004 is only based on t he data from 1/1/2004 to 9/7/2004. For the cross correlation analysis for this period, there is no obvious peak value, and the higher CCF of 0.2 extends over a lag interval from 130 to 160 days. This range (130160 days) is similar to the first major group of observed lag times (a little more than 150 days). However, these results are very different from the results of the cross correlation analysis of rainfall with re spect to spring disc harge for the 2004 wet season. This might be due to the rechar ge condition during the wet season, the length of which is 44 days, and very differ ent from the rest of year 2004. The average rainfall for the wet season from 7/26/2004 to 9/7/2004 is 0.813 inches/day and the average rainfall for the rest of 2004 (1/1/2004 to 7/25/2004) is only 0.128 inches/day. PAGE 143 145 From all years 2002, 2003 and 2005, the CCF was the highest in 2005 (peak value 0.43) when the average rainfall was al so the highest (0.221 inches) compared to the other two years. On the other hand, 2003 produced the lowest average rainfall and the lowest CCF peak. This could be ex plained by the different hydrodynamic characteristics of karst aquife rs for different recharge co nditions. If the recharge is low, most event water will be retarded by t he storage of water in the epikarst and the aquifer; giving in turn dampened cross correlation between recharge and spring discharge. As suggest above, less recharge should correspond to lower cross correlations. However in 2006, annual precip itation was also very low (0.079), and yet the CCF for rainfall versus spring discharges peaked at 0.375. This might because the cross correlation between re charge and spring discharges is also influence by antecedent rainfall events. If, for example, antecedent rainfall results in significant recharge and aquifer storage, subsequent precipitation events under these conditions could have a more pronounce affect on spring discharge and therefore increase the cro ss correlation between rainfall and spring discharge. As for the analysis of wet seasons in the year 2002 and 2005, there is no obvious cross correlation between rainfall and spring discharges for lag times less than 50 or 60 days. This might because it took more than 60 da ys for most rainfall events, during the wet seasons for these two years, to influence the Blue Spring discharge. For the wet seasons in the year 2004, there are two peak values of CCF. The first peak value is 0.465, which is rela tively high, for a lag time equal to 1 day; the second is 0.168, after 10 days. The high cross correlation between rainfall and PAGE 144 146 spring discharges for the wet season in the year 2004 might because of extremely high rainfall volumes (flood condition). During this kind of flood condition, some parts of the karst aquifer, such as the epik arst, might be flushed with water, and the rainfalls will travel to the spring by direct and fast channel flow, which is a different process than the normal recharge condition. The lag time of one day detected in the analysis of flood season in 2004 and long lag times at other times imply t hat there might be we ll developed conduit network, while, is delayed by epikarst at most time. (2) Cross Correlation of Rain fall with Groundwater Levels: Table 62 summarized results of cross correlations of rainfall versus groundwater levels at the four observation wells: V1091, V0867, V0083, L0059. Table 62. Results of cross correlation ana lysis of rainfall at Deland and groundwater levels Analyzed periods Period length (days) Averaged rainfall (in/day) Wells 1st lag (days) 1st peak 2nd lag (days) 2nd peak V1091 147 0.392 V0867 147 0.370 V0083 147 0.380 Period for four wells 5/24/2006 ~ 6/30/2006 403 0.172 L0059 Rise to a high value and keep level for a while About 0.47 147 0.407 V1091 72 0.379 146 0.234 The whole period 1/1/2004 ~ 6/30/2006 912 V0867 74 0.443 V1091 73 0.524 1/1/2004 ~ 12/31/2004 366 0.205 V0867 75 0.557 V1091 108 Yearly Analysis 1/1/2005~ 12/31/2005 365 0.188 V0867 69 V1091 1 0.479 8 0.418 7/26/2004~ 9/7/2004 44 0.813 V0867 1 0.447 8 0.455 Wet Season Analysis 7/6/2005 ~9/10/2005 68 0.303 No peak values of CCF in this period. PAGE 145 147 Compared to Blue Spring discharge data, groundwater data exhibit a dampened fluctuation. Groundwat er levels between the four observation wells are highly correlated. Peak values of cro ss correlation between wells are nearly 0.99, for lags equal to 0 days. Also, plots of cross correlation functions of groundwater levels at different wells are symmetrical and centered at k=0. This indicates that groundwater levels at these four observation wells are highly correlated and they react at the same time to a factor, which in ou r analysis might be rainfall. Cross correlations of rainfall versus gr oundwater levels at different wells are very similar for most analysis periods. T here is continuous data of groundwater levels at V1091 and V0867 from 1/1/2004 through 6/30/2006. For this period, the response time 74 days between a rainfall ev ent and groundwater at these two wells. This result is same as the analysis for the year 2004. From the cross correlation analysis of data from 2005, resu lts are different between wells V1091 and V0867. The lag time at V0867 is close to the 2004 result; whereas for V1091 the lag time from 2005 is much longer than from 2004. One possible explanation is that storativity is higher in the area around the well V1091. During 2004, there was abov e average rainfall, leavin g groundwater levels at elevated levels in the area around V1091 and in turn more responsive to rainfall events. In 2005, the rainfall was satisf ying the storage capacity of the aquifer which delayed aquifer responses to the effects of rainfall. In the wet season of 2004, there was c onsiderable rainfall. The lag time of rainfall with groundwater levels at the wells V1091 and V0867 were very short. During this flood season, it is speculated that aquifer conduits were filled with PAGE 146 148 recharge, and the groundwater in these conduits flowed to surrounding areas, therefore the groundwater leve ls responded quickly rainfall. During the wet season of 2005, the rainfa ll was much less. And there were no peaks in the CCR at wells V1091 and V0867 fo r the longest lag tested of 20 days. This means rainfall has little influence on groundwater levels at these two wells within an interval of 20 days during this period. In addition, the plot of cross correlations of rainfall versus groundwater levels at V1091 and groundwater levels at V0867 are very different, for reasons explained previously. The analysis of wet season in 2004 is ve ry different from the analysis results from the 2005 wet season. This might due to differences in flow regime between the two periods. After the above average wet season in 2004, karst channels were filled with water, and groundwater le vel respond quickly to changes in conduit flow. During the wet season in 2005, there is still more rainfall than the dry season, but no flood condition, consequently groundwater levels are under the influenced of slow matrix flow. For the cross correlation analysis of the period fr om 5/24/2005 through 6/30/2006, cross correlations of rainfall versus groundwater levels at the four different observation wells are very similar. Cross correlations rose to a high value and plateau for ~75 days making it difficult to assess a response time. However, there is a peak in the CCF plot for a lag equal to 147 days. This lag time is near the longer lag time analyzed from the cross correlation analysis of rainfall with Blue Spring discharge, which is 158 days. PAGE 147 149 (3) Cross Correlation of Groundwater Levels with Blue Spring Discharge: Cross correlations of groundwater levels at different wells with Blue Spring discharge are very similar. Groundwater levels at the four observation wells (V1091, V0867, V0083 and L0059) are highly co rrelated. Most peak values of cross correlation between groundwater levels at the four observation wells and Blue Spring discharge are near to 0.8 (Figures 327, 328, 329, 330). The as sociated lag times are equal to 0 days. One possible explanation is t hat the groundwater le vels at the four observation wells and Blue Spring discharge do not influence each other. They are controlled by different systems. 6.3 Calibration and Validation of th e Geomorphologic Karst Model of B lue Spring In this section, a geomorphologic kars t model is calibrated and validated to simulate Blue Spring discharge. From prev ious analyses, it can be seen that the geomorphologic karst model using common conduits for diffe rent flow components (Figure 51) is the most effective. The equation for this model is equation (58). To apply this model to Blue Spring area, specific information is needed. 1) The recharge area of the basin: Ve cchioli et al. (1990) estimated the approximate area of the Blue Spring basin to be 268 mi2, of which 46 mi2 was deemed to be discharge area with artesian fl ow. They calculated the approximate area of internally drained terrain in the spring basin to be 121 mi2 (Green 2008). Hence, for this analysis, the recharge ar ea of Blue Spring aquifer is estimated as 121 mi2, with rainfall distributed uniformly over the basin. 2) The flow travel lengt h of the highest order systemL :Unlike surface basins or the hypothetical aquifer, the length of karst conduits are not easily measured. PAGE 148 150 These lengths must be estimated. The distance between Deland and Blue Spring is approximately 5.6 miles. This value is us ed as the characteristic length of the conduit of the highest order in the equation. 3) Hortons numbers R Band R L: As explained above, Hortons numbers RB and RL for conduits cannot be determined dire ctly as is done for surface basins. These parameters need to be calibrated using observed data. 4) The mean length of the porous mediumLmwill be calibrated using observed data 5) The mean flow velocity of gr oundwater flow in porous medium Vm, and mean flow velocity in conduits Vc will be calibrated using observed data. 6) The effective rainfall ratio rt:Rainfall at Deland is used as the only source of recharge to the aquifer, but not all rainfall enters the aquifer. The ratio of recharge to the aquifer to total precipitation is t he effective ratio. This parameter will be calibrated using observed data. The observed Blue spring discharge is s hown in Figure 62. The data from 12/8/2001 to 4/14/2004 includes total of 859 days of data used to calibrate the model. The data from 4/15/2004 to 9/7/2004 was not used to calibrate parameters. This data was omitted because the equipment was damaged by hurricanes in September 2004, therefore t he quality of data during that time is suspect. Using the above data it was possible to calibrate the geomorphologic model. Rainfall data from 12/25/1986 to 4/14/2004 was used to generate spring discharge. The calibrated parameter values were: The mean length of the porous medium: mile Lm 0.3 PAGE 149 151 The velocity of the porous medium flow: daymile Vm /0027.0 The ratio of rainfall to the slow component: mile rt 3069.01 Then, the ratio of rainfall to the conduit flow: mile rt 017.01 As expected, the calibrated results for the R L and Vc are not unique. Table 63 shows the calibrated results of cVfor different values of R L. Figure 630 depicts a plot of data presented in Table 63. Figure 631 displays model simulated discharges of measured Blue Spring discharges from 12/8/2001 to 4/14/2004. The average relative error is 0.0228, which is acceptable. Table 63 Calibrated results of Vc for different R L R L Vc R L Vc R L Vc 3.5 0.03 2.8 0.0327 2.1 0.0381 3.4 0.0304 2.7 0.0333 2.05 0.0387 3.3 0.0307 2.6 0.034 1.9 0.0404 3.2 0.031 2.5 0.0346 1.8 0.0419 3.1 0.0314 2.4 0.0353 1.7 0.0434 3.0 0.0317 2.3 0.0361 1.6 0.0452 2.9 0.0322 2.2 0.0371 1.5 0.0473 To validate the calibrated geomorphologic kars t model, it was used to predict spring discharges from 3/21/2005 to 6/30/2006 for direct comparison to measured Blue Spring discharges. For this simulation De land rainfall data from 12/25/1986 to 6/30/2006 were used to generate the spring discharge. Figure 632 shows the comparison. The average relative error is 0. 0441. This error is tolerable, and it can be said that the calibrated geomorphologic karst model can be used to simulate Blue Spring discharge. PAGE 150 152 Figure 630. Plot of Vc versus RL Figure 631. Comparison of observed spring discharge with results from the calibrated geomorphologic karst model in the calibration period PAGE 151 153 Figure 632. Comparison of observed spring discharge with results from the calibrated geomorphologic karst m odel in the validation period 6.4 Analysis of the Travel Time of Rainfall in the Geomorphologic Models As shown in the Figure 51, in a geomorphologic karst model with same conduits network for different flow com ponents, the fast com ponent is the conduits flow. The slow component is the combinatio n of the flow through the porous to reach the conduits network and the conduits flow. Then the mean travel time of slow component can be calculated as the summation of the rainfall spends in porous medium with the mean travel time in the c onduits network. The mean travel time in porous medium can simply calculated by t he calibrated values of Lm and Vm. The mean travel time of the conduits network c an be estimated as similar way in section 4.5. The calculation is shown in Table 64. In Table 64, column (19) is the es timated mean time spent by slow flow component in porous medium to reach the conduits network. Column (20) is the estimated mean travel time of the conduits flow, which is the fast flow component. PAGE 152 154 Column (21) is the mean travel time of the slow component, which is equal to the summation of the column (20) and column (19). 6.5 Application of Laplace Transforms to Estimate the Travel Time The slow flow component is simulated by a fourth order geomorpholog ic model with equation of IUH shown as (326), and th e central normalize mo ment of the IUH can be calculated as equation (413). Substituting expressions (52) and (318) into (413). The central moment of IUH of the slow flow component can be calculated. The fast flow component is simulated by a third order geomorphologic karst model with equation of IUH shown as (54). The central moment of this IUH can be developed as: 2 4 4 4 3 3 2 2 1 3 4 4 2 4 3 2 3 2 2 2 1 ,0 ,1 ,12 2 C C CC C C CC m m ut t t (61) Where, 2 1 21**0 A C (62a) 2 2 4 2 43 3 3 2 22**0 **0 A C (62b) 2 3 4 3 4 4 3 3 3 4 43 3 4 4 23 232**0**0 **0*0 A A C (62c) 2 4 4 3 4 2 43 3 4 3 24*0 **0 **0 A C (62d) PAGE 153 155 Table 64. Analysis mean travel time of different flow components mile/day miles miles miles m iles miles miles miles miles 02 03 04 miles days days days (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17) (18) (19) (20) (21) RL Vc L4 L3 L2 Lm path1 path2 path3 path4 p12 p13 p14 p1 p2 p3 p4 Lc Tm Tc Tt 1.5 0.05 5.60 3.73 2.49 3.00 11. 82 8.09 9.33 5.60 0. 57 0.29 0.14 0.50 0. 08 0.29 0.14 9.94 111 1.11 210.14 1321.25 1.6 0.05 5.60 3.50 2.19 3.00 11. 29 7.79 9.10 5.60 0. 55 0.29 0.16 0.48 0. 07 0.29 0.16 9.50 111 1.11 210.20 1321.31 1.7 0.04 5.60 3.29 1.94 3.00 10. 83 7.54 8.89 5.60 0. 53 0.30 0.17 0.46 0. 07 0.30 0.17 9.13 111 1.11 210.35 1321.46 1.8 0.04 5.60 3.11 1.73 3.00 10. 44 7.33 8.71 5.60 0. 51 0.31 0.18 0.44 0. 07 0.31 0.18 8.81 111 1.11 210.30 1321.41 1.9 0.04 5.60 2.95 1.55 3.00 10. 10 7.15 8.55 5.60 0. 49 0.31 0.20 0.43 0. 07 0.31 0.20 8.54 111 1.11 211.32 1322.43 2.05 0.04 5.60 2.73 1.33 3.00 9. 66 6.93 8.33 5.60 0. 47 0.32 0.21 0.41 0. 06 0.32 0.21 8.22 111 1.11 212.32 1323.43 2.1 0.04 5.60 2.67 1.27 3.00 9. 54 6.87 8.27 5.60 0. 47 0.32 0.22 0.40 0. 06 0.32 0.22 8.11 111 1.11 212.97 1324.08 2.2 0.04 5.60 2.55 1.16 3.00 9. 30 6.76 8.15 5.60 0. 44 0.32 0.24 0.38 0. 06 0.32 0.24 7.90 111 1.11 213.05 1324.16 2.3 0.04 5.60 2.43 1.06 3.00 9. 09 6.66 8.03 5.60 0. 42 0.33 0.25 0.37 0. 06 0.33 0.25 7.74 111 1.11 214.40 1325.51 2.4 0.04 5.60 2.33 0.97 3.00 8. 91 6.57 7.93 5.60 0. 41 0.33 0.26 0.36 0. 05 0.33 0.26 7.59 111 1.11 215.11 1326.22 2.5 0.03 5.60 2.24 0.90 3.00 8. 74 6.50 7.84 5.60 0. 40 0.33 0.27 0.34 0. 05 0.33 0.27 7.46 111 1.11 215.66 1326.78 2.6 0.03 5.60 2.15 0.83 3.00 8. 58 6.43 7.75 5.60 0. 38 0.33 0.29 0.33 0. 05 0.33 0.29 7.34 111 1.11 215.99 1327.10 2.7 0.03 5.60 2.07 0.77 3.00 8. 44 6.37 7.67 5.60 0. 37 0.33 0.30 0.32 0. 05 0.33 0.30 7.24 111 1.11 217.31 1328.42 2.8 0.03 5.60 2.00 0.71 3.00 8. 31 6.31 7.60 5.60 0. 36 0.33 0.31 0.31 0. 05 0.33 0.31 7.14 111 1.11 218.33 1329.44 2.9 0.03 5.60 1.93 0.67 3.00 8. 20 6.27 7.53 5.60 0. 34 0.33 0.32 0.30 0. 05 0.33 0.32 7.05 111 1.11 218.97 1330.08 3 0.03 5.60 1.87 0.62 3.00 8. 09 6.22 7.47 5.60 0. 33 0.33 0.33 0.29 0. 04 0.33 0.33 6.97 111 1.11 219.86 1330.98 3.1 0.03 5.60 1.81 0.58 3.00 7. 99 6.18 7.41 5.60 0. 32 0.33 0.34 0.28 0. 04 0.33 0.34 6.90 111 1.11 219.60 1330.71 3.2 0.03 5.60 1.75 0.55 3.00 7. 90 6.15 7.35 5.60 0. 31 0.33 0.36 0.27 0. 04 0.33 0.36 6.83 111 1.11 220.23 1331.34 3.3 0.03 5.60 1.70 0.51 3.00 7. 81 6.11 7.30 5.60 0. 30 0.33 0.37 0.26 0. 04 0.33 0.37 6.76 111 1.11 220.34 1331.45 3.4 0.03 5.60 1.65 0.48 3.00 7. 73 6.08 7.25 5.60 0. 29 0.33 0.38 0.25 0. 04 0.33 0.38 6.71 111 1.11 220.60 1331.71 .5 0.03 5.60 1.60 0.46 3.00 7. 66 6.06 7.20 5.60 0. 28 0.33 0.39 0.25 0. 04 0.33 0.39 6.65 111 1.11 221.75 1332.86 PAGE 154 156 Substituting expressions (6 2) and (53) into (61). The central moment of IUH of the fast flow component c an be calculated. Table 65 s hows the calculated results. Table 65 Central moments of IUH for t he two flow components calculated from the common conduit network model for Blue Spring RL Vc Central moment of IUH2 (days) Central moment of IUH1 (days) 3.5 0.03 221.73 1332.837 3.4 0.0304 220.58 1331.689 3.3 0.0307 220.31 1331.32 3.2 0.031 220.21 1331.425 3.1 0.0314 219.58 1330.687 3 0.0317 219.84 1330.95 2.9 0.0322 218.94 1330.05 2.8 0.0327 218.30 1329.41 2.7 0.0333 217.28 1328.395 2.6 0.034 215.96 1327.07 2.5 0.0346 215.63 1326.745 2.4 0.0353 215.08 1326.189 2.3 0.0361 214.37 1325.479 2.2 0.0371 213.01 1324.123 2.1 0.0381 212.27 1323.381 2.05 0.0387 211.60 1322.709 1.9 0.0404 211.28 1322.388 1.8 0.0419 210.26 1321.371 1.7 0.0434 210.31 1321.424 1.6 0.0452 210.15 1321.266 1.5 0.0473 210.10 1321.212 PAGE 155 157 CHAPTER 7 CONCLUSIONS AND RECOMMENDATIONS The objective of this study was to develop simple analyt ical models for investigating and simulating groundwater fl ow in karst aquifers using readily available data. Because karst aquifers are extremely heterogeneous, conduits are indeterminate, and recharge is spatially variable, developing a physicallybased numerical model of karst systems is difficu lt at best. In addition, historical applications of time series analysis indicate the importance of considering multiple flow components when developing a karst flow model. In this research, simple analytical geomorphologic karst model we re developed for investigating and simulating groundwater flow in karst aqui fers. These models were tested using synthetic data generated from a numerical Dual Conductivity Model (DCM), and one specific model was successfully calibra ted and validated using measured spring discharges from Blue Spring, Florida. 7.1 Conclusions 1) The geomorphologic karst model wa s developed to simulate spring discharge. The model simulates the process of recharge and flow through karst aquifer as a semiMarkov process, and it ex presses the aquifer response to rainfall as a function of Conduit Horton numbers RL and RB and a mean fl ow velocities porous the medium and through the conduit net work (e.g., slow and fast). Because the conduit network is indeterminate, condui t Horton numbers and velocities must be estimated through model calibration using available rainfall and spring discharge data. PAGE 156 158 2) Successful calibration and validati on of the geomorphologic karst model using data from two hypothetical karst systems and Blue Spring suggest this tool can provide reasonably reliable predi ctions of spring discharge. 3) A comparison of simulation results generated from fourth and fifth order geomorphologic models showed no difference. Hence, fourth order modeling is sufficient. 4) Improved model predictions were obtained when at least two flow components (e.g., fast and slow) were consider ed in parallel. More realistic results were observed when the model was formulated such that the slow flow component shared the same conduits as the fast flow component and flow through the porous medium and the conduit network was treated in series. 5) It was found that the value of RB has little influence on the results. 6) The relationship between RL and v is defined by the first normalized temporal moment derived from the IUH of each flow component. 7) In the common conduit geomorphologic model, the ratio of Lm to Vm is unique. 8) Time series analyses of Blue Spring discharge data and groundwater levels and rainfall indicate: a) significant groundwater storage from autocorrelations of spring discharge and groundwater levels (a sl ow dissipation in autocorrelations over time); b) low cross correlations betw een rainfall and Blue Spring discharge for lag times less 2 months suggest a low degree of karsification; c) on more than one occasion cross correlation peaks were observed at 56 days and 158 days that indicate more than one flow component may exist; and d) cross correlation analysis PAGE 157 159 of rainfall with Blue Spring discharge over diff erent seasons indicate different aquifer responses to different antecedent recharge; e) the lag time of one day detected in the analysis of flood season in 2004 and long lag times at other times imply that there might be well developed c onduit network, while, is delayed by epikarst at most time. 7.2 Recommendations 1) Use a calibrated geomorphic model to create a physically bas e numerical karst flow model where parameter such as RL, Vc, Vm, Lm, rt1 and rt2 are used to define the conduit network and flow c onditions in the numerical model. 2) Develop more complex geomorphologic models that describe multiple nested flow components which share a co mmon conduit network (e.g., slow, fast, and faster flows components). 3) Apply geomorphologic models to c haracterize tracer travel times. PAGE 158 160 APPENDIX A PARAMETERS FOR EQUATION 340 Where, lam1= 1 lam2= 2 ,lam3= 3 ,lam4= 4 ,lam_star= 17a=1 17b=lam_star^2*(lam1^3*p15lam1*p14*lam4*lam2lam1*p15*lam4*lam2lam1*p15*lam4*lam3lam1*p13*lam3*p35*lam2lam1*p13*lam3*p35*lam4lam1*p14*lam4*lam3lam1*p12*lam2* p23*lam3*p35lam1*p12*lam2*p24*lam4lam1*p15*lam2*lam3lam1*p13*lam3*p34*lam4lam1*lam3*p12*lam2*p25lam1*p12*lam2*p25*lam4+lam3*lam2*l am4*p15+lam3*lam2*lam4*p14+lam3* lam2*lam4*p12*p24+lam3*lam2*lam4*p13* p35+lam3*lam2*lam4*p12*p23*p34+ lam3*lam2*lam4*p12*p25+lam3*lam2*la m4*p12*p23*p35+lam3*lam2*lam4*p13*p34 +lam1^2*lam3*p15+lam1^2*p15*lam4+lam 1^2*p15*lam2+lam1^2*p14*lam4+lam1^2 *p12*lam2*p25+lam1^2*p13*lam3*p35)/(la m_starlam1)^2/(lam1lam4)/(lam1lam2)/(lam1+lam3) 17c=(lam3*p25*lam4+lam3*p25*lam2lam 3*lam4*p23*p35lam3*p24*lam4+ lam3*lam2*p23*p35lam3*p23*p34*lam4 +lam2*p25*lam4+lam2*p24*lam4lam2^2*p25)*lam1*p12*lam_star^2/(lam1lam2)/(lam2+lam_star )^2/(lam3lam2)/(lam4+lam2) 17d=(lam3*p35lam4*p35lam4*p34)*(p13*lam3+lam2*p13+lam2*p12*p23)* lam1*lam_star^2/(lam1+lam3)/(lam_starlam3)^2/(lam4+lam3)/(lam3lam2) PAGE 159 161 17e=(lam3*lam2*p14lam3*lam4*p14+lam3*lam2*p13*p34+ lam3*lam2*p12*p23*p34+lam3*p24*lam2*p 12lam3*lam4*p13*p34+lam4^2*p14p12*lam2*p24*lam4p14*lam4*lam2)*lam1*lam_star^2/(lam1lam4)/(lam4+lam3)/(lam_starlam4)^2/(lam4+lam2) 17 f =lam_star*lam1*(lam_star^2*la m2*p15lam_star^2*lam2*p12*p25lam3*lam2*lam4*p15lam3*lam2*la m4*p14lam3*lam2*lam4*p12*p24lam3*lam2*lam4*p13*p35lam3*lam2*lam4*p12*p23*p34lam3*lam2*lam4* p12*p25lam3*lam2*lam4*p12*p23* p35lam3*lam2*lam4*p13*p34+ p15*lam_star^3+lam_star*p12*lam2*p25*la m4+lam_star*lam3*lam4*p13*p35+lam_s tar*p12*lam2*p24*lam4+lam_star*lam3*lam 4*p13*p34+lam_star*lam3*lam2*p12*p2 3*p35+lam_star*lam3*lam2*p13*p35+lam_star*lam3*p25*lam2*p12lam_star^2*p14*lam4lam_star^2* lam3*p15lam_star^2*p15*lam4+ lam_star*p14*lam4*lam2+lam_star*lam 2*lam4*p15+lam_star*lam3*p15*lam4 +lam_star*lam3*lam2*p15+lam_star *lam3*lam4*p14lam_star^2*lam3 *p13*p35)/(lam_starlam1)/(lam2+lam_st ar)/(lam_starlam4) /(lam_starlam3) 17g =lam1*(lam1*lam2^2*lam3^2*la m4^2*p144*lam2^2*lam_star^4* lam4*p15+2*lam2^2*lam_star^3*lam4 ^2*p12*p254*lam2*lam_star^4* lam4*p14*lam1lam4^2*lam3^2*lam 2^2*p13*p34*lam1+2*lam2*lam1* lam_star^5*p154*lam2*lam_star^4* lam4^2*p154*lam2*lam_star^4* lam4^2*p14+2*lam2^2*lam_star^3*p14*la m4^2+8*lam2*lam_star^5*p15*lam4+2*la m2^2*lam_star*lam4^2*p13*p34*lam3^2lam2^2*lam1*lam_star^4*p15+ PAGE 160 162 2*lam2^2*lam_star*lam3^2*lam4^2*p15+2* lam_star*lam1*lam3*lam2^2*lam4^2*p14 +2*lam_star*lam1*lam3^2*lam2*lam4^2*p14 +2*lam2^2*lam_star*lam3^2*lam4^2*p1 2*p25lam3^2*lam2^2*lam4^2*p35*p13*lam1lam2^2*lam1*p15*lam4^2 *lam3^2+2*lam2^2*lam_star*lam4^2* p13*p35*lam3^2+2*lam1*lam2^2*lam4 *lam_star*p13*lam3^2*p35+2*lam1*lam2^2*lam4*lam_star*p12*p23*lam3^2*p35+2*l am1*lam2^2*lam_star*lam3^2*lam4*p15+2*lam1*lam2^2*lam_star*lam3^2*lam4*p1 4+2*lam1*lam2^2*lam_star*lam3^2*lam4*p12*p24+2*lam1*lam2^2*lam_star*lam3^2 *lam4*p12*p23*p34+2*lam1*lam2^2*lam_star *lam3^2*lam4*p12*p25+2*lam1*lam2^ 2*lam_star*lam3^2*lam4*p13*p34+2*lam_ star^3*p13*lam3^2*p34*lam4^2+2*lam1*l am_star^5*lam4*p15+2*lam2^2*lam_star^ 3*lam4^2*p15+2*lam2^2*lam_star^3*p24* p12*lam4^2lam_star^4*lam4^2*p14*lam1lam1*lam_star^4 *lam4^2*p153*lam2*lam_star^4*lam4 ^2*p12*p244*lam2*lam1*lam_star^4* p15*lam4+2*lam2^2*p24*p12*lam1*lam_star^3*lam4lam3^2*lam1*lam2^2 *lam4^2*p12*p24+2*lam1*lam_star* lam3^2*lam2*lam4^2*p12* p25+ 2*lam2*lam_star*lam1*p15*lam4^2*lam3^2+2*lam2^2*lam_star*p12*p24*lam4^2*la m3^2+2*lam1*lam4^2*lam_star*lam3^2*lam2*p13*p35+2*lam1*lam4^2*lam_star*la m3^2*lam2*p13*p34+2*lam1*lam2^2*lam4^2*lam_star*p13*lam3*p35+2*lam1*lam2 ^2*lam4^2*lam_star*p12*p23* lam3*p35+2*lam2^2*lam_star*p14*lam4^2*lam3^2+2*l am1*lam2^2*lam_star*lam3*lam4^2*p12*p2 3*p34+2*lam1*lam2^2*lam_star*lam3*la m4^2*p13*p34+2*lam3^2*lam2*lam4^2*lam_star*p12*p24*lam1+2*lam3^2*lam2*la m4^2*lam_star*p12*p23*p34*lam1+2*lam3^2*lam2*lam4^2*lam_star*p12*p23*p35*l am1+2*lam2^2*lam_star*lam4^2*p12*p23*p 34*lam3^2+2*lam2^2*lam_star*lam4^2* p12*p23*p35*lam3^2+2*lam1*lam2^2*lam4^ 2*p12*p24*lam_star*lam3+2*lam2^2*la PAGE 161 163 m_star*lam1*p15*lam4^2*lam3+2*lam2^2*l am_star*lam1*p12*p25*lam4^2*lam3+2*l am2*lam1*lam_star^3*lam4^2*p15+6*lam3*lam2*lam_star^5*p12*p25+4*lam3*lam2 *p12*p23*p34*lam_star^3*lam4^23*lam3* lam_star^4*lam4^2*p13*p343*lam3^2 *lam2*lam_star^4*p12*p25+2*lam2^2 *lam_star^5*p12*p25lam1*lam2^2*lam4 ^2*p12*p23*lam3^2*p354*lam2*lam4* lam1*lam_star^4*p12*p25lam2^2*p24*p12*lam1*lam_star^2*lam4^2+ 2*lam2*p24*p12*lam1*lam_star^3*lam4^2lam2^2*lam1*lam_star^4*p12*p254*lam2^2*lam_star^4*p12*p25*lam4+2*lam 2*lam1*lam_star^3*lam4^2*p12*p25lam2^2*lam1*p12*p25*lam4^2*lam_star^ 2+2*lam4^2*lam_star^5*p153*lam1* p24*lam2*p12*lam_star^4*lam4+2*lam3*lam 1*lam_star^5*p15+2*lam3^2*lam_star^ 5*p13*p35lam3^2*lam1*lam_star^4* p153*lam_star^2*lam4^2*lam3 *lam2^2*p13*p35+6*lam_star^3*lam4^2*lam3*lam2*p13*p35+6*lam_star^3 *lam4*lam3*lam2^2*p13*p35*lam_star^ 2*lam4^2*lam3*lam2^2*p12*p23*p35+ 4*lam_star^3*lam4^2*lam3*lam2*p12*p23* p35+6*lam_star^3*lam4*lam3*lam2^2 *p12*p23*p358*lam_star^4*lam4*lam3 *lam2*p13*p343*lam_star^2*lam4^2* lam3*lam2^2*p13*p34+4*lam_star^3*lam4* lam3*lam2^2*p13*p34+2*lam_star^3*lam 3^2*lam4^2*p13*p35+2*lam_star^3* lam3^2*lam2^2*p13*p354*lam3^2 *lam2*p13*p35*lam_star^43*lam4^2* p13*p35*lam_star^4*lam3+4*lam2* lam_star^5*lam4*p12*p24+8*lam3^2*l am4*lam2*lam_star^3*p154*lam3* lam_star^4*p12*lam2^2*p2512*lam3*lam_star^4*p14*lam4*lam2+6*lam2^2* lam_star^3*p14*lam4*lam33*lam2^2 *p12*p24*lam4*lam_star^2*lam3^2+ 4*lam_star^3*p12*lam2*p24*lam4*lam 3^28*p24*lam2*p12*lam_star^4*lam4 *lam3+6*lam3*lam_star^5*p13*p35*lam2 +4*lam3*lam_star^5*lam4*p13*p34+ PAGE 162 164 6*lam_star^5*lam3*p14*lam44*lam 3*lam_star^4*p15*lam4^2+8*lam4* lam_star^3*p15*lam2^2*lam3+8*lam4^2*lam_star^3*p15*lam2*lam3+ 6*lam4^2*lam_star^3*lam3*p12*lam2*p25+8*lam2*lam_star^3*p14*lam4^2*lam3+6*l am_star^3*lam3*lam2^2*lam4*p12*p24+6* lam_star^3*lam3*lam2*lam4^2*p12*p24+ 8*lam3*lam4*lam_star^3*p12*lam2^2*p254*lam2^2*p15*lam4^2 *lam3*lam_star^216*lam3*lam2*la m_star^4*p15*lam4+4*lam3*lam2* lam_star^5*p12*p23*p35+8*lam3^2*lam 4*lam2*p13*p35*lam_star^34*p15 *lam4^2*lam3^2*lam_star^2*lam28*lam3*lam_star^4*p12*lam2*p23*p35 *lam45*lam3*lam2*p12*p23*p34*lam_star^4*lam43*lam3*lam_star^4 *lam2^2*p13*p35+6*lam4*lam2*lam_star^5*p12*p254*lam_star^6*p15*lam44*lam2*lam_star^6*p154*lam3*lam_sta r^6*p15+2*lam1*lam_star^5*p14*lam4 +2*lam_star^3*p15*lam4^2*lam3^2+2*lam2^2*lam_star^3*lam3^2*p154*lam3*lam2^2*lam_star^4*p15+8*la m3*lam_star^5*p15*lam44*lam3^2*p15 *lam_star^4*lam44*lam_star^4*p14* lam4^2*lam3+2*lam_star^3*p14* lam4^2*lam3^23*lam_star^4*p14*lam4 *lam3^23*lam2^2*p14*lam_star^4 *lam4+6*lam2*p14*lam_star^5*lam43*lam4^2*p12*lam2*p23*lam3^2* p35*lam_star^24*lam4*p12*lam2^2* p23*lam3^2*p35*lam_star^2lam3^2 *lam2^2*lam1*lam4^2*p12*p25+6*lam_ star^3*p13*lam3*p34*lam4^2*lam2 +6*lam_star^3*p13*lam3^2*p34*lam4*lam2+ 2*lam_star^3*p12*lam2^2*p23*lam3^2* p35+2*lam_star^3*lam3^2*p12*lam2^2*p253*lam3^2*lam2^2*lam4*p12* p23*p34*lam_star^23*lam3^2*la m2*lam4^2*p12*p24*lam_star^2+ 6*lam_star^3*p12*lam2*p23*lam3^2* p35*lam44*lam3^2*lam2^2*lam4* PAGE 163 165 p12*p25*lam_star^23*lam3*lam2^2 *lam4^2*p12*p23*p34*lam_star^23*lam3^2*lam2*lam4^2*p12*p23*p34*lam_ star^2+6*lam_star^3*lam3^2*p12* lam2*p25*lam43*lam3^2*lam2*lam4^2* p12*p25*lam_star^23*lam_star^2* lam4*lam3^2*lam2^2*p13*p344*lam_star^2*lam4^2*lam3^2*lam2*p13*p354*lam2^2*p14*lam4^2*lam_star^2*lam34*lam_star^2*p14*lam4^2*lam3^2 *lam2+6*lam_star^3*p14*lam4*lam3^2*lam23*lam2^2*p14*lam4*lam_star^2 *lam3^24*lam_star^2*lam4*lam3^2 *lam2^2*p13*p354*lam_star^2*lam4^2 *lam3^2*lam2*p13*p34lam2^2*lam1* lam_star^2*p14*lam4^2+2*lam2^2 *p14*lam1*lam_star^3*lam44*lam3*lam1 *lam_star^2*lam2*lam4^2*p144*lam3* lam_star^2*lam2^2*lam1*p14*lam44*l am3*lam_star^2*lam2^2*lam1*p12*p24 *lam4+6*lam3*lam1*lam_st ar^3*p12*lam2*p24*lam44*lam3*lam1*lam_star^4 *p14*lam4+2*lam3*lam1*lam4^2*lam_ star^3*p15+2*lam3*lam1*lam2^2* lam_star^3*p15+2*lam3*lam1*lam_star^3 *p14*lam4^2+2*lam3*lam1*lam_star^5*p1 3*p354*lam3*lam1*lam2*lam_star^4* p154*lam3*lam1*p15*lam_star^4*lam44*lam3*lam1*lam2*lam4^2*p12*p24*lam_ star^2+8*lam3*lam1*lam_star^3*p12* lam2*p25*lam44*lam3*lam1*lam2*lam4 ^2*p12*p25*lam_star^2+8*lam3*lam1* lam_star^3*p14*lam4*lam24*lam3*lam1*lam_star^2*lam4^2*lam2*p13*p343*lam3*lam1*lam_star^2*lam2^2*la m4*p12*p23*p34*lam3*lam1*lam_star^2 *lam2^2*lam4*p13*p34+2*lam3*lam1*lam2 ^2*lam_star^3*p12*p23*p35+8*lam3*lam 1*lam2*lam_star^3*p15*lam44*lam3*lam1*lam2*lam4^2*p15*lam_star^24*lam3*lam1*lam2*p13*p35*lam_star ^4+2*lam3*lam1*lam2^2*p13*p35* lam_star^3+2*lam3*lam1*lam2^2*la m_star^3*p12*p253*lam3*lam1*lam2 *lam4^2*p12*p23*p35*lam_star^2+4*lam 3*lam1*lam2*p12*p23*p34*lam_star^3 PAGE 164 166 *lam44*lam3*lam1*lam2^2*lam4*p12* p23*p35*lam_star^24*lam3*lam1* lam2^2*lam4*p13*p35*lam_star^24*lam 3*lam1*lam2^2*lam4*p15*lam_star^2 +2*lam3*lam1*lam_star^3*p13*p35*la m4^2+8*lam3*lam1*lam2*p13*p35* lam_star^3*lam43*lam3*lam1*lam2 *lam_star^4*p12*p23*p35+6*lam3*lam1* lam_star^3*lam4*p13*p34*lam24*lam3*la m1*lam2*lam4^2*p13*p35*lam_star^23*lam3*lam1*lam2*p12*p23*p34*lam_star ^2*lam4^2+6*lam3*lam1*lam2*lam4 *lam_star^3*p12*p23*p354*lam3*lam1*la m2^2*lam4*p12*p25*lam_star^2+ 2*lam3*lam1*lam_star^3*lam4^2*p13*p344*lam3*lam1*lam_star^4*p13*p35* lam43*lam3*lam1*lam_star^4*lam4*p13* p344*lam3*lam1*lam2*lam_star^4* p12*p254*lam3^2*lam4*p13*p34*lam1*lam_star^2*lam24*lam3^2*lam1*p15* lam_star^2*lam4*lam24*lam3^2*lam1 *p14*lam4*lam2*lam_star^23*lam3^2 *lam1*p12*lam2*p24*lam4*lam_star^24*lam3^2*lam1*p12*lam2*p25*lam4 *lam_star^2lam3^2*lam1*lam_star^2* p14*lam4^2lam3^2*lam1*lam_star^2* p13*p34*lam4^23*lam3^2*lam4*lam1*lam_star^2*lam2*p12*p23*p34+2*lam3^2 *lam1*lam2*lam_star^3*p12*p23*p35+2*lam3 ^2*lam1*lam2*p13*p35*lam_star^3+2*l am3^2*lam1*lam2*lam_star^3*p12* p25lam3^2*lam1*lam2^2*p35*p23* p12*lam_star^24*lam3^2*lam1*lam2* lam4*p12*p23*p35*lam_star^2lam3^2 *lam1*lam2^2*p13*p35*lam_star^2l am3^2*lam1*p15*lam4^2*lam_star^2 +2*lam3^2*lam1*lam_star^3*p14*lam4+2*lam3^2*lam1*lam2*lam_star^3*p15lam3^2*lam1*p13*p35*lam4^2*lam_sta r^2+2*lam3^2*lam1*lam_star^3* p13*p35*lam4+2*lam3^2*lam1*lam_star^3* lam4*p13*p34+4*lam3^2*lam2*p12*p23* p34*lam_star^3*lam4lam3^2*lam1*p13*p35*lam_star^43*lam3^2* lam_star^4*lam4*p13*p34+2*lam3^2*lam1* lam4*lam_star^3*p15lam3^2*lam1* PAGE 165 167 lam2^2*p15*lam_star^2lam3^2*lam1*la m2^2*p25*p12*lam_star^23*lam3^2* lam2*lam_star^4*p12*p23*p354*lam3^2 *lam1*lam2*p13*p35*lam_star^2 *lam4+2*lam1*lam2*lam_st ar^5*p12*p25lam2^2*lam1*p15*lam4^2*lam_star^2+ 2*lam2^2*lam1*lam_star^3*lam4*p154*lam2^2*p15*lam3^2*lam_star^2*lam43*lam2*lam_star^6*p12*p25+8*lam3*lam_s tar^5*p15*lam23*lam3*lam_star^6 *p13*p354*lam3^2*lam2*lam_star^4* p15+6*lam4*lam_star^5*p13*lam3*p354*lam4^2*lam_star^2*lam3*p12*lam2^2 *p254*lam_star^2*lam3*lam2^2* lam4^2*p12*p2412*lam3*lam2*lam_st ar^4*lam4*p13*p354*lam3^2*p13*p35 *lam4*lam_star^43*lam2*lam_star^4* p12*p25*lam4^23*lam_star^6*p14*lam4 +2*lam_star^5*p14*lam4^2+2*lam2^2*la m_star^5*p15lam1*lam_star^6*p15+ 2*lam3^2*lam_star^5*p1512*lam3*lam_star^4*p12*lam2*p25*lam4 +2*lam_star^7*p15lam3^2*lam2^2*p12* p23*p34*lam1*lam4^2+2*lam2*p14* lam1*lam_star^3*lam4^2+2*lam2^2*lam1*lam_star^3*lam4*p12*p253*lam3* lam2^2*p12*p23*p35*lam_star^43*lam2^ 2*lam_star^4*lam4*p12*p24+4*lam3 *lam2^2*p12*p23*p34*lam_star^3 *lam4)/(lam_starlam1)^2/ (lam_starlam3)^2/(lam_starl am4)^2/(lam2+lam_star)^2 27a =p23*p35+p25+p23*p34+p24 27b=0 27c=(lam2*p24*lam4lam3*lam4*p24l am3*lam4*p23*p34+lam2*p25*lam4lam3*lam4*p25lam3*lam4*p23*p35+ lam2*p25*lam3+lam2*p23*lam3*p35lam2^2*p25)*lam_star^2/(lam4lam2)/(lam_starlam2)^2/(lam2lam3) PAGE 166 16827d=(lam3*p35+lam4*p35+lam4*p34)*l am_star^2*p23*lam2/(lam_starlam3)^2/(lam4lam3)/(lam2lam3) 27e=lam2*lam_star^2*(p24*lam4lam3*p24p34* lam3*p23)/(lam_starlam4)^2/(lam4lam3)/(lam4lam2) 27 f =lam_star*lam2*(p25*lam_star^2lam _star*p25*lam4lam_star*p24*lam4lam3*lam_star*p23*p35am3*lam_star* p25+lam3*lam4*p23*p35+lam3*lam4*p24 +lam3*lam4*p25+lam3*lam4* p23*p34)/(lam_starlam3) /(lam_starlam2)/(lam_starlam4) 27g=(4*lam2*p23*lam3*p35*lam_star^2*l am43*lam2*p23*lam3*p34*lam4 *lam_star^24*lam2*p24*lam4*lam_star^2*lam34*lam2*p25*lam_star^2*lam4 *lam3lam2*lam_star^2*lam4^2*p24lam2*lam3^2*lam4^2*p24lam2*p25* lam_star^2*lam4^2+2*lam3^2*p24*lam 4^2*lam_star3*lam3^2*p24*lam4* lam_star^24*lam3*p24*lam_star^2*lam 4^2lam3^2*p25*lam2*lam_star^2lam3^2*p25*lam2*lam4^2+2*lam3^2*p25* lam4^2*lam_star4*lam3^2*p25*lam4 *lam_star^24*lam3*p25*lam_star^2*l am4^2lam2*lam3^2*lam4^2*p34*p23+ 2*lam2*lam3*lam4^2*lam_star*p24+2*lam2 *lam3^2*lam4*lam_star*p24+2*lam2*lam 3^2*lam4*lam_star*p34*p23+2*lam2*lam4^ 2*lam_star*p34*lam3*p23+2*lam2*p25*l am4^2*lam3*lam_star+2*p34*lam3^2*p23*lam4^2*lam_star3*p34*lam3^2 *p23*lam4*lam_star^23*p34*lam3*p23* lam_star^2*lam4^2p35*lam3^2*p23 *lam2*lam_star^2p35*lam3^2*p23*lam 2*lam4^2+2*p35*lam3*p23*lam2*lam4^2 PAGE 167 169 *lam_star+2*p35*lam3^2*p23*lam2*lam4*lam_star+2*p35*lam3^2*p23*lam4^2 *lam_star4*p35*lam3^2*p23*lam4*lam_star^23*p35*lam3*p23*lam_star^2 *lam4^2+2*lam3^2*p25*lam2*lam4*lam_s tarlam_star^4*lam2*p25+2*lam_star^3 *lam3^2*p23*p35+2*lam_star ^3*lam3*lam2*p25+2*lam_star^3*lam3*lam2*p23* p35+2*lam_star^3*lam3^2*p254*lam_st ar^4*lam3*p25+2*lam_star^3*p25* lam4^2+2*lam_star^3*p24*lam4^24*lam 4*lam_star^4*p25+6*lam_star^3*p24 *lam4*lam3+6*lam_star^3*p23*lam3*p35* lam4+2*lam_star^5*p253*lam_star^4 *lam3*p23*p35+2*lam4*lam_star^3*lam2 *p24+2*lam4*lam_star^3*lam2*p25 +8*lam4*lam_star^3*lam3*p253*lam4*lam _star^4*p24+4*lam4*lam3*lam_star^3 *p23*p34)*lam2/(lam_starlam4)^2/(lam_ starlam2)^2/(lam_starlam3)^2 37a=p35+p34 37b=0 37c=0 37d=(lam4*p35lam3*p35+lam4*p34)*lam_star^2/(lam4lam3)/(lam3+ lam_star)^2 37e =lam3*lam_star^2*p34/(lam_starlam4)^2/(lam4lam3) 37 f =(lam_star*p35lam4*p35lam4*p34)*l am3*lam_star/(lam_starlam4)/(lam3+lam_star) PAGE 168 17037g=lam3*(2*p35*lam4^2*lam_star+2*p34*lam4^2*lam_starlam3*p34*lam4^2p35*lam3*lam4^23*lam4*lam_star^2*p344*lam4*lam_star^2*p35+2*lam3* lam4*lam_star*p35+2*lam3*lam4*la m_star*p34lam3*lam_star^2*p35+ 2*lam_star^3*p35)/(lam_starl am4)^2/(lam3+lam_star)^2 47a=1 47b=0 47c=0 47d=0 47e=lam_star^2/(lam_starlam4)^2 47 f =lam4*lam_star/(lam_starlam4) 47g=(2*lam_starlam4)*lam4/(lam_starlam4)^2 57a=1; 57b=0; 57c =0; 57d =0; 57e =0; 57 f =lam_star; 57g =1 67a=1; 67b=0; 67c=0; 67d=0; 67e=0; 67 f =0; 67g=1 77a=1; 77b =0; 77c =0; 77d=0; 77e=0; 77 f =0; 77g=0 PAGE 169 171 LIST OF REFERENCES Agyei, E., Hatfield, K. 2005. Enhancing grad ientbased parameter estimation with an evolutionary approac. J. of Hydrol. 316, 266280 A.J. Long, R. G. Derickson, 2002. Linear systems analysis in a karst aquifer. J. of Hydrol. 219, 206217. Armstrong, B., Chan, D., 2003. Doline and a quifer characteristics within Hernando, Pasco, and northern Hillsborough counties. Karst Studies in West Central Florida. Chatfield, C., 2003. The analysis of time series: an introduction. Duan, Q., Sorooshian, S., Gupta,V.K., 1994. Optimal us e of the SCEUA global optimization method for calibrating wa tershed models. J. Hydrol. 158, 265 284. Duan, Q., Sorooshian, S., Gupt a, V., 1992. Effective and effi cient global optimization for conceptual rainfallrunoff models. Water Resour. Res. 28, 1015 Duan, Q., Gupta, V.K.,Sorooshian, S, 1993. Shuffled Complex evolution approach for effective and efficient global minimization J.Op timization Theor., 76, 501 521. DenicJukic, V., Jukic. D., 2003. Composite tr ansfer functions for karst aquifers. J. of Hydrol. 274, 8094. Feryberg, D. L., 1988. An ex ercise in ground water model calibration and prediction. Ground Water. 26, 350360. Green, R.T., 2008. Karst aquifer mode ling in northcentral florida. Hao, Y., Yeh, T.J., Gao, Z., Wang, Y., Zhao, Y., 2006. A gray system model for studying the response to climate change: The Liulin karst springs, China. J. of Hydrol., 328, 668676. Hao, Y., J, Yeh, T.J., Wang, Y., Zhao, Y., 2 007. Analysis of karst aquifer spring flows with a gray system decompositi on model. Ground water. 45, 4652. Larocque, M., Banton, O., Razack, M., 2000. Transientstate history matching of a karst aquifer ground water flow model Ground Water. 38, 939946. Larocque, M., Mangin ,A., Razack, M., Bant on, O., 1998. Contribution of correlation and spectral analyses to the regional st udy of a large karst aquifer (Charente, France). J. of Hydrol., 205, 217231. PAGE 170 172 Legates, D. R., 2002. Evaluati ng the use of goodnessoffit means in hydrologic and hydroclimatic model validatio n. J. of Hydrol. 261, 132149. Littlefield, J.R., Culbreth, M. A.,1984. Relationship of modern sinkhole development to large scalephotolinear features. In: Beck BF (ed) Sinkholes: their geology,engineering and env ironmental impact. Proceedings of the 1st multidisciplinary conference on sinkholes. 1984, A.A.Balkema, Rotterdam,189 Madsen, H. Kristensen, M. A multiobjective calibration framework for parameter estimation in the MIKE SHE integrated hydrologica l modeling system. ModelCARE 2002, Proceedings of the 4th Internationa l Conference on Calibration and Reliability in Gr oundwater Modelling, Prague, Czech. Massei, N., Dupont, J.P., 2006. Investigat ing transport properties and turbidity dynamics of a karst aquifer using correlation, spectral, and wavelet analyses. J. of Hydrol. 329, 244 257 Newman, M. A. 2001. Inverse characterizati on of subsurface contaminant mass flux: a threedimensional physical and numer ical modeling study Dissertation, University of Florida. Obarti, F. J., Garay, P., Morell, I., 1988. An attempt to karst cla ssification in spain based on system analysis Padilla, A., PulidoBosch, A., 1994. Study of hydrographs of karstic aquifers by means of correlation and crossspectral analysis. J. of Hydrol. 168, 7389. Painter, S.T., Sun, A., Green, R.T., 2004. Enhanced characterization and representation of flow through karst aquifers. Palmer, A.N. 1991. Origin and morphology of limestone caves. Geological Society of America Bulletin. (103), 121. Panagopoulos, G., Lambrakis, N. 2006. The c ontribution of time series analysis to the study of the hydrodynamic characteristics of the karst systems: Application on two typical karst aquifers of Greece. J. of Hydrol. 329, 368376 Pinault, J.L., Plagnes, V ., 2001. Inverse modeling of the hydrological and hydrochemical behavior of hydrosystem s: Characterization of karst system functioning. Water Resour. Res.37, 21912204. Rahnemaei, M., Zare, A. R., Nematollahi, H. 2005. Application of spectral analysis of daily water level and spring di scharge hydrographs data for comparing physical characteristics of karstic aquifers. J. of Hydrol. 311, 106. Rodrigueziturbe, I., Valdes,J.B., 1979. T he geomorphologic struct ure of hydrologic response. Water Resour. Res. 15, 14091420. PAGE 171 173 Rutledge, A. T., 1985. Groundwater hydrol ogy of Volusia county, Florida, with emphasis on occurance and movement of brackish water. U.S. GEOLOGICAL SURVEY. WaterResource s Investigations Report, 844206. Samani, N.,2001. Response of karst aquife rs to rainfall and evaporation, Maharlu Basin, Iran. J. of Cave and Karst Studies. 63, 3340. Sepulveda, N., 2001. Comparisons among gr ound water flow models and analysis of discrepancies in simulated transmissiviti es of the upper floridan aquifer in ground water flow model overlap areas. Shoemaker, W.B., OReilly, A.M., Sepulveda, N., Williams, S.A., Mo tz, L.H., Sun, Q., 2004. Comparison of estimated areas c ontributing recharge to selected springs in northcentral Florida by us ing multiple groundwater flow models. U.S. Geological Survey Openfile report 03228 Solomatine, D.P., Dibike, Y.B., 1999. Au tomatic calibration of groundwater models using global optimization techniques. Hydrol. Sci. j. 44(6), 879894. Trivedi, H.V., Singh, J.K., 2005. Applicat ion of Gray system theory in the development of a runoff prediction model. Biosystems Engineering, 92(4) 521526. Worthington, S. 1999. A comp rehensive strategy for under standing flow in carbonate aquifers. Speleogenesis and Evolution of Karst Aquifers, 30 37. Worthington, S.R.H., 2007. Groundwater residence times in unconfined carbonate aquifers. Journal of Cave and Karst Studies, 69 (1) 94102. PAGE 172 174 BIOGRAPHICAL SKETCH Lili Yu was born in Tianjin, China. She earned her Bachelor of Engineering (B.E.) in July 2002 from Hohai University in Nanjing, China, majoring in Water Conservancy and Hydropower Architectura l Engineering. She then completed her M.E. in Hydraulics and River Dynamics at Tianjin Un iversity, Tianjin, China, graduated in April 2005. In August 2005, she began the Ph.D. program in the Civil and Coastal Engineering Department at Univer sity of Florida, under the tutelage of Professor Kirk Hatfield. While at UF, she gained practi cal research experience while working on the project of investig ating groundwater flow in the Blue Spring area. 