PAGE 1
NUMERICAL AND ANALYTICAL MODELING OF DOUBLEPOROSITY AQUIFERS by John A. Powers A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in the Department of Geology in the University of South Florida April, 1992 Major Professor: Mark T. Stewart, Ph.D
PAGE 2
Graduate Council University of South Florida Tampa, Florida CERTIFICATE OF APPROVAL MASTER'S THESIS This is to certify that the Master's Thesis of John A. Powers with a major in Geology has been approved by the Examining Committee on December 12, 1991 as satisfactory for the thesis requirement for the Master of Science degree. Thesis Committee: Major s 'tewart, Ph.D Member: H. Len Vacher, Ph.D Ph.D<
PAGE 3
DEDICATION For Kristen and Sean. ii
PAGE 4
ACKNOWLEDGEMENTS Most importantly I wish to thank my wife, Kristen and my son, Sean. Kristen gave me her support, love and encouragement throughout my graduateschool education, while Sean made me smile and laugh as only a toddler can. Many thanks are due to my parents, for fostering in me a reverence for books and a love of learning that I have never forgotten. Much of what I have achieved in my short life I owe to them. I would like to thank Dr. Mark T. Stewart for his guidance and patience throughout this project. Dr. H. Len Vacher deserves thanks for his insightful review of this manuscript. Dr. Sam Upchurch is to be commended for his help on this project. Thanks are also extended to John R. Fuller. iii
PAGE 5
TABLE OF CONTENTS LIST OF TABLES LIST OF FIGURES ABSTRACT INTRODUCTION . . . . General Objectives . PREVIOUS INVESTIGATIONS . . Field Studies . . . . Modeling Techniques: Fractured Media ... Development of DoublePorosity Concept Analytical Models . . . . Numerical Models . METHODS Description of Numerical Code . Development of Numerical Model . . Boundary Conditions . . . Model Parameters . . . . RESULTS . . . SteadyState Model . Transient Model DISCUSSION . SteadyState Model Transient Model CONCLUSIONS . LIST OF REFERENCES iv Page v vi ix 1 1 4 5 5 6 7 8 11 14 14 14 20 20 23 23 32 47 47 49 59 61
PAGE 6
Table 1 Model Parameters LIST OF TABLES v Page 21
PAGE 7
LIST OF FIGURES Figure 1 Boundary conditions and grid discretization used in the model. The model is unconfined and consists of one layer with noflow boundaries on all sides and bottom. The model is 5,350 feet on each side and 100 feet thick. A vertical exaggeration of approximately lOX is used here. The arrow indicates direction of groundwater flow. The spring node is Page defined with a constant head of zero. 19 2 Equipotentials (feet) from the steadystate model utilizing a 1:10 block/fracture ratio of hydraulic conductivity. 24 3 Equipotentials (feet) from the homogeneous steadystate model. 25 4 Equipotentials (feet) from the steadystate model utilizing a 1:100 block/fracture ratio of hydraulic conductivity. 26 5 Equipotentials (feet) from the steadystate model utilizing a 1:500 block/fracture ratio of hydraulic conductivity. 27 6 Fluidvelocity vectors and equipotentials (feet) from the homogeneous, steadystate model. 28 7 Fluidvelocity vectors and equipotentials (feet) from the steadystate model utilizing a 1:10 block/fracture ratio of hydraulic conductivity. Fractures are shown. 29 8 Fluidvelocity vectors and equipotentials (feet) from the steadystate model utilizing a 1:100 block/fracture ratio of hydraulic conductivity. Fractures are shown. 30 vi
PAGE 8
9 Fluidvelocity vectors and equipotentials (feet) from the steadystate model utilizing a 1:500 block/fracture ratio of hydraulic conductivity. 31 10 Comparison of timedrawdown data from the homogeneous, transient model with a Theis (1935) type curve. Simulated values of T and S are 1.54E+006 ft2jday and 0.1, respectively. 33 11 Comparison of timedrawdown data from the transient model with a Theis (1935) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at early time. Simulated values of T for blocks and fractures are 1.54E+006 ft2jday and 1.54E+007 ft2jday respectively, and for s, 0.1 and 0.01, respectively. 34 12 Comparison of timedrawdown data from the transient model with a Theis (1935) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at late time. Simulated values of T for blocks and fractures are 1.54E+006 ft2jday and 1.54E+007 ft2jday respectively, and for s, 0.1 and 0.01, respectively. 35 13 Comparison of timedrawdown data from the transient model with a Theis (1935) type curve. A 1:500 block/fracture ratio of hydraulic conductivity is utilized. Match is at early time. Simulated values of T for blocks and fractures are 1.54E+006 ft2jday and 7.70E+008 ft2jday respectively, and for s, 0.1 and 0.01, respectively. 36 14 Comparison of timedrawdown data from the transient model with a Theis (1935) type curve. A 1:500 block/fracture ratio of hydraulic conductivity is utilized. Match is at late time. Simulated values of T for blocks and fractures are 1.54E+006 ft2jday and 7.70E+008 ft2jday respectively, and for s, 0.1 and 0.01, respectively. 37 15 Comparison of timedrawdown data from the transient model with a HantushJacob (1956) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at vii
PAGE 9
early time. simulated values of T for blocks and fractures are 1.54E+006 ft2fday and 1.54E+007 ft2fday respectively,and for s, 0.1 and 0.01, respectively. 39 16 Comparison of timedrawdown data from the transient model with a Hantush (1960) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at early time. Simulated values of T for blocks and fractures are 1.54E+006 ft2fday and 1.54E+007 ft2fday respectively, and for S, 0.1 and 0.01, respectively. 40 17 Comparison of timedrawdown data from the transient model with a Neuman (1975) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at late time. Simulated values of T for blocks and fractures are 1.54E+006 ft2fday and 1.54E+007 ft2fday respectively, and for S, 0.1 and 0.01, respectively. 18 Comparison of timedrawdown data from the transient model with a Moench (1984) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at late time. Simulated values of T for blocks and fractures are 1.54E+006 ft2fday and 1.54E+007 ft2fday respectively, and for s, 0.1 and 0.01, respectively. 42 19 Comparison of timedrawdown data collected from the Floridan aquifer (Richard c. Fountain and Assoc., Jan., 1975) with a Theis (1935) type curve. Match is at early time. 43 20 Comparison of timedrawdown data collected from the Floridan aquifer (Richard c. Fountain and Assoc., Jan., 1975) with a Theis (1935) type curve. Match is at late time. 44 21 Comparison of timedrawdown data collected from the Floridan aquifer (Geraghty and Miller, Inc., Sept., 1976) with a Theis (1935) type curve. 46 viii
PAGE 10
NUMERICAL AND ANALYTICAL MODELING OF DOUBLEPOROSITY AQUIFERS by John A. Powers An AbStract Of a thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in the Department of Geology in the University of South Florida April, 1992 Major Professor: Mark T. Stewart, Ph.D ix
PAGE 11
Groundwater flow in a doubleporosity aquifer can be simulated using MODFLOW, a finitedifference, numerical code (McDonald and Harbough, 1985). The simulated aquifer is unconfined, and a doubleporosity system is represented by discrete blockandfracture units. The simulated fractures are vertical, orthogonal, and interconnected. The steadystate results reveal that the ratio of block to fracture hydraulic conductivity controls the hydraulic behavior of the regional flow system. As the ratio decreases, the simulated effects of the fractures become more apparent. When the hydraulic conductivity of the fracture exceeds the block hydraulic conductivity by several orders of magnitude, local flow systems develop within each block. The blockandfracture model was stressed to simulate response to a pumping well. Simulated timedrawdown data were analyzed using standard analytical solutions including Theis (1935), and compared with pumping test data from the Floridan aquifer. Two Theis curves can be fitted to timedrawdown data from the simulated aquifer; one at early time and one at late time. At early time, the effect of fractures on the simulated pumping test is large because the transmissivity value obtained from the earlytime match greatly exceeds the transmissivity value used in the X
PAGE 12
simulation for the blocks. The value of transmissivity obtained from the latetime match, however, closely approximates the block transmissivity used in the simulation. The transmissivity values obtained from the earlyand latetime matches increase as the simulated fracture transmissivity is increased two orders of magnitude. The increase in the earlytime transmissivity is greater than the increase of transmissivity for the latetime match. Deviation from the Theis type curves increases as the blocktofracture ratio o f hydrauli c conductivity is decreased because the model is increasingly nonhomogeneous. The simulated response explains some aspects of some pumping tests performed in the Floridan aquifer. Other timedrawdown data taken from pumping tests of the Floridan aquifer, however, match one Theis curve and are interpreted as representing a homogeneous aquifer. Abstract approved: Major Mark T.' stewart, Ph. D Professor, Department of Geology Date of Apprdval xi
PAGE 13
1 INTRODUCTION General The purpose of this study is to numerically simulate a doubleporosity aquifer and determine the principal variables that control the hydraulic response of such an aquifer. Doubleporosity aquifers are heterogeneous and anisotropic systems consisting of two porosity regimes. Primary porosity exists within the block matrix of the aquifer and is considered original in that it formed during petrogenesis {Streltsova, 1976). Fractures are secondary features that superimpose hydraulic properties different from those of the matrix blocks on the overall system {Barenblatt, Zheltov and Kochina, 1960). The fractures are often enhanced by dissolution or weathering (Barenblatt, Zheltov and Kochina, 1960). Traditionally, fractures are conceptualized as zones of high permeability with little or no storage capacity, and blocks are characterized as having low permeability with high storage capacity (Pirson, 1953). Fractures are thought to serve as major conduits for groundwater flow, while blocks function as storage areas {Moench, 1984). The contribution of fractures to groundwater storage is considered to be small, and some studies have
PAGE 14
assumed fracture storage to be insignificant (Barenblatt, Zheltov and Kochina, 1960). 2 An extensive network of regional fractures has been identified in westcentral Florida (Vernon, 1951; Faulkner, 1973). The fractures are thought to serve as conduits of groundwater discharge to springs (Faulkner, 1973; Sinclair, 1978). Accordingly, regional aquifers in this area may be characterized by doubleporosity behavior. Analytical solutions exist for doubleporosity aquifers. Such solutions are not used in this study, because, in order to solve the governing equations simplifying assumptions are required. Very simple geometrical configurations for blockandfracture systems are often utilized (van GolfRacht, 1982). Other assumptions are: no groundwater flow in blocks, and flow only in fractures; the geometry of blockandfracture systems need not be specified (Warren and Root, 1963); storage contribution of fractures are negligible (Barenblatt, Zheltov and Kochina, 1960); groundwater flow to a production well is radial; and fractures are horizontal (for example, Boulton and Streltsova, 1978). Because doubleporosity aquifers are difficult to investigate analytically, the focus of this thesis is on numerical techniques. It must be recognized that numerical models have limiting assumptions just as do analytical solutions,
PAGE 15
although the assumptions are of a different type and not nearly as constraining (Mercer and Faust, 1980) 3 The approach used in this study is to create block regimes and fracture regimes in a simple, numerical model. The model consists of one layer. Variable discretization is used to divide the model grid into narrow fractures and larger blocks. Actual field data are not utilized because the goal of this study is to investigate general relationships and not a specific field condition. The simulated aquifer is assumed to be unconfined and isotropic. Doubleporosity relationships are incorporated into the model by spatially defining a network of blocks and fractures with separate hydraulic characteristics. Inflow to the simulated aquifer is through uniform recharge, discharge is by a regional spring or wells. Both steadystate and transient conditions are investigated. The steadystate case reveals how the simulated aquifer behaves regionally in response to uniform recharge and simulation of a regional freshwater spring. The transient case reveals more localized effects, such as the response of the aquifer to pumping wells. Timedrawdown data from the transient case are qualitatively compared to selected, timedrawdown data collected from pumping tests of the Floridan aquifer. Standard analytical solutions are applied to both model and
PAGE 16
4 field timedrawdown curves. Analytical solutions applied in this study include homogeneous, transient flow to a pumping well (Theis, 1935), recharge from vertical infiltration (Hantush and Jacob, 1955), recharge from release of storage from semiconfining layers (Hantush, 1960), delayed drainage (Neuman, 1975), and block/fracture flow (Moench, 1984). Objectives The major objectives for this thesis are: {1) for both the transient and steadystate cases, construct and compare simulations of an unconfined system under doubleporosity and homogeneous conditions, utilizing variablegrid discretization; (2) evaluate possible interpretations of timedrawdown data from the heterogeneous simulations by applying standard analytical solutions; and (3) qualitatively compare timedrawdown data from the numerical solutions to selected, timedrawdown data collected during pumping tests of the Floridan aquifer.
PAGE 17
5 PREVIOUS INVESTIGATIONS Field Studies Aerial photographs taken from high altitude reveal the presence of largescale, linear features throughout the United States (Parizek, 1976). These linear features are unaffected by local topography and are thought to represent the surface expression of zones of increased concentration of fractures (Parizek, 1976). Within the photolinears there are often narrow zones where porosities and permeabilities can be 10 to 1000 times that of the surrounding intergranular blocks (Parizek, 1976). In Florida, the work of Stewart and Wood {1990) and Moore and Stewart {1983) has shown that fractures and secondary dissolution features are a significant component of the Floridan aquifer. Geophysical field work in the Tertiary aquifer of Florida has shown that fracture traces can be as wide as 100140 meters with a highporosity zone 40 meters wide located within the larger fracture trace (Stewart and Wood, 1990). Fracture zones associated with fracture traces usually have high hydraulic conductivities and are ideal places to install water wells due to increased flow of ground water (Parizek, 1976). Bulk porosity decreases away from the inner high porosity zone, most likely due to cementation
PAGE 18
(Stewart and Wood, 1990). In Florida, faults and fractures have developed in two principal directions. The primary system strikes northwestsoutheast and a secondary system trends northeastsouthwest (Vernon, 1951). They intersect one another at nearly ninety degrees and are often found to coincide with photolinears (Vernon, 1951). Modeling Techniques: Fractured Media 6 Fractured aquifers can be simulating utilizing the singleporosity approach or the doubleporosity approach. The singleporosity method can involve (Ward and others, 1987): a) modeling a fractured aquifer by assuming a hydraulically equivalent porous medium, b) individually simulating each fracture to be a discrete hydraulic conduit or c) a combination of a and b. Mathematical studies (Snow, 1965) have been attempted to determine appropriate parameters needed to describe singleporosity media. Snow (1965) developed a permeability tensor to characterize the equivalent permeability of fractured rocks from fracture geometry alone. This approach assumes that the extreme heterogeneity of fractured systems can be approximated by a single tensor. This approach is useful, but only if the density of fractures is large enough. As the fracture density increases, the system can be equivalent to a homogeneous medium. According to Long and others (1982), there may not exist a permeability tensor that can correctly characterize a fractured aquifer at any scale. Reasons for
PAGE 19
7 this are that the representative elementary volume required to describe a fractured aquifer may actually be larger than the total volume of the aquifer itself, or the fractures may be hydraulically isolated from one another so that, at any scale, the singleporosity approach is invalid. Individual fractures have also been analyzed using the singleporosity method (Witherspoon and others, 1981) under the assumption that the contribution of individual fractures in a system can be summed to arrive at the total permeability of the fractured rock mass. Network analysis (Ward and others, 1987) is used to separately model each individual fracture, but the amount and cost of the data required for this approach can be formidable. Flow in isolated fractures is investigated (Maini, 1971) by varying the conditions governing the flow regime of the fracture; however, any contribution to the flow regime by the nonfractured portion is ignored. Development of DoublePorosity Concept Muskat (1937) stated that fractures in certain carbonate aquifers serve as conduits of groundwater flow, with the block matrix serving as the primary source of groundwater storage. Fractures in doubleporosity aquifers contribute little to primary porosity but increase the permeability of the overall system by several orders of magnitude (Streltsova, 1976). Fractures provide the flow
PAGE 20
8 channels by which ground water moves through the system, but the storage of ground water is primarily due to the porosity of the block medium. Barenblatt, Zheltov and Kochina, (1960) were among the first to propose an analytical solution to the hydraulic aspects of groundwater flow in doubleporosity media. Transfer of fluid from the block matrix to fractures is dependent upon and proportional to head differences between the primary porosity of the blocks and secondary porosity of the fractures (Pollard, 1959). In addition, Barenblatt, Zheltov and Kochina (1960) stated that viscosity of the fluid and the geometry of blocks and fractures also play an important role in the transfer of the fluid. Barenblatt, Zheltov and Kochina (1960) assumed that the storage contribution of fractures to the overall system is so small that it can be essentially ignored, an assumption that may not be theoretically sound. Analytical Models Since Barenblatt, Zheltov and Kochina (1960), many analytical approaches have been utilized that assumed various geometrical configurations as solutions to fluid flow in doubleporosity media. Warren and Root (1963) assumed pseudosteadystate flow conditions while Kazemi (1969) assumed transient flow conditions. Pseudosteadystate flow is a mathematical approximation of transientflow. The advantage of pseudosteadystate flow conditions over transient flow conditions is that the mathematical
PAGE 21
complexity of the analytical solution is considerably reduced because the hydraulic head in the blocks is treated as spatially constant and therefore, divergence of flow within the blocks can be ignored. The geometry of the blocks need not be specified for pseudosteadystate flow. Warren and Root's (1963) model assumed that fluid flux occurs in response to differences in the average hydraulic head of the blocks and the fractures. Moench (1984) examined the model of Warren and Root (1963) and concluded that assumptions inherent in pseudosteadystate flow conditions are not theoretically sound unless there is a fracture skin of very low permeability surrounding the blocks. The fracture skin impedes the free exchange of fluid from block to fracture. Fracture skins can exist in carbonate aquifers and are essentially thin mineral accretions that are deposited at the block/fracture interface (Moench, 1984). The hydraulic conductivity of a fracture skin is usually much lower than that of the block it surrounds (Moench, 1984). According to Moench (1984), when a fracture skin is present most of the change in hydraulic head will occur across the fracture skin; as a result, the head gradient in the block will be small, justifying pseudosteadystate conditions. A feature in many analytical solutions of radial flow in doubleporosity media is the assumption that fracture storage is negligible and can be ignored (Barenblatt, 9
PAGE 22
10 Zheltov and Kochina, 1960). The omission of fracture storage, however, has no sound theoretical basis (Moench, 1984). Pumping tests reveal that fractures contribute the bulk of fluid flow to the well in early time. At late time, fluid flow from the blocks dominates (Moench, 1984). Fracture storage, although small, is significant. Boulton and Streltsova (1978) developed new equations for drawdown for doubleporosity media. Transient flow to a pumping well in a watertable aquifer is assumed. The model consists of two layers that have different hydraulic properties and are referred to as the block and fissure. Type curves for the drawdown in fissure and block were obtained by theoretically varying the hydraulic properties of the two layers. Flow in a confined aquifer was also analytically investigated by Boulton and Streltsova (1977a,b) using the same block and fissure geometry used for unconfined conditions. Type curves are presented assuming either the well is uncased in both the block and fissure (Boulton and Streltsova 1977a), or the well is cased only in the fracture (Boulton and Streltsova, 1977b). Flow is radial in both analytical models with vertical flow occurring in both the block and fissure for Boulton and Streltsova, (1977a) while only occurring in the block for Boulton and streltsova, (1977b). Both models assume block and fissure are compressible.
PAGE 23
11 Wang and Narasimhan (1985) constructed a statistical, numerical model to investigate what happens when the fractures of an unconfined system are partially dewatered. The effective hydraulic conductivity of a fracture drops dramatically as fluid saturation decreases and the amount of air in the fracture increases. Capillary forces hold water at the blockfracture interface, reducing fluid flow between matrix blocks. Only where the fluid phases of blocks intersect can water negotiate a fracture and move between blocks. Numerical Models Jones (1985) describes a numerical model of the Northern Withlacoochee Basin of the Southwest Florida Water Management District used to estimate the effects of future regional groundwater development. The model results were used to estimate the effect of withdrawals on Rainbow Springs and Silver Springs. The model was calibrated to steadystate, predevelopment conditions. Once calibrated, the model was stressed by pumping wells to assess the effects of groundwater development on the Surficial Aquifer, Floridan aquifer and major springs in the basin. The model has two layers, the upper one representing the Surficial Aquifer and the lower one representing the Floridan aquifer. Leakance between the upper and lower layers was determined by a headdependant flux term in the
PAGE 24
12 finitedifference equations. Flow in the confining layer (Intermediate aquifer) separating the Surficial aquifer from the Floridan aquifer was not directly simulated. Both aquifers were modeled assuming homogeneous and isotropic conditions. Initial estimates of transmissivity for the lower layer were too low as most of the data were extrapolated from outside the basin. In order to calibrate the model, large transmissivity values characteristic of fractures were necessary to simulate known spring discharges and the observed head values. By calibrating to averageannual steadystate conditions, the contribution of block flow to the system was inadvertently ignored. When the steadystate model was used as the starting point for verification of a transient simulation, the historical behavior of heads in the basin was not accurately predicted. The historical volume of spring flow was also not accurately predicted. During periods of low precipitation, simulated spring discharges approached zero. The natural springs exhibit nearly constant discharges throughout the year. Guvanasen and others (1987) also modeled the Floridan aquifer in the Northern Withlacoochee basin for the Southwest Florida Water Management District. The model boundaries were modified from that of Jones (1985) in order to remove boundary problems resulting from sparse data.
PAGE 25
13 Instead of uniformly varying transmissivity throughout the basin, statistical estimation was used. In this manner, spatial variation of transmissivity and other hydrologic parameters in the calibrated model was estimated. Local variation of transmissivity is much larger in the Guvanasen and others (1987) model than the Jones (1985) model. Guvanasen and others (1987) first simulated the Northern Withlacoochee basin as equivalent porous media, as was done in the Jones (1985) model. Unusually large values of hydraulic conductivity were needed to calibrate the model. A discrete fracture network was then identified from the work of Vernon (1951) and simulated in the model by the use of line elements. The simulated fractures substituted for the high values of hydraulic conductivity used in the equivalent porous model. The results of the simulation indicate that the orientation and connectivity of the fracture network has a larger effect on the behavior of the simulated springs then the hydraulic conductivity and storativity. Also, not all fractures must be simulated in order to achieve a calibrated model. Fractures with orientations parallel to regional flow directions have the greatest hydraulic effects, while fractures oriented perpendicular to flow have little effect on flow. The fractures also were never revealed at a regional scale by deflections of contours of hydraulic equipotentials.
PAGE 26
14 METHODS Description of Numerical Code The modular, threedimensional groundwater flow model (MODFLOW) developed by McDonald and Harbaugh (1984) for the United States Geological Survey was employed in this study to simulate a doubleporosity aquifer. MODFLOW is a blockcentered, finitedifference computer code that uses the strongly implicit procedure or the slicesuccessive overrelaxation procedure to solve the cellular equations of groundwater flow. The code is designed to be modular in structure so it can be altered and recompiled "in parts" to meet the demands of the user without changing the entire code. Most of the packages of the code, however, were not needed for this project. A doubleporosity aquifer was simulated utilizing the basic, blockcentered flow, well, recharge and stronglyimplicit procedure packages. Transient and steadystate conditions can be specified, and confined or unconfined conditions can be simulated. Development of Numerical Model Both steadystate and transient conditions are investigated. The steadystate model investigates the regional response of the aquifer to uniform recharge and a
PAGE 27
15 regional spring; the transient model emphasizes localized effects, such as the aquifer response to pumping wells. Doubleporosity characteristics are incorporated into the steadystate and transient models through spatial variation of hydraulic conductivity. The block matrix is simulated by much lower hydraulic conductivities, compared to the fracture matrix (block/fracture ratios of hydraulic conductivity: 1:10, 1:100, 1:500). Approximate logarithmic progression in hydraulic conductivity was implemented in order to determine at what level of hydraulic conductivity contrast the system no longer can be interpreted to be a homogeneous aquifer. In addition, blocks are defined as having an orderofmagnitude larger storage capacity than fractures. Hydraulic conductivity in the model varies in the horizontal plane but is constant in the vertical plane. All fractures in the system are vertical and intersect at ninety degrees, which approximates the orientation of fracture traces found in the Floridan aquifer (Faulkner, 1973). Because MODFLOW is a finitedifference code, all cells are rectangular and ground water must flow perpendicular to cell faces; as a result, fractures are parallel to the principal axes of the model. Anisotropic effects are not investigated.
PAGE 28
16 Fracture zones have been defined as areas of increased fracture density with each zone consisting of many individual fractures that collectively comprise a fracture trace (Parizek, 1976). Accordingly, it is difficult to model a fracture zone as it exists in the field due to the scarcity of data defining individual fractures (Guvanasen and others 1987) and memory limitations of computers. Simplifying assumptions, therefore, must be made. Fracture zones are simulated in this study as uniform zones of higher hydraulic conductivity. These simulated fracture zones have widths of tens of meters, as suggested by Stewart and Wood (1990). Simulations of heterogeneous and homogeneous aquifers are compared. The homogeneous simulation acts as a control or basis of comparison for the heterogeneous simulations. For comparison purposes, both have the same grid dimensions and discretization scheme. Geraghty and Miller's Aquifer Test Solver (AQTESOLV; Duffield & Rumbaugh, 1989) is used to interpret the simulated pumpingtest data, as well as field timedrawdown data. AQTESOLV enables the utilization of a variety of analytical solutions for many types of aquifer tests. Aquifer parameters can be automatically estimated using the Marquardt nonlinear leastsquares technique for best match points. In addition, AQTESOLV permits the interactive
PAGE 29
17 matching of type curves to pumpingtest data on the computer screen. Discharge from unconfined aquifers can occur in a variety of ways: spring flow, pumping wells, downward leakance, evapotranspiration and lateral discharge. Downward leakance is not simulated, because the simulated aquifer is assumed to rest on an impermeable base. Lateral discharge is not required because the simulated basin is internally drained. Evapotranspiration is not simulated because model recharge is assumed to be effective recharge. Spring flow is the only component of discharge for the steadystate model used here. In westcentral Florida, groundwater discharge occurs primarily through freshwater springs (Faulkner, 1970). There are two firstorder springs located in the Northern Withlacooche Basin on the west coast of Florida: Rainbow Springs and Silver Springs (Meinzer, 1923). According to Faulkner (1970), 92 percent of the groundwater discharge from the Northern Withlacooche Basin is through these springs. Silver Springs discharges an average of 522 million gallons per day while Rainbow Springs discharges an average of 461 million gallons per day (Jones, 1985). One spring was simulated in this study by specifying a node with a constant head of zero, forcing all flow within the model to discharge through the spring cell.
PAGE 30
18 A single discharge well is used to investigate pumping effects in the transient model. MODFLOW simulates a pumping well by placing an extra flux term in the flow equation representing the cell being stressed. It is important to note that the flux due to the "well" is applied or averaged throughout the entire cell and not at a point. The well is simulated as fully penetrating the aquifer. Variablegrid discretization is used to divide the steadystate model into a series of intergranular blocks separated by singlecell fractures. Each block is composed of many cells. The steadystate model consists of 47 columns and 47 rows. Figure 1 illustrates the grid discretization. Fractures are highlighted. It was the original intent of this study to use the same discretization scheme and basin size for the transient and steadystate models. Because of barrierboundary effects in the transient model, however, grid dimensions had to be expanded. This distance proved to be so large that not enough computer memory could be allocated to keep the same discretization scheme used in the steadystate model. The discretization scheme remained the same within the core of the transient model; larger nodes were used towards the model boundaries. A fracture width of 50 feet was used for both the transient and steadystate models. The transient model grid consists of 67 columns and 67 rows.
PAGE 31
........__,,, Block Fracture Spring Node Figure 1. Boundary conditions and grid discretization used in the model. The model is unconfined and consists of one layer with noflow boundaries on all sides and bottom. The model is 5,350 feet on each side and 100 feet thick. A vertical exaggeration of approximately lOX is used here. The arrow indicates direction of groundwater flow. The spring node is defined with a constant head of zero. r\0
PAGE 32
20 Boundary Conditions In this study the modeled system is that of an internally drained basin with noflow boundaries on all sides. Flux out of the basin is accomplished either by a node representing a freshwater spring (steadystate case) or discharging well(s) (transient case). The spring node is defined with a constant head of zero. All regional flow out of the system must eventually exit this cell. Boundary conditions utilized in this study are depicted in Figure 1. The bottom of the simulated basin is a noflow boundary. Model Parameters Hydraulic conductivity, aquifer recharge, initial elevations of the watertable surface and aquifer bottom are specified in the steadystate model. Flux into the steadystate model is by aquifer recharge while flux out is through a constanthead node with zero head (spring node). Starting heads for all cells are zero. All model parameters are included in Table 1. Specific yield is used in the transient model. Recharge is not included because the objective is to determine the aquifer response to a pumping well in terms of drawdown. Groundwater flux in the transient model is due to release from storage; flux out of the model is through a discharge node representing a pumping well. A simulated
PAGE 33
Table 1. Model parameters. SteadyState Model Block Hydraulic Conductivity Fracture Hydraulic Conductivity Block/Fracture Ratio Recharge Basin Thickness 1:1 1:10 1:100 1:500 Transient Model Block Hydraulic Conductivity Fracture Hydraulic Conductivity Block/Fracture Ratio Block Specific Yield Fracture Specific Yield Discharge Well Basin Thickness 1:1 1:10 1:100 1:500 15.4 ft/d 15.4 ft/d 154 ft/d 1540 ft/d 7700 ft/d 0.025 ft/d 100 ft 15400 ft/d 15400 ft/d 154000 ft/d 1540000 ft/d 7700000 ft/d 0 1 0.01 4.0E+07 cubic ft/d 100 ft 21
PAGE 34
22 spring was not included in the transient model because it is a steadystate regional feature. The hydraulic conductivity for blocks and fractures in the transient model are three orders of magnitude larger than values used in the steadystate model in order to achieve greater drawdown. Note that while the absolute values of transmissivity are much larger in the transient model, ratios of block/fracture hydraulic conductivity remain the same. The fracture transmissivity for the block/fracture ratio of 1:10 falls within the range of published values for westcentral Florida, and the fracture transmissivity for ratios of 1:100 does not greatly exceed the published values (Jones, 1985; Geotrans, Inc., 1988). The fracture transmissivity for the 1:500 simulation exceeds the published values by an order of magnitude (Jones, 1985). A large increase in transmissivity will shift a loglog plot of drawdown data along the y axis but will not affect the shape of the plot.
PAGE 35
23 RESULTS SteadyState Model The pattern of equipotentials for a simulation using a 1:10 block/fracture ratio of hydraulic conductivity (Figure 2) is similar to a simulation of a homogeneous aquifer (Figure 3). However, heads in the fractured aquifer are about 40% lower than in the homogeneous case. As the ratio of block/fracture h ydraulic conductivity is lowered two orders of magnitude, nonhomogeneous effects first appear farthest upgradient from the outflow node. Equipotentials refract sharply across the fractures (Figure 4). At very low ratios of block/fracture hydraulic conductivity (1:500), local flow systems develop within the blocks (Figure 5). Fluidvelocity vectors of the homogeneousaquifer simulation are uniform in direction and magnitude towards the spring (Figure 6). When fracture hydraulic conductivity is increased one order of magnitude (Figure 7), velocities increase in the fractures. Vectors are no longer uniform in magnitude with respect to location in the flow field. When fracture hydraulic conductivity is increased two orders of magnitude, local flow (Figure 8) is divergent toward fractures. At very low ratios, localized flow (Figure 9) in
PAGE 36
24 Figure 2. Equipotentials (feet) from the steadystate model utilizing a 1 : 10 block/fracture ratio of hydraulic conductivity.
PAGE 37
.. """ L55ss so Figure 3. Equipotentials (feet) from the homogeneous, steadystate model. 25 lOll lOll
PAGE 38
26 .... .. .. I I I I I / _/ ""TO / L/ r,o/ 9 / ____/ v / / ..... \ !((d r \ \ i Figure 4. Equipotentials (feet) from the steadystate model utilizing a 1:100 block/fracture ratio of hydraulic conductivity.
PAGE 39
27 Figure 5. Equipotentials (feet) from the steadystate model utilizing a 1:500 block/fracture ratio of hydraulic conductivity.
PAGE 40
T apom 1snoaua6omoq UIO.:lJ pue g a.Inf>""F.!l .. ,.. I .. .... _. __ \ .._ .... ........__ ..... .. ... 1 I .. .. _. ........... .. ... ...... _,. ... .0. ... ........................... .. ... I .. ... """"'"" ... .. .. .. ... I ... ................................... ... ' I ... .,.,.,. .,. "' "' "' ... ,, "' I ..... .... ................ ... ... I ,.,., "' "' "' I I I 1 1 l ..... .... ",, ... ' ' I ''", "' .. I' I' I' I I 1 .... .... ''''' ' I ,.,.,. I' I' ,. "'' I 1 .... ',, ... ' ' ttOt J 1 I ? ,. /' I' I I I I I I I ' ''''' ' dc 1 1,." I I I I I I I t ''''' ' ''''' '\ '\ S'_ o">.1 It T ' "'' '\ 'Z', Ill T T T "'' ' "''' '\ ' '" '" "' "; m>ZL . 1 1 1 1 1 , , , 1 , I Ill T T t t t n t t ' ' "' , , , , , , , I I I Ill I I I I T t \\ t ' ' '" , 1 , 1 f tr I I I I Ill I I I I T t tT t ' ' n ,,, , , , , I I I 1 I t I t t t I t !f t I I 1 t n 1 I , , , , , , I t I Ill I I t I t !f t t t t I !f t t t 1 , 1 , I II I 1 I I I I Ill I I I I t !! t t t t t t t t 1 , , , , , , I I I I fl I I t I t !! I t t T ' 1T I I I I I I , , , I , 1 1 II 1 I I I 1 t 1T t I t 1 1T t 1 I ff I I 1 1 , 1 1 I 1 f 1 tf I 1 erorft t t t tt I t t t I t Itt t I f I If I I f I 1 f Iff f f , I If I f I I I I ttl I I t I t ttl !0!0 f f f If t I I 1 f ff 1 I I I II I I I t t t tIt I I I t n t I 1 f f I f , ff 1 f I I f t1Si f I I t t t!! I I t t t It t I I 1 f f f f ff f f I I I II I t t I I t If t t t t t I If t I If 1 I f 1 If, f f f f Ill f f f I I I!! I t t I t I It If I I I f f 1 If I f f I I I I l l f f I t t If t t t t t t t" I I I I I , I If I I I I I I I I I I I t 1 1 1 If I t t t 1 If t I I I I I I I I I ft I I I I I I I f f t I f t If I I I t t If' .... If I f I I f I f f I f f I f I If f I I I I If t I I I t I t I t Iff I I f I I If I f f I f I It I 1 I f I 1 If I t I I t I I II I f f f I f f I If I I I I f I If I I I I I t I!! I I I t I I If Ill I I f I I I I I I t I t I If t I I I t t t tf t t t t I t I I I I If I f fl I f f I tf I I I I I I I ft I I I t I I I I f t I I I Iff I f , I !! I I f f I I !! I t t I t I If t f Iff I I I I I If I I I I I Iff I I I I t I II I t t t I I I t I II I I I I I fl I I I f I If I I f t I t If I I I I t I tIt f f I , 1 t ft I I f I ttl I I I I I I If I I I I I Itt .... . f t f t , I I II I t I I I fl I f I t I t !It I t t 1 I t t 1 'I , , , , , : I : I f II f t ) t t t t t t , . .... ' "' ' I 'I . ' '; t 9 I l .... .... ... ..... .. sz
PAGE 41
uMoqs JO ot:t tapom L .. l ll+C\ .. C ... j ... OJ ... ...... I mu ..... ..... ..... .... ..... ..... ..... .... .. 6Z
PAGE 42
JO Oot:t tapom a I
PAGE 43
31 Figure 9. Fluidvelocity vectors and equipotentials (feet) from the steadystate model utilizing a 1:500 block/fracture ratio of hydraulic conductivity. Fractures are shown.
PAGE 44
some instances is opposite in direction to regional flow, and block flow is radial toward fractures. Transient Model 32 Timedrawdown data from the simulation of a homogeneous aquifer show an excellent fit with a type curve of the Theis (1935) analytical solution (Figure 10). Values of transmissivity and specific yield correspond closely to actual simulated values. The transient model approaches steadystate conditions after ten days of pumping. Timedrawdown data from the simulation of a heterogeneous aquifer utilizing a block/fracture ratio of 1:10 deviate from the Theis (1935) type curve. Theis type curves can be matched to earlytime (Figure 11) and latetime (Figure 12) portions of the timedrawdown plots. For each plot, the earlytime value for transmissivity exceeds the transmissivity value for late time, and specific yield for the l atetime match exceeds that for the earlytime match. Theis type curves are also matched to early(Figure 13) and latetime (Figure 14) portions of timedrawdown plots in which the block/fracture ratio of hydraulic conductivity is decreased to 1:500. Again, the earlytime value for transmissivity exceeds the value for late time, and specific yield for the latetime match exceeds that for the earlytime match. The difference between earlyand latetime values of hydraulic conductivity and specific
PAGE 45
r M l 0 'tl l rti Q 'tl 0,) +.l u 0,) ,.. 0 u 100. 10. 1 .. 1 .. 002 33 = 0.1.1.45 0.1 1. 10. Tin1e (dags) Figure 10. Comparison of timedrawdown data from the homogeneous, transient model with a Theis (1935) type curve. Simulated values of T and s are 1.54E+006 ft2fday and 0.1, respectively.
PAGE 46
100. 10. 1. 0.1 1.E002 34 = 8.89245 0.1 1. 10. Tin1.e (days) Figure 11. Comparison of timedrawdown data from the transient model with a Theis (1935) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at early time. Simulated values of T for blocks and fractures are 1.54E+006 ft2jday and 1.54E+007 ft2jday respectively, and for s, 0.1 and 0.01, respectively.
PAGE 47
s:: l 10. l lti Sot Cl 't:l OJ u OJ Sot Sot 0 u 1. 0.1 1.E002 35 0.1 1. 10. Ti:me Cdags) Figure 12. Comparison of timedrawdown data from the transient model with a Theis (1935) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at late time. Sifulated values of T for blocks and fractures are 1. 54E+006 ft /day and 1. 54E+007 ft2 /day respectively, and for s, 0.1 and 0.01, respectively.
PAGE 48
36 roo 10. +l = 'H 0 oo ....., = 00 00 00 "" 0 l 00 0 0 0 '0 0 0 l 0 0 ltj 0 1. 0 0 "" o/ Q '0 a/" Q) +l / u Q) /' .. 1 0 u 0.1 1.E002 0.1 1. 10. Tin:te Cdags) Figure 13. Comparison of timedrawdown data from the transient model with a Theis (1935) type curve. A 1:500 block/fracture ratio of hydraulic conductivity is utilized. Match is at early time. Simulated values of T for blocks and fractures are 1.54E+006 ft2fday and 7.70E+008 ft2fday respectively, and for s, 0.1 and 0.01, respectively.
PAGE 49
t: l 0 'tl l lti 0 'tl OJ +l 0 OJ 0 u 10. 1. 0.1 1.E002 0.1 1. 10. Tin1e Cdags) 37 Figure 14. Comparison of timedrawdown data from the transient model with a Theis (1935) type curve. A 1:500 block/fracture ratio of hydraulic conductivity is utilized. Match is at late time. Sifulated values of T for blocks and fractures are 1.54E+006 ft fday and 7.70E+008 ft2/day respectively, and for s, 0.1 and 0.01, respectively.
PAGE 50
yield increases as the block/fracture ratio of hydraulic conductivity decreases. 38 Timedrawdown data from a simulation of a heterogeneous aquifer, utilizing a block/fracture ratio of 1:10, are investigated using two, leakyaquifer analytical solutions (Hantush and Jacob, 1955; Hantush, 1960). Type curves from the Hantush and Jacob (1955) and Hantush (1960) analytical solutions can be matched to the simulated, timedrawdown data at early time (Figures 15 and 16 respectively). Neuman's {1975) analytical equations are also applied to the simulated, timedrawdown data utilizing a block/fracture ratio of 1:10. A plot of the Neuman type curve matches the timedrawdown data from the simulation of a heterogeneous aquifer {Figure 17) in a manner similar to the type curve from the Theis analytical solution. The timedrawdown data are also compared to type curves from Moench's {1984) analytical solution (Figure 18) for doubleporosity aquifers. The match between the Moench type curve and the timedrawdown data from the simulation of a heterogeneous aquifer appears to be similar to the match between the Theis type curve and timedrawdown data. Timedrawdown data from the Floridan aquifer in Manatee County, are shown in Figures 19 and 20 {Richard c. Fountain and Assoc., Jan., 1975). Separate Theis type curves match the pumpingtest data at early and late time. Other,
PAGE 51
100. 10. 1. = 2.6275E+006 = 0 .09385 1 .E005 0.1 1.E002 0.1 1. 10. Tirne (dags) 39 Figure 15. Comparison of timedrawdown data from the transient model with a HantushJacob (1955) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at early time. Simulated values of T for blocks and fractures are 1.54E+006 ft2/day and 1.54E+007 ft2/day respectively, and for s, 0.1 and 0.01, respectively.
PAGE 52
= = = i.E005 10. 1. / 0.1 1.E002 0.1 ooo o o o o_ 1. Tirne (dags) 40 o oo o o oo o 10. Figure 16. Comparison of timedrawdown data from the transient model with a Hantush (1960) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at early time. Simulated values of T for blocks and fractures are 1.54E+006 ft2jday and 1.54E+007 ft2jday respectively, and for S, 0.1 and 0.01, respectively.
PAGE 53
100. 10. 1. 0.1 1.E002 = 0.02076 = 0.00.1372 0.1 1. Tin1e (dags) 41 10. Figure 17. Comparison of timedrawdown data from the transient model with a Neuman (1975) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at late time. Simulated values of T for blocks and fractures are 1.54E+006 ft2/day and 1.54E+007 ft2/day respectively, and for s, 0.1 and 0.01, respectively.
PAGE 54
100. 10. 1. 0.1 1.E002 = 1.546E+664 = 6.661193 = 1.8378E+666 1.8071E665 = 8. 6. 0.1 1. Tin1.e (da gs) 42 10. Figure 18. Comparison of timedrawdown data from the transient model with a Moench {1984) type curve. A 1:10 block/fracture ratio of hydraulic conductivity is utilized. Match is at late time. Simulated values of T for blocks and fractures are 1.54E+006 ft2fday and 1.54E+007 ft2/day respectively, and for s, 0.1 and 0.01, respectively.
PAGE 55
r"' ..., '14 ..., t: l 0 'C l It! Q 10. 1 0.1 1.E002 1.E004 = 43 0 0 0 00 0 0 0 0;:;>, 9' / J l I 9 1.E002 1. Ti:me Cda'gs) Figure 19. Comparison of timedrawdown data collected from the Floridan aquifer (Richard c. Fountain and Assoc., 1975) with a Theis (1935) type curve. Match is at early time.
PAGE 56
44 = 4.5922E004 f"'' .., 1 .. """ 9:'o7 l 0 0 'tl 0 l rti r .... 0.1 0 0 1 .. E002 1 .. E004 1 .. E002 1 .. Ti1ne Cdags) Figure 20. Comparison of timedrawdown data collected from the Floridan aquifer {Richard c. Fountain and Assoc., 1975) with a Theis {1935) type curve. Match is at late time.
PAGE 57
45 pumpingtest data taken from the same area of Florida (Geraghty and Miller, Sept. 1976), however, are matched by a single Theis type curve from early to late time (Figure 21).
PAGE 58
46 1. 0.1 10. 100. 1000. 10000 Tin'le (min) Figure 21. Comparison of timedrawdown data collected from the Floridan aquifer (Geraghty and Miller, Inc., 1976) with a Theis (1935) type curve.
PAGE 59
47 DISCUSSION Steadystate Model The steadystate model was found to be sensitive to the ratio of block/fracture hydraulic conductivity. When fracture hydraulic conductivity is one order of magnitude larger than block hydraulic conductivity (Figure 2), the heterogeneity in the regional flow system is not apparent. The pattern of the equipotentials is similar to that of the simulation of a homogeneous aquifer (Figure 3). They are similar in that both sets of equipotentials are essentially concentric with respect to the outflow node. Heads are lower in the heterogeneous case, however, reflecting the greater transmissivity of the system as a result of the presence of fractures. At this level of fracture density, the contrast in hydraulic conductivity between blocks and fractures is not large enough to reveal the presence of fractures in the regional flow system. Without additional knowledge of the system, examination of contoured heads would indicate the system is homogeneous. As fracture hydraulic conductivity is raised two orders of magnitude compared to the blocks, the regional flow system (Figure 4) can no longer be interpreted to be a
PAGE 60
homogeneous aquifer. Nonhomogeneous effects first appear farthest upgradient from the outflow node. Contours refract across the fractures. The contrast in hydraulic conductivity is now large enough to reveal the presence of fractures away from the spring where the groundwater gradient is small. 48 When block hydraulic conductivity is much smaller than fracture hydraulic conductivity (block/fracture ratio of 1:500), the heterogeneity of the aquifer is readily apparent as localized areas of divergent flow coincide with each block (Figure 5). At this level of heterogeneity, each individual block supports a local flow system within the regional aquifer. Block flow is outward toward the fractures. Notice that close to the outflow node the local flow systems are not apparent. The groundwater gradient is too high. Also, along the sides of the model, local flow cells appear to be bisected by the basin boundary and in addition, the direction of groundwater flow in the corners of the basin is outward. Both areas reveal the interaction of a bounded basin with the heterogeneity of the aquifer. Fluidvelocity vectors further define how the behavior of the system depends upon the ratio of block/fracture hydraulic conductivity. Every particle of water i n the system must e xit the model through the outflow node. For the simulation of a homogeneous aquifer (Figure 6), vectors
PAGE 61
49 are uniform in direction towards the spring. Also, velocities increase towards the outflow node. When the block/fracture ratio of hydraulic conductivity is decreased to 1:10, vector velocities increase in the fractures (Figure 7). Vectors are no longer uniform in magnitude and direction with respect to location in the flow field. As the ratio of block/fracture hydraulic conductivity is decreased further (1:100), local flow (Figure 8) is increasingly divergent toward the fractures. At a 1:500 block/fracture ratio of hydraulic conductivity, local flow directions (Figure 9) are, in some instances opposite in direction to regional flow. Directions of flow in blocks are outward toward the fractures. This effect is less pronounced near the outflow node. Also, the difference between maximum and minimum velocities increases linearly as the ratio of block/fracture hydraulic conductivity decreases. Transient Model Transient simulations were completed to assess the response of the system to a pumping well. Timedrawdown data from the simulation of a heterogeneous aquifer were plotted and interpreted by standard analytical solutions. The objective was to determine how the data would be interpreted assuming limited information was available about the aquifer as is almost always the case for field tests.
PAGE 62
50 The graphical analysis of pumpingtest data is difficult due to the nonuniqueness of interpretation. The timedrawdown response for many different types of aquifers is similar, for example: slow drainage of unconfined aquifers due to vertical anisotropy and confined, doubleporosity aquifers. The groundwater scientist must decide which analytical model to apply assuming that assumptions in the analytical solution are also true for the field data. Choice of analytical method is problematic unless the hydraulic characteristics of the aquifer are well known. The Theis analytical solution {Theis, 1935) is often the first choice when analyzing pumpingtest data because deviations from the Theis type curves are well understood and documented. It is the simplest analytical model to apply. If pumpingtest data deviate from the Theis type curves in a systematic way, interpretations can be made of the aquifer. Pumpingtest data from the simulation of a homogeneous aquifer are matched by a single Theis type curve from early to late time as shown in Figure 10. Values of transmissivity and specific y ield correspond closely to the modeled values. Pumpingtest data from the simulation of a heterogeneous aquifer (block/fracture ratio of 1:10) are not matched by the Theis type curve very well. Two Theis type curves can be fitted to the timedrawdown data, one at early time and one at late time, as shown in Figures 11 and 12,
PAGE 63
51 respectively. For the earlytime match, the simulated drawdown exceeds the values predicted by the Theis equation at late time. When a plot of a Theis type curve is matched to the latetime portion of the simulated drawdown, earlytime data exceed the values predicted by the Theis equation. The calculated transmissivity for the earlytime match exceeds the latetime match, and specific yield for late time exceeds that of early time. The difference between transmissivity values determined from earlyand latetime matches increases as the block/fracture ratio of hydraulic conductivity is decreased and fewer data fall on the Theis type curves. Deviation from the Theis curves is greater when the block/fracture ratio is decreased to 1:500 (Figure 13, earlytime match, Figure 14, latetime match). The earlytime response from the simulation of a heterogeneous aquifer can be attributed to the fractures while at late time, the flow system is dominated by the blocks. Type curves from other available analytical solutions do not successfully match the timedrawdown data from the simulation of a heterogeneous aquifer. The characteristics of the timedrawdown data are not duplicated by recharge from vertical infiltration (Hantush and Jacob, 1955), recharge from release from storage of semiconfining layers (Hantush, 1960), delayed drainage (Neuman, 1975), assumption of homogeneous conditions (Theis, 1935) or Moench's model of a doubleporosity aquifer (1984). A plot of the time
PAGE 64
52 drawdown data from the simulation of a heterogeneous aquifer does resemble the theoretical response of a finite, homogeneous aquifer to a pumping well (Freeze and Cherry, 1979). At late time the simulated drawdown are greater than values predicted by the Theis analytical equation, suggesting a change in transmissivity. The theoretical response of a semiconfined aquifer to vertical infiltration is a deviation from the Theis type curve at late time during a pumping test (Hantush and Jacob, 1955). Observed drawdowns are less than the values predicted by the Theis equation at late time. In contrast, the simulated, timedrawdown data exceed the values predicted by the Theis equation at late time. The application of the Hantush and Jacob (1955) analytical equations, as shown in Figure 15, to the simulated, timedrawdown data result in estimates of transmissivity and storativity similar to the match with the Theis type curve. The governing analytical equations reduce to the Theis equation because the simulation is nonleaky. Similar results are obtained when the timedrawdown data from a simulation of a heterogeneous aquifer are compared to Hantush's (1960) analytical solution. Hantush (1960) accounted for storage in the semiconfining layer in addition to vertical infiltration from overlying aquifers. Values of transmissivity and storativity resulting from an
PAGE 65
53 earlytime plot of a type curve (Figure 16) are similar to estimates from plots of the Theis and Hantush and Jacob type curves. The factor B is very small indicating again, that the Hantush (1960) equations reduce to the more general Theis equation. Neuman (1975) proposed an analytical solution to predict the response of the water table to pumping stresses in unconfined aquifers. The type curves have three segments. The first segment (early time) is similar to the response in confined aquifers where water is released instantaneously from elastic storage. Drawdown follows a Theis type curve during early time. During the second segment (transition) there is a decrease in the rate of drawdown relative to the earlytime portion of the type curve. Drawdowns from the third segment (late time) again conform to a Theis curve. Timedrawdown data from a simulation of a heterogeneous aquifer (1:10 block/fracture of hydraulic conductivity) are not duplicated by the Neuman (1975) analytical solution. During early time, drawdown follows a Theis curve (Figure 11). As time progresses there is an increase, rather than a decrease, in the rate of drawdown relative to the Theis curve. Values of specific yield and storativity derived from the Neuman solution (Figure 17) do not correspond well to values used in the simulation of a heterogeneous aquifer. Storativity (S) from the Neuman solution should correspond to fracture specific
PAGE 66
54 yield, but instead, is similar to the value for block specific yield used in the simulation of a heterogeneous aquifer. Also, specific yield (Sy) derived from the Neuman solution is similar to the value for fracture specific yield instead of block specific yield. In addition, the factor (B) is small enough that the governing analytical equations essentially reduce to the Theis solution. Because the numerical model in this study is an attempt to simulate a doubleporosity aquifer, the simulated, timedrawdown data are compared to an analytical solution for doubleporosity aquifers. The Moench analytical model (1984) was selected as it is a synthesis of the most widely used analytical techniques for doubleporosity aquifers. Type curves qualitatively appear to be similar to Neuman's family of type curves for unconfined aquifers, because there are three segments. Earlytime drawdown conforms to a Theis type curve and represents the response of fracture(s) to the pumping well. During early time, water is derived primarily from storage in the fractures. At late time, fluid is derived from storage in the block matrix. Between early and late time, the timedrawdown curve flattens due to flow from blocks to fractures as a result of a head differential between the two aquifer matrices and leakance from the blocks.
PAGE 67
55 Two schools of thought describe the transition portion of the timedrawdown curves. One approach assumes that flow occurs under fullytransient conditions (for example, Kazemi, 1969) and the other approach assumes pseudosteadystate conditions (Warren and Root, 1963). The mathematics of pseudosteadystate analytical solutions are simpler than that of transient analytical solutions. Moench (1984) has shown that the pseudosteadystate model is not valid because it ignores divergent flow in blocks. Field data exist that support both models. Moench (1984) maintains that the pseudosteadystate model is a special case within the transient model. Moench's analytical solution (1984) assumes transient flow conditions coupled with a fracture skin at the block/fracture interface. Transfer of fluid between the block and fracture is impeded at the fracture skin. The solution assumes the block matrix is composed of horizontal slabs interspaced with horizontal fractures. In contrast, fractures in the simulation of a heterogeneous aquifer are vertical. The application of Moench's curves to the timedrawdown data from the simulation of a heterogeneous aquifer results in a poor solution (Figure 18). Fracture hydraulic conductivity (K) from the plot of a Moench type curve closely corresponds to the block hydraulic conductivity used in the simulation. Fracture storage (Ss) derived from the Moench solution exceeds the block storage (Ss') from the
PAGE 68
56 Moench solution. In addition, hydraulic conductivity (K') for the block (Moench solution) exceeds the value for fracture hydraulic conductivity (K, Moench solution) by two ordersof magnitude. In addition, Theis curves do not match the timedrawdown data from the simulation of a heterogeneous aquifer in the same manner as described in Moench (1984). As time progresses there is an increase, rather than a decrease, in the rate of drawdown relative to the Theis curve. In other words, the earlytime Theis curve underpredicts the simulated drawdown but the Moench analytical solution requires the earlytime Theis curve to overpredict. The application of analytical equations that describe doubleporosity aquifers to the timedrawdown data from the simulation of a heterogeneous aquifer does not appear to be valid. Timedrawdown data from the simulation of a heterogeneous aquifer can be compared to field pumpingtest data collected from the Floridan aquifer. Timedrawdown data from the Floridan aquifer are matched to earlyand latetime, Theis type curves, as shown in Figures 19 and 20 (Richard c. Fountain and Assoc. Jan., 1975). In this area of Florida, the Floridan aquifer is very well confined (Ryder, 1985) resulting in minimal leakance from overlying aquifers. The field drawdown data deviate upward from the Theis type curve at late time in a similar manner as the timedrawdown data from the simulation of a heterogeneous
PAGE 69
57 aquifer. In contrast, other pumpingtest data taken from the same area of Florida, (Geraghty and Miller, Sept. 1976) fit the Theis type curve from early to late time (Figure 21). Field pumpingtest data that deviate from a Theis type curve are probably collected from observation wells situated in or near vertical fracture zones, thus explaining their upward deviation from a Theis type curve at late time. Depending upon location, pumpingtest data exist from the Floridan aquifer that are indicative of homogeneous and heterogeneous conditions. Therefore, the magnitude of local fracturezone density may be important in determining response to pumping. What is unclear, however, is how close to a fracture zone an observation well must be in order for timedrawdown data to deviate from a Theis type curve. The fracture density in the model is large enough that, regardless of position of the observation well, the resultant timedrawdown curves are similar. Results from the simulation of a heterogeneous aquifer indicate, however, that if a pumping well is located in a fracture, drawdown will be slightly greater for observation wells situated in fractures rather than blocks. In contrast, if a pumping well is located in a block, observed drawdown is solely dependant upon distance from the pumping well. Whether or not the observation well is located in a fracture or block is not a factor. If block size could be increased relative
PAGE 70
58 to fracture size, the position of the observation well might be important as timedrawdown curves taken from blocks would be different from data collected in fractures. The timedrawdown curve from an observation well located in a block would probably match a Theis type curve from early to late time, while data from observation wells nearer to fractures would deviate from the Theis type curve.
PAGE 71
59 CONCLUSIONS Under steadystate conditions, the simulations of a heterogeneous aquifer demonstrate that the ratio of block/fracture hydraulic conductivity is an important characteristic of doubleporosity aquifers. The recognition of a doubleporosity aquifer in the field is directly dependent upon the quality and quantity of field data. The density of available field data is rarely large enough to adequately describe the fracture matrix on a regional basis. The simulation of fractures in a regional model can be problematic, unless the focus is placed on local areas within the regional basin where field data describing fractures are more available. The timedrawdown data from the simulation of a heterogeneous aquifer do not fit the following analytical models well: recharge from vertical infiltration; recharge from release of storage from semiconfining layers; slow drainage due to vertical anisotropy; assumption of homogeneous conditions; or doubleporosity. The timedrawdown curves are similar, however, to the theoretical response of a finite aquifer to a pumping well.
PAGE 72
60 The application of the Theis analytical equation (1935) to timedrawdown data from the simulation of a heterogeneous aquifer yields two matches; one at early and one at late time. Drawdown at early time consists of fluid derived primarily from fracture storage while release from storage from blocks accounts for drawdown at late time. A greater deviation from the Theis type curves occurs at higher ratios because the model is increasingly nonhomogeneous. The transient model explains certain aspects of some pumping tests performed in the Floridan aquifer. Other timedrawdown data, however, taken from pumping tests of the Floridan aquifer match one Theis type curve and are interpreted as representing a homogeneous aquifer. Field data that do not conform to Theis type curves are probably collected from observation wells situated in or near fracture zones in the Floridan aquifer, thus explaining their deviation.
PAGE 73
61 LIST OF REFERENCES Barenblatt, G.I., I.P. Zheltov, and N. Kochina, 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. Journal of Applied Mathematics and Mechanics., 24 (5), pp. 12861303. Boulton, N.S. and T.D. Streltsova, 1977a. Unsteady flow to a pumped well in a twolayered waterbearing formation. Journal of Hydrology, 35, pp. 245256. Boulton, N.S. and T.D. Streltsova, 1977b. Unsteady flow to a pumped well in a fissured waterbearing formation. Journal of Hydrology, 35, pp. 257270. Boulton, N.S. and T.D. Streltsova, 1978. Unsteady flow to a pumped well in a fissured aquifer with a free surface level maintained constant. Water Resources Research, 14 (3), pp. 527532. Boulton, N .S. and T.D. streltsova, 1978. Unsteady flow to a pumped well in an unconfined fissured aquifer. Journal of Hydrology, 37, pp. 349363. Duffield, G.M. and J.O. Rumbaugh, 1989. Geraghty and Miller's AQTESOLV: Aquifer Test Solver. Geraghty and Miller Modeling Group, Reston, Virginia. Faulkner, G.L., 1973b. Geohydrology of the CrossFlorida Barge Canal Area with special reference to the Ocala vicinity: u.s. Geological Survey Water Resources Investigations Report 173, 117 p .. Freeze, R.A. and J.A. Cherry, 1979. Groundwater. Prentice Hall, Inc.
PAGE 74
62 Geotrans, Inc., 1988. Hydrologic investigation of the northern portion of the Southwest Florida Water Management District: Northern District Model Project, Phase Two Report. Unpublished report submitted to Southwest Water Management District. Geotrans, Inc., Herndon, Virginia. Geraghty and Miller, Inc., Sept. 1976. Aquifer performance test, Consumptive Use Permit Application: Phillips Petroleum. Unpublished report submitted to the south West Water Management District. Geraghty and Miller, Inc., Tampa, Florida. Guvanasen, v., D.S. Ward, P.F. Anderson, R. Schneider, P.N. Sims and L.A. Plante, 1987. Hydrologic investigation of the northern portion of the Southwest Florida Water Management District: Northern District Model Project, Phase One Report. Unpublished report submitted to Southwest Water Management District. Geotrans, Inc., Herndon, Virginia. Hantush, M.S., 1960. Modification of the theory of leaky aquifers. Journal of Geophysical Research, 65 (11), pp. 37133725. Hantush, M.S. and C.E. Jacob, 1955. Nonsteady radial flow in an infinite leaky aquifer. Am. Geophys. Union Trans., 36, pp. 95100. Jones, K.C., 1985. Northern Withlacoochee Hydrologic Investigation; Levy, Marion, Lake, and Alachua Counties, Florida. southwest Florida Water Management District, Florida. Kazemi, H., 1969. Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution. Soc. Pet. Eng. Jour., 9, pp. 451462. Long, J.c.s., J.S. Remer, C.R. Wilson, and P.A. Witherspoon, 1982. Porous media equivalents for networks of discontinuous fractures. Water Resources Research. 18 (3), pp. 645658. Maini Y.N.T., 1971. Insitu hydraulic parameters in measurement and interpretation. PhD dissertation, Imperial College, London, England.
PAGE 75
63 McDonald, M.G. and A.W. Harbough, 1984. A Modular ThreeDimensional Finite Difference Groundwater Flow Model: u.s. Geological Survey Open File Report 83875. Mercer, J.W. and C.R. Faust, 1980. Groundwater Modeling. National Water Well Association, Dublin, Ohio. Meinzer, O.E., 1927. Large springs in the United States. U.S. Geological Survey Water Supply Paper 557. Moench, D.M., 1984. Doubleporosity models for a fissured groundwater reservoir with fracture skin. Water Resources Research, 20 (7), pp. 831846. Moore, D.M. and M. Stewart, fracture zones in a karst aquifer. 61, pp. 325340. Geophysical signatures of Journal of Hydrology, Muskat, M., 1937. The flow of homogeneous fluids through porous media. J.W. Edwards, Inc. Ann Arbor, Mich .. 763 p. Neuman, S.P., 1975. Analysis of pumping test data from anisotropic unconfined aquifers considering delayed yield. Water Resources Research, 11 (2), pp. 329342. Parizak, R.R, 1976. On the significance of fracture traces and lineaments in carbonate and other terranes: in, Karst Hydrology and Water Resources, proceedings U.S.Yugos. Syrnp.1, Dubrovinik, 1975, Water Resources Publications. Pirson, S.J., 1953. Performance of fractured oil reservoirs. Bull. Arner. Assoc. Petrol. Geol., 37(2), pp. 232244. Pollard, P., 1959. Evaluation of acid treatments from pressure buildup analysis. Trans. Arner. Instit. Mining Eng., vol. 216, pp. 3843. Richard c. Fountain and Associates, January, 1975. Aquifer performance test, consumptive Use Permit Application: Beker Phosphate Company. Unpublished report for the South West Florida Water Management District. Richard C. Fountain and Associates, Tampa, Florida.
PAGE 76
64 Ryder, P.D., 1985. Hydrology of the Floridan Aquifer system in westcentral Florida: U.S. Geological Survey Professional Paper 1403F, 63 p. Sinclair, W.C., 1985. Preliminary evaluation of the watersupply potential of the springriver system in the WeekiWachee area and the lower Withlacoochee River, westcentral Florida: U.S. Geological Survey Water Resources Investigations Report 7874. Snow, D.T., 1965. A Parallel plate model of fractured permeable media, Ph.D dissertation, University of California, Berkeley, California. Stewart, M. and J. Wood, 1990. Geologic and geophysical character of fracture zones in a tertiary carbonate aquifer, Florida. In: Geotechnical and Environmental Geophysics, S.T. Ward (editor), Society of Exploration Geophysics, Investigations in Geophysics No. 5, vol. 2, pp. 235244. Streltsova, T.D., 1976. Hydrodynamics of groundwater flow in a fractured formation. Water Resources Research, 12(3), pp. 405414. Theis, C.V., 1935. The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage, Trans. Am. Geophys. Union. Annual Meeting, 16th, pp. 519524. van GolfRacht, 1982. Fundamentals of Fractured Reservoir Engineering. Elsevier, Amsterdam. Vernon, R.O., 1951. Geology of Citrus and Levy Counties, Florida. Florida Geological survey Bulletin 33, 256 p. wang, J.S.Y. and T.N. Narasimhan, 1985. Hydraulic mechanisms governing fluid flow in a partially saturated, fractured, porous medium. Water Resources Research, vol. 21, no. 12, pp. 18611874.
PAGE 77
65 Warren, J.E. and P.J. Root, 1963. The behavior of naturally fractured reservoirs. Soc. Pet. Eng. Jour., vol. 3, pp. 245255. Witherspoon, P.A., Y.W Tsang, J.C.S. Long, and J. Noorishad, 1981. New approaches to problems of fluid flow in fractured rock masses, Proc. 22nd U.S. Symp. Rock Mechs., Boston, Mass., pp. 120. Ward, D.S., D.C. Skipp, D. Griffin, and M.D. Barcelo, 1987. Dualporosity and discrete fracture simulation of groundwater flow in westcentral Florida, Proc. Fourth International Conference on the Use of Models,Assoc. of Ground Water Scientists and Engineers, Indianapolis, Indiana.
