U.S. Department of the Interior

U.S. Geological Survey

Open-File Report 2013–1086

Prepared in cooperation with the Miami-Dade County Water and Sewer Department

and Department of Regulatory and Economic Resources

Estimation of Capture Zones and Drawdown at the

Northwest and West Well Fields, Miami-Dade County,

Florida, Using an Unconstrained Monte Carlo Analysis:

Recent (2004) and Proposed Conditions

i

Prepared in cooperation with the Miami-Dade County Water and Sewer Department 

and Department of Regulatory and Economic Resources 

 

Estimation of Capture Zones and Drawdown at the 

Northwest and West Well Fields, Miami-Dade County, 

Florida, Using an Unconstrained Monte Carlo Analysis: 

Recent (2004) and Proposed Conditions 

By Linzy Brakefield, Joseph D. Hughes, Christian D. Langevin, and Kevin Chartier

Open-File Report 2013–1086

U.S. Department of the Interior 

U.S. Geological Survey
ii

U.S. Department of the Interior 

KEN SALAZAR, Secretary 

U.S. Geological Survey 

Suzette M. Kimball, Acting Director

U.S. Geological Survey, Reston, Virginia: 2013 

For more information on the USGS—the Federal source for science about the Earth,

its natural and living resources, natural hazards, and the environment—visit 

http://www.usgs.gov or call 1–888–ASK–USGS

For an overview of USGS information products, including maps, imagery, and publications,

visit http://www.usgs.gov/pubprod

To order this and other USGS information products, visit http://store.usgs.gov

Suggested citation:

Brakefield, L., Hughes, J.D., Langevin, C.D., and Chartier, K., 2013, Estimation of capture zones and 

drawdown at the Northwest and West Well Fields, Miami-Dade County, Florida, using an unconstrained 

Monte Carlo analysis: recent (2004) and proposed conditions: U.S. Geological Survey Open-File 

Report 2013–1086, 124 p.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply 

endorsement by the U.S. Government. 

Although this information product, for the most part, is in the public domain, it also may contain

copyrighted materials as noted in the text. Permission to reproduce copyrighted items must be 

secured from the copyright owner. 
iii

Acknowledgments 

The authors would like to acknowledge assistance from Miami-Dade County Water and Sewer 

Department and the Department of Regulatory and Economic Resources employees Virginia Walsh, 

Maria MacFarland, Wilbur Mayorga and Rashid Istambouli for their assistance and support with data and 

review of the analyses and report. The authors would like to thank Dave Welter of the South Florida 

Water Management District for his insightful reviews, which led to improvement of the manuscript. 

The authors would also like to acknowledge the helpful technical reviews of Kevin Cunningham 

and Jeffrey Starn of the U.S. Geological Survey, which greatly improved the clarity and thoroughness of 

the report. Conversations with Kevin Cunningham, in particular, greatly strengthened the Monte Carlo 

analysis presented in Appendix 1. Appreciation is also extended to Dorothy Payne, Robert Renken, Eve 

Kuniansky, Michael Deacon and Kim Swidarski for their editorial and technical reviews, which greatly 

improved the flow, organization, and clarity of the report. 

iv

Contents 

Acknowledgments.............................................................................................................................................................iii

Contents ..............................................................................................................................................................................iv

Abstract ................................................................................................................................................................................1

Introduction .........................................................................................................................................................................2

Purpose and Scope ........................................................................................................................................................3

Approach ..........................................................................................................................................................................3

Previous Studies .............................................................................................................................................................4

Description of Study Area.................................................................................................................................................5

Hydrostratigraphy ...........................................................................................................................................................6

Aquifer Properties ..........................................................................................................................................................6

Water-Budget Components..........................................................................................................................................8

Ground-Water Withdrawals .....................................................................................................................................8

Surface-Water/Groundwater Interaction ..............................................................................................................8

Rainfall, Evapotranspiration, and Runoff ................................................................................................................9

Simulation of Groundwater Flow near the Northwest and West Well Fields, 1996–2004 ...................................10

Model Construction......................................................................................................................................................10

Spatial and Temporal Discretization .....................................................................................................................10

Assignment of Aquifer Parameters .......................................................................................................................12

Boundary Conditions ................................................................................................................................................12

Biscayne Bay .........................................................................................................................................................13

Groundwater Withdrawals for Public Supply ..................................................................................................13

Canals ......................................................................................................................................................................13

Recharge ................................................................................................................................................................15

Evapotranspiration ................................................................................................................................................15

Quarry Lakes ..........................................................................................................................................................15

Initial Conditions........................................................................................................................................................16

Model Fit.........................................................................................................................................................................16

Comparison of Simulated and Observed Water Levels .....................................................................................17

Comparison of Simulated and Observed Net Exchange between Canals and the Aquifer ........................18

Comparison of Dry and Wet Season Groundwater-Level Altitudes in Miami-Dade County.......................19

Estimation of Well-Field Capture Zones and Drawdown ..........................................................................................19

Methods..........................................................................................................................................................................19

Lake Representation—Bulk Hydraulic Conductivity Calculation.....................................................................20

Determination of Drawdown and Travel-Time Capture Zones.........................................................................21

Evaluation of Uncertainty Using Monte Carlo Analysis.....................................................................................22

Accounting for Uncertainty in Spatial Distribution of Hydraulic Conductivity ..........................................22

Accounting for Uncertainty in Effective Porosity............................................................................................24

Selection of Stress Period for Steady-State Particle-Tracking Simulation ...................................................24

Simulated Time of Travel and Drawdown for Dry, Average, and Wet Hydrologic Conditions.......................25

Base-Case Deterministic Simulations ..................................................................................................................25

Monte Carlo Simulations .........................................................................................................................................25
v

Two Proposed Scenarios: Potential Effects of Seepage Mitigation and Lake Expansion on Capture Zones 

and Drawdown..................................................................................................................................................................29

Proposed Scenario 1: Investigation of Effects of Seepage Mitigation...............................................................29

Proposed Scenario 2: Investigation of Effects of Lake Expansion ......................................................................30

Model Limitations .............................................................................................................................................................31

Summary and Conclusions .............................................................................................................................................34

Selected References .......................................................................................................................................................37

Appendices ........................................................................................................................................................................43

Figures 

Figure 1. Locations of all base maps used throughout this report. Base maps are referred to in each map 

figure caption ........................................................................................................................................... 58

Figure 2. Model grid extent, study area, hydrogeologic sections A–A′ and B–B′, location of production well 

fields and location of canals, southeastern Florida ................................................................................. 59

Figure 3. Study area and existing Wellhead Protection Areas for production well fields in 

Miami-Dade County ................................................................................................................................. 60

Figure 4. Hydrogeologic cross section with subset section A–A′ ............................................................................ 61

Figure 5. Hydrogeologic cross section B–B′ across central part of the study area, including the 

Northwest Well Field................................................................................................................................ 62

Figure 6. Altitude of the base of the Biscayne aquifer and its western extent throughout most of 

Miami-Dade County ................................................................................................................................. 63

Figure 7. Daily municipal supply well withdrawals from the Biscayne aquifer in Miami-Dade County 

for January 1, 1996, through December 31, 2004 ................................................................................... 64

Figure 8. Primary and secondary canal system in study area................................................................................. 65

Figure 9. Location of active and inactive model cells, and type of MODFLOW boundary condition........................ 66

Figure 10. Altitude of the top of the Biscayne aquifer model................................................................................. 67

Figure 11. Hydraulic conductivity field for the base-case simulation and locations of transmissivity 

estimates derived from aquifer performance tests ............................................................................... 68

Figure 12. Locations of lakes in 1999 (transient simulation) and 2004 ................................................................. 69

Figure 13. Daily observed and model-simulated water levels for selected monitoring wells in the study area...... 70

Figure 14. Locations of groundwater observation data used to calibrate model and mean error calculated 

from observed and simulated daily water levels .................................................................................. 71

Figure 15. Selected surface-water basins (shaded) in simulated area, locations of production wells, 

and canals ........................................................................................................................................... 72

Figure 16. Monthly model-simulated and observed net exchange for composite C-4 and C-6 basins, 

and mean error, mean absolute error and the coefficients of determination for simulated versus 

observed net exchange for each basin................................................................................................ 73

Figure 17. Average water-level differences for the Biscayne aquifer between May 1993 and May 1998 ............. 74

Figure 18. Average water-level differences for the Biscayne aquifer between November 1993 and 

November 1998 ................................................................................................................................... 75

Figure 19. Hydraulic conductivity field used in transient and deterministic base-case steady-state 

simulations in the study area for the base-case simulation, and for three selected 

additional realizations .......................................................................................................................... 76
vi

Figure 20. Log-normal probability distribution of effective porosity used for sampling of random 

values for Monte Carlo analysis........................................................................................................... 77

Figure 21. Tailwater stages at S-337 and G-119 with selected days representing dry, average, and 

wet hydrologic conditions identified with black circles.......................................................................... 78

Figure 22. Simulated 0.25-foot and 0.1-foot drawdown contours, and particle endpoints for travel times 

of 10, 30, 100, and 210 days for the base-case simulation of dry conditions....................................... 79

Figure 23. Simulated 0.25-foot and 0.1-foot drawdown contours, and particle endpoints for travel

times of 10, 30, 100, and 210 days for the base-case simulation of average conditions .................... 80

Figure 24. Simulated 0.25-foot and 0.1-foot drawdown contours, and particle endpoints for travel 

times of 10, 30, 100, and 210 days for the base-case simulation of wet conditions............................. 81

Figure 25. Simulated particle endpoint locations for travel times of 10, 30, 100, and 210 days for 

dry conditions for 10,000 realizations .................................................................................................. 82

Figure 26. Simulated particle endpoint locations for 210 days of travel for dry conditions for 

100 realizations around the Northwest Well Field, and an example of the radial class division 

described by Varljen and Shafer.......................................................................................................... 83

Figure 27. Median and 95-percent confidence intervals for the 210-day capture zones for the 

Northwest Well Field and West Well Field for steady-state, dry conditions, and particle 

endpoints for 210 days of travel for the base-case simulation of dry conditions .................................. 84

Figure 28. Median and 95-percent confidence intervals for the 100-day capture zones for the 

Northwest Well Field and West Well Field for steady-state, dry conditions, and 

particle endpoints for 100 days of travel for the base-case simulation of dry conditions...................... 85

Figure 29. Median and 95-percent confidence intervals for the 30-day capture zones for the 

Northwest Well Field and West Well Field for steady-state dry conditions, and particle 

endpoints for 30 days of travel for the base-case simulation of dry conditions .................................... 86

Figure 30. Median and 95-percent confidence intervals for the 10-day capture zones for the 

Northwest Well Field and West Well Field for steady-state, dry conditions, and 

particle endpoints for 10 days of travel for the base-case simulation of dry conditions........................ 87

Figure 31. Median and 95-percent confidence intervals for the 210-day capture zones for the 

Northwest Well Field and West Well Field or steady-state, average conditions, and 

particle endpoints for 210 days of travel for the base-case simulation of average conditions.............. 88

Figure 32. Median and 95-percent confidence intervals for the 100-day capture zones for the 

Northwest Well Field and West Well Field for steady-state, average conditions, and 

particle endpoints for 100 days of travel for the base-case simulation of average conditions.............. 89

Figure 33. Median and 95-percent confidence intervals for the 30-day capture zones for the 

Northwest Well Field and West Well Field for steady-state, average conditions, and 

particle endpoints for 30 days of travel for the base-case simulation of average conditions................ 90

Figure 34. Median and 95-percent confidence intervals for the 10-day capture zones for the 

Northwest Well Field and West Well Field for steady-state, average conditions, and 

particle endpoints for 10 days of travel for the base-case simulation of average conditions................ 91

Figure 35. Median and 95-percent confidence intervals for the 210-day capture zones for the 

Northwest Well Field (NWWF) and West Well Field (WWF) for steady-state, wet conditions, 

and particle endpoints for 210 days of travel for the base-case simulation of wet conditions .............. 92

Figure 36. Median and 95-percent confidence intervals for the 100-day capture zones for the 

Northwest Well Field and West Well Field for steady-state, wet conditions, and particle 

endpoints for 100 days of travel for the base-case simulation of wet conditions.................................. 93

Figure 37. Median and 95-percent confidence intervals for the 30-day capture zones for the 

Northwest Well Field and West Well Field for steady-state, wet conditions, and particle 

endpoints for 30 days of travel for the base-case simulation of wet conditions.................................... 94
vii

Figure 38. Median and 95-percent confidence intervals for the 10-day capture zones for the 

Northwest Well Field and West Well Field for steady-state, wet conditions, and particle 

endpoints for 10 days of travel for the base-case simulation of wet conditions.................................... 95

Figure 39. Upper bounds of the 95-percent confidence intervals (CI) for the 210-day capture zones for 

the Northwest Well Field and West Well Field for steady-state, wet and dry conditions ...................... 96

Figure 40. Composite areal extent of the upper bounds of the 95-percent confidence intervals for 

the 210-day capture zones for the Northwest Well Field and West Well Field for 

steady-state, wet and dry conditions.................................................................................................... 97

Figure 41. Median and 95-percent confidence intervals for 0.1-foot and 0.25-foot drawdowns for 

the Northwest Well Field and West Well Field for steady-state, dry conditions.................................... 98

Figure 42. Median and 95-percent confidence intervals for 0.1-foot and 0.25-foot drawdowns for 

the Northwest Well Field and West Well Field for steady-state, average conditions............................ 99

Figure 43. Median and 95-percent confidence intervals for 0.1-foot and 0.25-foot drawdowns for 

the Northwest Well Field and West Well Field for steady-state, wet conditions ................................. 100

Figure 44. Locations of the recharge canal west of the Northwest Well Field, the proposed slurry 

wall between the L-31N Levee and the L-31N Canal, and the Horizontal Flow Barriers 

west of the West Well Field for Scenario 1 ........................................................................................ 101

Figure 45. Median and 95-percent confidence intervals for the 210-day capture zones for the 

Northwest Well Field for steady-state, dry conditions for seepage mitigation scenario 

and without seepage mitigation, and final particle endpoints for both conditions beyond 

210 days ............................................................................................................................................ 102

Figure 46. Upper bound of the 95-percent confidence intervals for the 210-day capture zones and 

upper bound of the 95-percent CI for the 0.1-foot drawdown for the West Well Field for 

steady-state, dry conditions for seepage mitigation scenario and without seepage mitigation ......... 103

Figure 47. Upper bound of the 95-percent confidence intervals for the 100-day capture zones for 

the West Well Field (WWF) for steady-state, dry conditions for seepage mitigation 

scenario and without seepage mitigation........................................................................................... 104

Figure 48. Locations of existing lakes, projected expansion of existing lakes and new proposed lakes 

for lake-expansion scenario............................................................................................................... 105

Figure 49. Median and upper bound of the 95-percent confidence intervals for the 210-day capture 

zones and upper bound of the 95-percent CI for the 0.25-foot drawdown for the Northwest 

Well Field for steady-state dry conditions for the lake expansion scenario and without 

lake expansion................................................................................................................................... 106

Figure 50. Upper bound of the 95-percent confidence intervals for the 210-day capture zones and 

upper bound of the 95-percent CI for the 0.1-foot drawdown for the West Well Field for 

steady-state, dry conditions for the lake expansion scenario and without lake expansion................. 107

Figure 1-1. Simulated 500-day particle pathlines and head contours for one-layer conceptual Model A 

and three-layer conceptual model B .................................................................................................. 108 

Figure 1-2. Model grid, locations of constant-head cells, river cells and pumping wells for Model A and

Model B ............................................................................................................................................. 109

Figure 1-3. Median and 95-percent confidence intervals for 1-foot drawdown for Model A and

Model B ............................................................................................................................................. 110

Figure 1-4. Median and upper bound 95-percent confidence intervals for 210-day capture zones 

for Model A and Model B ................................................................................................................... 111
viii

Figure 3-1. Simulated water levels for base-case, steady-state average conditions with model 

boundary as described in the main body of the report, with model boundary extended

west of the western canals and general head boundary conductance of 269 square

feet per day, and with model boundary extended west of the western canals and

general head boundary conductance of 26,910 ft2/d ......................................................................... 112 

Figure 4-1. Simulated steady-state groundwater head and backward-tracked particle pathlines

for the model representing lakes with high hydraulic conductivity and the model

using the LAK Package to represent Lakes....................................................................................... 113 

Figure 5-1. Daily observed and model-simulated water levels for monitoring wells in the study

area not shown in figure 13................................................................................................................ 114

Figure 6-1. Backward-tracked particle pathlines from a pumping well through lakes without

residence-time post-processing and with residence-time post-processing ........................................ 117

Figure 6-2. Map showing locations and numbering of lakes referred to in paper by Chin and

others (2010) ..................................................................................................................................... 118 

Tables 

Table 1. Monitoring well names, mean errors of head, mean absolute errors of head, and the 

coefficients of determination of simulated vs. observed values of head................................................. 120

Table 2. Statistics of minimum value and maximum value of hydraulic conductivity for all realizations, 

prior to calculation of bulk hydraulic conductivities................................................................................. 121

Table 4-1. Simulated groundwater inflow, groundwater outflow, and evaporation lake water-budget

terms for lakes simulated using the high-hydraulic-conductivity approach and the 

LAK Package......................................................................................................................................... 122 

Table 6-1. Residence times used in model documented in this report and those calculated by

Chin and others (2010)........................................................................................................................... 124
ix

Conversion Factors, Abbreviations, and Datums 

Multiply By To obtain

Length

inch (in.) 2.54 centimeter (cm)

inch (in.) 25.4 millimeter (mm)

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Area

square mile (mi2) 2.590 square kilometer (km2) 

Flow rate

foot per day (ft/d) 0.3048 meter per day (m/d)

cubic foot per day (ft3/d) 0.02832 cubic meter per day (m3/d)

million gallons per day (Mgal/d) 0.04381 cubic meter per second (m3/s)

inch per year (in/yr) 25.4 millimeter per year (mm/yr)

Rainfall rate

inch per day (in/d) 0 .0254 meter per day (m/d)

Hydraulic conductivity

foot per day (ft/d) 0.3048 meter per day (m/d)

Transmissivity*

foot squared per day (ft2/d) 0.09290 meter squared per day (m2/d) 

Leakance

foot per day per foot [(ft/d)/ft] 1 meter per day per meter

inch per year per foot [(in/yr)/ft] 83.33 millimeter per year per meter

[(mm/yr)/m]

APT Aquifer performance test

CDM Camp, Dresser & McKee, Inc.

CI Confidence interval

ENP Everglades National Park

HFB Horizontal flow barrier package of MODFLOW

MAE Mean absolute error

ME Mean error

NWWF Northwest Well Field 

RER Regulatory and Economic Resources (Miami-Dade County)

SGSIM Sequential Gaussian simulation method

SFWMD South Florida Water Management District

WHPA Wellhead protection area

WWF West Well Field
x

Altitude is distance above the vertical datum. 

Vertical coordinate information is referenced to the North American Vertical Datum of 1988 (NAVD 88) unless stated to be 

referenced to the National Geodetic Vertical Datum of 1929 (NGVD 29).

Horizontal coordinate information is referenced to the North American Datum of 1983 (NAD 83).

*Transmissivity: The standard unit for transmissivity is cubic foot per day per square foot times foot of aquifer thickness 

[(ft3

/d)/ft2]ft. In this report, the mathematically reduced form, foot squared per day (ft2/d), is used for convenience.
1

Estimation of Capture Zones and Drawdown at the 

Northwest and West Well Fields, Miami-Dade County, 

Florida, Using an Unconstrained Monte Carlo 

Analysis: Recent (2004) and Proposed Conditions 

By Linzy Brakefield, Joseph D. Hughes, Christian D. Langevin, and Kevin Chartier

Abstract 

Travel-time capture zones and drawdown for two production well fields, used for drinking-water 

supply in Miami-Dade County, southeastern Florida, were delineated by the U.S Geological Survey 

using an unconstrained Monte Carlo analysis. The well fields, designed to supply a combined total of 

approximately 250 million gallons of water per day, pump from the highly transmissive Biscayne 

aquifer in the urban corridor between the Everglades and Biscayne Bay. A transient groundwater flow

model was developed and calibrated to field data to ensure an acceptable match between simulated and 

observed values for aquifer heads and net exchange of water between the aquifer and canals. Steadystate conditions were imposed on the transient model and a post-processing backward particle-tracking 

approach was implemented. Multiple stochastic realizations of horizontal hydraulic conductivity, 

conductance of canals, and effective porosity were simulated for steady-state conditions representative 

of dry, average and wet hydrologic conditions to calculate travel-time capture zones of potential source 

areas of the well fields. Quarry lakes, formed as a product of rock-mining activities, whose effects have 

previously not been considered in estimation of capture zones, were represented using high hydraulicconductivity, high-porosity cells, with the bulk hydraulic conductivity of each cell calculated based on 

estimates of aquifer hydraulic conductivity, lake depths and aquifer thicknesses. A post-processing 

adjustment, based on calculated residence times using lake outflows and known lake volumes, was 

utilized to adjust particle endpoints to account for an estimate of residence-time-based mixing of lakes. 

Drawdown contours of 0.1 and 0.25 foot were delineated for the dry, average, and wet hydrologic 

conditions as well. In addition, 95-percent confidence intervals (CIs) were calculated for the capture 

zones and drawdown contours to delineate a zone of uncertainty about the median estimates. 

Results of the Monte Carlo simulations indicate particle travel distances at the Northwest Well 

Field (NWWF) and West Well Field (WWF) are greatest to the west, towards the Everglades. The manmade quarry lakes substantially affect particle travel distances. In general near the NWWF, the capture 

zones in areas with lakes were smaller in areal extent than capture zones in areas without lakes. It is 

possible that contamination could reach the well fields quickly, within 10 days in some cases, if it were 

introduced into lakes nearest to supply wells, with one of the lakes being only approximately 650 feet 

from the nearest supply well. 
2

In addition to estimating drawdown and travel-time capture zones of 10, 30, 100, and 210 days 

for the NWWF and the WWF under more recent conditions, two proposed scenarios were evaluated 

with Monte Carlo simulations: the potential hydrologic effects of proposed Everglades groundwater 

seepage mitigation and quarry lake expansion. The seepage mitigation scenario included the addition of 

two proposed anthropogenic features to the model: (1) an impermeable horizontal flow barrier east of 

the L-31N canal along the western model boundary between the Everglades and the urban areas of 

Miami-Dade County, and (2) a recharge canal along the Dade-Broward Levee near the NWWF. Capture 

zones and drawdown for the WWF were substantially affected by the addition of the barrier, which 

eliminates flow from the western boundary into the active model domain, shifting the predominant 

capture zone source area from the west more to the north and south. The 95-percent CIs for the 210-day 

capture zone moved slightly in the NWWF as a result of the recharge canal. The lake-expansion 

scenario incorporated a proposed increase in the number and surface area of lakes by an additional 

25 square miles. This scenario represents a150-percent increase from the 2004 lake surface area near 

both well fields, but with the majority of increase proposed near the NWWF. The lake-expansion 

scenario substantially decreased the extent of the 210-day capture zone of the NWWF, which is limited 

to the lakes nearest the well field under proposed conditions. 

Introduction 

The karstic Biscayne aquifer in Miami-Dade County has high permeability, is shallow and 

minimally confined, and thus, there are strong interactions between surface-water features and the 

aquifer, making the aquifer sensitive to activities at the surface. Regional groundwater flow across the 

aquifer beneath Miami-Dade County (fig. 1) is generally from the wetlands in the western part of the 

county toward Biscayne Bay to the southeast. Groundwater in the aquifer in the north-central part of the 

county flows predominantly through stratiform, high-permeability flow zones characterized by 

macroporosity features that provide pathways for rapid transport of groundwater, and result in some of 

the highest transmissivity values ever recorded. The county pumps approximately 400 million gallons 

per day (Mgal/d) from the aquifer, and it is a designated sole source of drinking water (Federal Register 

Notice, 1979). Some public supply well fields that pump from the Biscayne aquifer are in proximity to 

canals that were built to control inundation and aquifer recharge, and to quarry lakes, and there is 

concern about the risk that contamination of these surface-water features poses to the Biscayne aquifer. 

To better protect the groundwater supply, the U.S. Geological Survey, in cooperation with Miami-Dade 

County, initiated in 2008 a hydrologic investigation of the Biscayne aquifer to evaluate the interactions 

between groundwater and surface water.

Of particular concern for Miami-Dade County are the Northwest Well Field (NWWF) and West 

Well Field (WWF), located in the northwestern part of the urbanized part of the county (fig. 2). The 

NWWF is the largest supply well field for drinking water purposes in Miami-Dade County and is 

designed to withdraw 225 Mgal/d of water. The WWF is located southwest of the NWWF and has a 

capacity to withdraw 25 Mgal/d. Both well fields are near canals, and the NWWF is near numerous 

quarry lakes. Several proposed construction projects may affect groundwater flow fields at and around 

the well fields. The county has established “Wellhead protection areas” (WHPAs) (fig. 3) around the 

well fields to limit activities that may adversely affect the productivity of the well fields and the 

groundwater quality. The county uses results of groundwater flow models, including estimated capture 

zones and drawdown, to establish the WHPAs. There is considerable uncertainty in these estimations, 

and quantification of the uncertainty may assist the county in protection of groundwater that supplies the 

well fields. Furthermore, evaluating the uncertainty in hydraulic properties, canal exchanges, and lake 
3

residence times may assist the county in understanding the risk to the groundwater supply from 

contamination of local surface and subsurface areas near these well fields. 

Purpose and Scope

This report documents estimates of well-field capture zones and drawdown using a stochastic 

approach near the NWWF and WWF based on current understanding of the local hydrogeology and 

hydrology. To estimate well-field capture zones and drawdown, a groundwater model was constructed 

that allows for an evaluation of uncertainty in the results based on known uncertainties in the 

conceptualization of the flow system, groundwater exchanges with the canals and lakes, and hydraulic 

parameters. This report documents the development of the flow model for transient and steady-state 

conditions, including an unconstrained Monte Carlo analysis used with the steady-state simulations to 

estimate time-based capture zones and drawdown contours with 95-percent confidence intervals. The 

model was also used to investigate effects on capture zones and drawdown by proposed anthropogenic 

alterations to the hydrologic system. One scenario considers the effects of seepage mitigation, including 

a barrier to the groundwater flux from the west of the WWF, and implementation of a recharge canal 

west of the NWWF. A second scenario highlights potential effects from proposed lake expansion in the 

Lake Belt area.

Approach

Hydrologic analyses based on groundwater flow models have been used for delineating the areas 

that provide recharge to a supply well (Starn and Brown, 2007; Starn and others, 2000; Franke and 

others, 1998; Snyder and others, 1998; Barlow, 1995; Reilly and Pollock, 1993). Hydrologic analyses 

are often used to determine the impact of a municipal supply well on the surrounding aquifer and the 

spatial extent of this impact. One common method is to delineate travel-time capture zones. Capture 

zones provide insight into the lateral extent of the recharge area for a supply well, but because they are 

defined by specified times of travel, the zones may not necessarily encompass the entire contributing 

area for a given supply well. Additionally, drawdown, defined in this study as the difference in waterlevel altitude caused solely by well-field withdrawals, can be used to indicate how a well field will 

affect groundwater flow patterns. In this investigation, travel-time capture zones of 10, 30, 100, and 

210 days were delineated, as well as 0.25-foot (ft) and 0.1-ft drawdown contours. 

Groundwater models that provide estimates of capture zones and drawdown typically use a 

deterministic approach in which a single result is produced based from a single set of model input 

parameters. Uncertainty in simulated results cannot be quantified with this approach. The hydrologic 

analysis presented herein uses an unconstrained Monte Carlo analysis to delineate the 95-percent

confidence interval (CI) associated with simulated capture zones and drawdown contours. With the 

Monte Carlo analysis, capture zone and drawdown contours are calculated for many different parameter 

sets that are equally plausible, in a statistical sense. Each parameter set and associated model result is 

called a realization. The approach is considered unconstrained in that no attempt was made to calibrate 

each realization. Reference to this approach is simply referred to as “Monte Carlo analysis” throughout 

the report, although in all instances the Monte Carlo simulations are unconstrained. Uncertainty in 

capture zones and drawdown can be assessed by evaluating the ensemble, which is the collection of all 

realizations. Monte Carlo methods have been used previously for these types of analyses (Varljen and 

Shafer, 1991; Bair and others, 1991; Franzetti and Guadagnini, 1996; Starn and others, 2000; Starn and 

Brown, 2007; Starn and others, 2010). 
4

The numerical groundwater flow model used for this study was calibrated for a 9-year 

simulation period, and then steady-state simulations of dry, average, and wet conditions were used along 

with the Monte Carlo analysis to estimate the median and 95-percent CI for both capture zones and 

drawdowns. Well-field withdrawal rates used in the steady-state simulations were based on updated 

design capacities provided by Miami-Dade County and differ from those used in previous analyses. The 

hydraulic effect of quarry lakes is a prominent component of the steady-state simulations presented 

herein. The unconstrained Monte Carlo analysis was also applied to proposed hydrologic scenarios in 

order to demonstrate how well-field capture zones and drawdown might be affected by possible 

alterations to the existing water management system in Miami-Dade County. The approach and methods 

discussed herein are also transferable to other studies that involve estimation of well-field capture zones 

and drawdown. 

Previous Studies 

The existing Wellhead Protection Areas (WHPAs) for the NWWF (fig. 3) were delineated using 

results from a one-layer groundwater model developed by Camp Dresser & McKee Inc. (CDM). The 

CDM model was constructed using a modified version of the Prickett-Lonnquist Aquifer Simulation 

Model (PLASM) (Camp Dresser & McKee, 1985; Prickett and Lonnquist, 1971). The PLASM code is 

similar to MODFLOW in that it is based on a finite-difference approximation of the groundwater flow 

equation. Lakes were not represented in the CDM model and the NWWF withdrawal rate was set at 

220 Mgal/d (Camp Dresser & McKee, 1985). The effects of lakes are included in the present analysis 

and the maximum design-capacity pumping rate for the NWWF is 225 Mgal/d. 

The existing WHPAs for the WWF (fig. 3), based on travel-time capture zones and drawdown 

delineated by a separate groundwater flow model than the model used for the NWWF, are based on the 

projected 1992 withdrawal rate of 140 Mgal/d (Department of Environmental Resources Management, 

1992). Presently, the withdrawal rate for the WWF is only 25 Mgal/d from three wells in the southern 

part of the well field. Therefore, results from the numerical simulations described herein are not directly 

comparable to the existing capture zone and drawdown delineations of the WWF (fig. 3).

In 2001, CH2M Hill developed a MODFLOW model and used particle tracking to delineate 

capture zones of 180 and 230 days of travel to determine rock-mining setback distances for ongoing 

mining near the NWWF (CH2M Hill, 2001). The NWWF withdrawal rate used in the CH2M Hill model 

was 235 Mgal/d. Effects of lakes on the capture zones were investigated with the model by representing 

the lakes as high-porosity and high hydraulic-conductivity zones having a depth of 70 ft below land 

surface. The CH2M Hill model was based on the North Miami-Dade County Sub-Regional 

Groundwater Flow Model developed by the South Florida Water Management District (SFWMD). The 

CH2M Hill model contained multiple layers and incorporated more accurate, updated estimates of 

Biscayne aquifer transmissivities compared to those used in the 1985 CDM model. The travel-time 

capture zones of the CH2MHill model results were more spatially extensive west of the NWWF because 

of the updated transmissivities, presence of lakes in the eastern part of the well field, and a greater wellfield withdrawal rate. Although the CH2M Hill model provides a more updated assessment of 

groundwater flow and times of travel near the well field than offered by the 1985 CDM model, the 

existing Miami-Dade WHPA maps are based on capture zone and drawdown results from the 1985 

CDM model. 

Since the determination of the existing WHPAs for the NWWF and WWF, more recent 

subsurface investigations in the area have characterized and mapped highly porous, stratiform, flow 
5

zones in the Biscayne aquifer that may play an important role in groundwater flow and contaminant 

transport in the study area (Cunningham and others, 2006a, 2009; Renken and others, 2005, 2008). 

These investigations focused primarily on the upper part of the Tamiami Formation, Fort Thompson 

Formation, and Miami Limestone. The purpose of these investigations was to characterize the lithologic 

variations, cyclostratigraphy, aquifer properties, and identify groundwater flow zones in the karstic 

limestone that composes the Biscayne aquifer. The increased potential for risk of contamination because 

of these flow zones near supply wells was investigated through the use of multiple tracer tests in 2003 

and 2004 near the NWWF (Renken and others, 2005, 2008; Shapiro and others, 2008; Harvey and 

others, 2008). Results indicated that flow near the NWWF is dominated by megaporosity flow zones, 

and that the potential transport of contaminants could be more rapid near this well field than previously 

thought. Shoemaker and others (2008a) used a version of MODFLOW capable of accounting for 

turbulent groundwater flow (MODFLOW-Conduit Flow Process (CFP) – Shoemaker and others, 2008b) 

and a multi-layer model to investigate the potential for (1) turbulent flow in highly permeable flow 

zones of the Biscayne aquifer in the Lake Belt area and (2) the effects of turbulent flow on hydraulic 

heads and parameter sensitivities. 

Guha (2008) developed a multi-layer model with stochastic MODPATH 5.0 particle tracking of 

the NWWF area to estimate a 60-day travel-time capture zone for WHPA refinement purposes. The 

Monte Carlo analysis used 300 spatially-constant realizations sampling a uniform distribution of 

hydraulic conductivity for an assumed occurrence of a single flow zone in the Biscayne aquifer. Lakes 

were modeled as high hydraulic conductivity cells using a hydraulic conductivity of 1.00 × 106 feet per 

day (ft/d) (305 km/d). Chin and others (2010) used an advective particle-tracking method with a 

random-walk dispersive component to estimate travel-time capture zones in locations containing lakes 

using a residence-time approach. This method was applied to Guha’s (2008) model to estimate capture 

zones for 10, 30, 100, and 210 days of travel. This approach is compared to the approach utilized in the 

present investigation and discussed in this report. 

Description of Study Area 

The study area is in the north-central part of Miami-Dade County and includes the NWWF and 

WWF (fig. 2). Much of the study area consists of protected wetlands (in particular the Pennsuco 

Wetlands, fig. 3), borrow-pit rock mines, upland forests, and agricultural lands. The southeastern 

portion of the study area consists of residential and commercial areas. The study area lies just east of the 

protected wetlands of Everglades National Park (ENP) and Water Conservation Area 3B. The study area 

includes the western edge of the urbanized part of Miami-Dade County. The study area is used not only 

for limestone mining, but it also serves as an important source of drinking water supply, and as a 

development buffer area between the urban parts of the county and ENP and Water Conservation Area 

3B. The study area is dominated by the 89-square-mile (mi2

) Lake Belt area (figs. 2 and 3), named for 

its characteristic presence of extensive quarry “lakes” formed through rock-mining activities. These 

man-made water bodies are the major open-water features in the study area and are referred to simply as 

lakes in the remainder of this report.
6

Hydrostratigraphy

The focus of this investigation is the Biscayne aquifer, the primary source of municipal water 

supply for Miami-Dade County and a designated sole source of drinking water (Federal Register Notice, 

1979) (figs. 4 and 5). The Biscayne aquifer is the highly permeable, uppermost part of the surficial 

aquifer system (Parker, 1951), which is composed of lithostratigraphic units of varying permeability 

from land surface downward to the top of the intermediate confining unit, at a depth between about 160 

and 200 ft below NGVD 29 (Fish and Stewart, 1991; Reese and Cunningham, 2000). In this 

investigation, the Biscayne aquifer is defined as that part of the surficial aquifer system composed of the 

modern peats, Miami Limestone, Fort Thompson Formation, and the Tamiami Formation where it is 

sufficiently permeable (based on borehole geophysical logs and flowmeter measurements published by 

Cunningham and others, 2004a, 2006b) (figs. 4 and 5). The Biscayne aquifer is underlain by a confining 

to semiconfining unit. 

Parker (1951) indicated that the Biscayne aquifer is unconfined. Fish and Stewart (1991) 

clarified this by noting that over most of the area “(except near well fields or margins of waterconservation areas), water levels at depth are almost identical to local water-table elevations.” Fish and 

Stewart (1991) observed that variations in the vertical distribution of water-level responses, during the 

early stages of pumping in wells open to highly permeable zones overlain by less permeable sands, may 

indicate semiconfined behavior. For well fields along the margins of water conservation areas, Fish and 

Stewart (1991) also noted that water levels measured at depth may differ from those observed at the 

water table. For example, near the Levee 30 Canal west of the NWWF, Sonenshein (2001) mapped 

vertical water-level differences on the order of 0.1 ft or less. 

In Miami-Dade County, the Biscayne aquifer is thickest at the coast and thins toward the west 

(fig. 4). In the study area, the base of the Biscayne aquifer generally dips eastward (fig. 6). The areal 

extent, bottom altitudes, and characteristics of the Biscayne aquifer are discussed at length in Fish and 

Stewart (1991). This highly transmissive aquifer, in the study area, comprises multiple vertical 

lithofacies successions (Cunningham and others, 2006a, b) arranged within high-frequency cycles 

(fig. 5). Vertical lithofacies successions record the prevailing water depth during a single cycle of 

relative rise and fall in sea level (Kerans and Tinker, 1997), and within the Biscayne aquifer in the study 

area, porosity and permeability are typically linked to the vertical arrangement of lithofacies within each 

succession that composes a high-frequency cycle (Cunningham and others, 2006a, b). In particular, 

megaporosity (also referred to as “touching vug porosity”), which is associated with highly permeable 

flow zones in the Biscayne aquifer, correlates with cyclostratigraphic sequences across the study area 

(fig. 5). 

Highly porous, groundwater flow zones have been characterized and mapped throughout much 

of the study area and identified in other parts of the Biscayne aquifer in Miami-Dade County 

(Cunningham and others, 2004a, b, 2006a, b, 2009; Renken and others, 2005, 2008; Shapiro and others, 

2008; Cunningham and Sukop, 2011). The cumulative thickness of the flow zones across the entire 

thickness of the aquifer at the NWWF is estimated to be about 29.5 ft (9 m) (Renken and others, 2008). 

Aquifer Properties

Transmissivity values for the Biscayne aquifer derived from aquifer performance tests (APTs) 

vary widely throughout Miami-Dade County. Fish and Stewart (1991) reported a range in transmissivity 

values for the surficial aquifer, which includes the Biscayne and gray limestone aquifers, of about 49 to 

2.9 × 106 square feet per day (ft2/d); few of these estimates were made near the NWWF and WWF. 
7

Transmissivity values reported for the Biscayne aquifer, however, are generally higher than those 

reported for the gray limestone aquifer (Fish and Stewart, 1991; Reese and Cunningham, 2000). Results 

of an APT at the NWWF (Renken and others, 2008) yielded estimates of transmissivity ranging from 

1.67 × 106 ft2/d to 3.07 × 106 ft2/d. These data were acquired as part of a pathogen and chemical 

transport study to characterize and delineate flow through flow zones in the Biscayne aquifer near the 

NWWF. The flow zones constitute the main pathways of groundwater flow in the Lake Belt area (figs. 2 

and 5) and therefore, are an important component of the hydrogeologic framework in regards to both 

groundwater flow and contaminant transport (Cunningham and others, 2006a, b; Renken and others, 

2005, 2008; Shapiro and others, 2008). Delayed aquifer responses during these aquifer performance 

tests were attributed to the dual-porosity nature of the Biscayne aquifer. The flow was conceptualized as 

initial flow and storage depletion from high-porosity zones, and subsequent release of storage from less 

transmissive overlying and underlying units (Renken and others, 2005, 2008). These delayed aquifer 

responses are early-time phenomena and strongest near supply wells. 

An important aquifer property controlling contaminant transport rates through an aquifer is 

effective porosity. Although aquifer bulk porosity is defined as the percentage of the total aquifer 

volume composed of voids, or rock pores, effective porosity is usually smaller because it only includes 

those voids that are interconnected; therefore, bulk porosity is always an upper bound to any estimates 

of effective porosity. Cunningham and others (2009) describe the karstic Biscayne aquifer as a dualporosity system consisting of (1) matrix porosity (characterized by interparticle and separated vugs), 

which provides groundwater storage, and (2) megaporosity (characterized by macroporous touching 

vugs), which creates stratiform, areally extensive, groundwater-flow pathways, and less common 

bedding-plane and cavernous vugs and thin vertical solution pipes. Reported values of bulk helium 

porosity measured in the study area range from 5.5 to 50.2 percent (Cunningham and others, 2004a, 

2006b). Notably, a computer rendering of a megaporous outcrop limestone sample of the Biscayne 

aquifer from the Atlantic Coastal Ridge in Miami-Dade County has a megapore system dominated by 

burrow porosity calculated to be 78 percent (Cunningham and Sukop, 2011). 

Estimates of effective porosity are commonly made from the mean arrival time of a tracer 

breakthrough curve. A field tracer test conducted at the NWWF using rhodamine water-tracing dye and 

deuterium oxide yielded effective porosity estimates ranging from 4 percent to 41 percent. The 4 percent 

value considers all megaporosity high-flow zones, which is considered to be the lower range for 

effective porosity in this area, whereas the 41 percent value considers transport solely through what was 

hypothesized to be the main flow zone at an altitude of -31.5 ft (-9.6 m) (Renken and others, 2008). This 

approximation is based on the travel time of the peak breakthrough concentration, as the mean travel 

time could not be calculated for this test because it was prematurely terminated. Calibrated numerical 

models of the Biscayne aquifer in Miami-Dade County, which can also be used to estimate effective 

porosity, have incorporated values of 0.20 for effective porosity (Merritt, 1996, 1997; Langevin, 2001, 

Chin and others, 2010; Lohmann and others, 2012). 

Estimates of hydraulic storage properties are needed to simulate transient responses. The 

Biscayne aquifer is generally unconfined, thus, specific yield, which is a measure of the ratio of the 

volume of water an aquifer will yield by gravity drainage to the volume of the aquifer, is the storage 

property required. Specific yield is commonly assigned a value between 0.20 and 0.25 based on (1) the 

quantification of aquifer responses to known rainfall events and (2) estimates from calibrated models 

(Merritt, 1996, 1997; Langevin, 2001, Lohmann and others, 2012). 
8

Water-Budget Components

The hydrology of southern Florida is characterized by distinct wet and dry seasons, high 

evapotranspiration rates, and surface-water management features designed to control aquifer heads. In 

general, the regional groundwater flow direction in the Biscayne aquifer is from west to east, with 

eventual discharge into Biscayne Bay. Groundwater withdrawals from the Biscayne aquifer are 

considerable, and therefore substantially affect the groundwater flow field. The high transmissivities of 

the Biscayne aquifer facilitate rapid responses of groundwater levels to changes in rainfall and managed 

canal stages.

Ground-Water Withdrawals 

A principal hydrologic stress affecting flow through the groundwater system is well-field 

withdrawal. Groundwater is withdrawn from the NWWF and WWF for residential as well as 

commercial and industrial uses. The NWWF consists of 15 wells that have a maximum design-capacity 

pumping rate of 225 Mgal/d, making it the largest well field used for drinking water supply in Florida. 

The WWF consists of three wells that have a maximum design-capacity pumping rate of 25 Mgal/d. 

These design capacities are much larger than the actual and historical pumping rates. The actual 

volumetric daily pumping rates for the NWWF over the period of January 1996 to December 2004 

ranged from 18.5 to 166 Mgal/d (V. Walsh, Miami-Dade Water and Sewer Department, unpub. data, 

2008). Volumetric daily pumping rates for the WWF, between 1997 (when the well field became 

operational) and December 2004, ranged from 0 to 24.6 Mgal/d (V. Walsh, Miami-Dade Water and 

Sewer Department, unpub. data, 2008). 

There are many other supply well fields withdrawing groundwater from the Biscayne aquifer in 

Miami-Dade County (fig. 2). Total municipal groundwater withdrawals in Miami-Dade County for the 

1996-2004 period generally range from 350 to 440 Mgal/d (fig. 7). Groundwater withdrawals from the 

Biscayne aquifer for agricultural purposes are not accounted for in this investigation, but are estimated 

to be approximately 46.5 Mgal/d (Marella, 2009). Most agricultural groundwater withdrawals are made 

in the southern part of Miami-Dade County outside of the study area. 

Surface-Water/Groundwater Interaction 

The interconnected network of water-management features, including canals and control 

structures, is a unique aspect of the study area and southeastern Florida. These features have a strong 

influence on the hydrologic system. Water-management canals include primary canals, managed by the 

SFWMD, and secondary and tertiary canals maintained by the county or by local entities. Structures 

located along these canals are used to manage canal stages by controlling flow between canals. During 

the dry season, canal structure operations are designed to control water levels in the canals to maintain 

aquifer heads for water supply and to minimize saltwater intrusion near the coast. During the wet 

season, operations are designed mainly to provide flood control, whereby runoff from urban and 

agricultural areas collects in canals and is eventually discharged to Biscayne Bay. Flow within the canal 

network is generally from west to east, from primarily wetland areas through more urbanized parts of 

the county, with eventual discharge of surplus water to the ocean.

Groundwater in the Biscayne aquifer is closely coupled to the surface-water system because of 

the high transmissivity of the aquifer and its direct connection to the canals. Near the NWWF and WWF 

for example, the secondary Snapper Creek Extension Canal, (fig. 8), which extends from north to south 
9

directly east of the NWWF, is used to maintain heads in the aquifer during the dry season to prevent 

westward flow from urban contamination sources east of the canal. Thus, water levels in canals 

substantially influence the local aquifer heads and flow field. Canals can act as both a source of recharge 

to the aquifer (mostly in the dry season) and discharge (mostly in the wet season). In natural systems, 

baseflow refers to the portion of streamflow that is not runoff, that results from seepage of water from 

the aquifer into the stream. When groundwater levels fall below the stream stage, there can be losses of 

surface water into the aquifer. For purposes of this report, “net exchange” is defined as the gain or loss 

of water to the canal from the aquifer. Actual volumes or rates of recharge from canals and discharge to 

canals within a reach are difficult to quantify and are currently not well understood but are a large 

component of the water budget. Estimates of net exchange can be calculated using structure discharge 

data and are a surrogate for canal recharge and discharge to the aquifer. 

Rainfall, Evapotranspiration, and Runoff 

Aquifer recharge is typically the difference between rainfall and the sum of evapotranspiration 

and surface-water runoff. Rainfall is the largest source of water to the Biscayne aquifer in Miami-Dade 

County, with an estimated average of 53.5 inches per year (in/yr) during January 1, 1996, through 

December 31, 2004. This estimate is based on average annual precipitation data collected from rainfall 

stations throughout the county and downloaded from the SFWMD DBHYDRO environmental database. 

A substantial part of this rainfall, however, never directly recharges the aquifer because it is lost to 

(1) evaporation from ponded surface water (including lakes, canals, and wetlands), (2) plant 

transpiration, and (3) direct runoff to surface-water canals, some of which is eventually released into 

Biscayne Bay. 

Evapotranspiration, defined as the loss of water to the atmosphere from evaporation and plant 

transpiration, is considerable in southeastern Florida. The average evapotranspiration rate was 

40.4 in/yr during 1996 to 2004 for site 8 in German (2000), which is dry part of the year and not 

constantly wet or over open water (Lohmann and others, 2012). Based on this rate, a maximum of 

approximately 76 percent of the annual rainfall is lost through evapotranspiration. 

Runoff refers to excess water flowing into surface-water bodies. Runoff occurs during highprecipitation events and depends on land-use type and soil moisture content. Precipitation falling on 

impervious surfaces present in highly urbanized areas becomes runoff that is routed to nearby drains, 

ditches, streams, or canals. Some portion of runoff into surface-water canals may eventually recharge 

the Biscayne aquifer; while the rest is discharged to Biscayne Bay. Precipitation falls primarily in 

wetland, wooded, or grassy undeveloped areas near the NWWF and WWF, which are located in the 

western, less urbanized part of Miami-Dade County. Precipitation not lost to evapotranspiration or that 

does not recharge the aquifer is routed to the C-4 and C-6 Canals (Cooper and Lane, 1987). Estimates of 

runoff in Miami-Dade County are difficult to quantify and therefore not available. Chin and Patterson 

(2005) found that the ratio of runoff to rainfall in the C-103/C-16 basin was in relatively close 

agreement with the percentage of impervious area. This correlation is likely a result of the high 

permeability of soils and underlying aquifer materials and observed limited runoff from pervious areas. 

Runoff is not accounted for in this investigation because the assumption is that direct runoff is minimal 

in the study area. This assumption is consistent with the findings of Chin and Patterson (2005). 

Impervious (urban) areas compose about 11 percent of the study area surrounding the NWWF and 

WWF, with most of the impervious area located in the southeastern part of the study area, outside the 

area directly affected by the NWWF and WWF.
10

Simulation of Groundwater Flow near the Northwest and West Well 

Fields, 1996–2004 

A two-dimensional, transient MODFLOW-2005 model (Harbaugh, 2005) was developed to 

simulate groundwater flow in the Biscayne aquifer of Miami-Dade County. The model was calibrated to 

observed groundwater levels and estimated net exchange during 1996-2004, and was shown to 

reasonably represent groundwater flow patterns in the area. The model was adjusted to represent steadystate dry, average and wet conditions, and used to estimate travel-time capture zones and drawdown 

contours for the NWWF and WWF. Backward particle-tracking simulations run using MODPATH 5.0 

(Pollock, 1994), along with a lake residence-time post-processing adjustment, were used to develop 

travel-time capture zones; drawdown contours were calculated by subtracting water-level altitudes for 

simulations with and without groundwater withdrawals at the NWWF and WWF (and other major 

supply well fields nearby). Because of uncertainty in the hydraulic conductivity distribution, canal 

conductances, and representative values of effective porosity, the steady-state conditions were simulated 

using a stochastic approach, specifically, an unconstrained Monte Carlo analysis. 

Model Construction 

The model was constructed using MODFLOW-2005 (Harbaugh, 2005) to simulate the Biscayne 

aquifer flow system in Miami-Dade County for 1996-2004. The model accounts for daily canal-stage 

fluctuations; rainfall; evapotranspiration; water levels in Biscayne Bay, Card Sound, and Barnes Sound; 

and well-field withdrawals. The model was calibrated to observed historic Biscayne aquifer 

groundwater levels and net exchanges.

The model extends from Biscayne Bay, Card Sound, and Barnes Sound east and southeast of the 

study area to the L-31N and L-30 Canals in the western part of the study area (figs. 8 and 9). The 

northern model boundary is the C-9 Canal and the southern boundary is the C-111 Canal (fig. 9). These 

canals, as well as Biscayne Bay to the east and the L-30 and L-31N to the west, were chosen because 

they are hydraulic head boundaries and have daily observation data that can be used to specify the 

model boundary conditions. Stage and flow data for all of the canals making up the border of the active 

model domain can be used to specify best-estimates for head-dependent boundaries, or in the case of 

Biscayne Bay, time-varying specified-head conditions. While the Biscayne aquifer extends beyond 

these canals beneath wetlands to the west, the amount and variability of underflow is not known. Thus 

the area outside the northern, western, and southern head-dependent boundaries along the canals was not 

simulated and any exchange of water between these boundary canals and aquifer represents both 

underflow and canal exchange (fig. 9). 

Spatial and Temporal Discretization 

To simulate accurately groundwater heads, horizontal gradients, and flow near the NWWF and 

WWF, a uniform grid with a column and row width of 164 ft (50 m) was used throughout the model 

domain. This partitioning resulted in a model domain with 1,730 rows and 930 columns, yielding a total 

of 748,563 active model cells covering an approximate area of 1,550 mi2

. 

Multiple flow zones are known to exist in the Biscayne aquifer at some locations, but the aquifer 

was simulated using a single model layer for several reasons. Although the flow zones have been 

mapped over most of the study area, the connectivity of identified flow zones is unknown, they have not 
11

been mapped over the entire model domain, and there are only a few monitoring wells that measure 

water levels of specific flow zones. Furthermore, no APT data exist in the study area that provide 

vertical and horizontal hydraulic conductivity estimates of the individual flow zones, and no estimates 

exist of the spatial distribution of effective porosity. Simulation of the system using multiple layers to 

represent multiple flow zones would introduce substantial uncertainty into the analysis. A Monte Carlo 

analysis, which simulates the flow system thousands of times using thousands of possible realizations of 

the hydraulic properties, incorporates the uncertainty into simulated well-field capture zones and 

drawdown patterns. Furthermore, for the analysis to be tractable with available computational resources, 

a single realization had to run in less than a minute, and simplification of a one layer model helped 

control simulation run times. Simulations of groundwater flow in unconfined aquifers are often 

represented in two dimensions, instead of three, if the Dupuit assumptions hold (Anderson and 

Woessner, 1992). The Dupuit assumptions are that the hydraulic gradient is equal to the slope of the 

water table and that equipotentials are vertical. No aquifer can precisely meet these assumptions, but in 

most areas and for most times, the Dupuit assumptions are satisfied in the Biscayne aquifer. For these 

reasons, the model uses a single-layer to simulate two-dimensional groundwater flow in the Biscayne 

aquifer. 

Test simulations were conducted to evaluate the differences between simulating the aquifer 

using one or multiple layers (Appendix 1). An analysis was designed to identify any differences 

between the simulated hydraulic responses and travel times of a three-layer model with explicit 

parameters and a one-layer model with bulk hydraulic properties. This analysis demonstrated that 

estimates of travel-time capture zones and drawdowns from the one-layer model were slightly more 

spatially extensive than those same results from the three-layer model, but that the results were 

generally similar. Results from this investigation are consistent with Sepulveda and Kuniansky (2010) 

that demonstrated flow in the unconfined surficial aquifer system in Florida could be simulated 

reasonably well as a one-layer unit with calculated composite horizontal hydraulic conductivity based 

on a thickness-weighted, arithmetic mean of hydraulic conductivity values of multiple lithologic units.

The top of the model was set at land surface obtained from two different topographic datasets: a 

topographic contour land-surface map provided by ENP and a topographic map from NOAA (NGDC) 

Coastal Relief Model data (http://www.ngdc.noaa.gov/mgg/bathymetry/relief.html, June 2006). These 

datasets were combined and the centroid of each model grid cell was used to sample this coverage to 

determine the upper altitudes for each cell. The data and interpolations used for the land-surface 

altitudes are the same as those used in the Biscayne model developed by Lohmann and others (2012). 

Interpolated land-surface altitudes along the Biscayne Bay shoreline were set to a 0.0 ft altitude to 

reflect NAVD 88 (fig. 10). Bottom cell altitudes were assigned using a map of the base of the Biscayne 

aquifer (Fish and Stewart, 1991, p. 41; fig. 6). Fish and Stewart (1991) developed their map and cross 

sections of hydraulic conductivity using an integrated dataset of lithologic sample and core 

descriptions, APT results, and drilling flow rates throughout Miami-Dade County. 

Most temporal data used for model input or calibration were measured on a daily basis, with the 

exception of evapotranspiration rates, which were averaged from data collected at 15-minute intervals. 

One-day stress periods were used for the January 1996 – December 2004 period, yielding 3,288 

MODFLOW stress periods. 
12

Assignment of Aquifer Parameters 

The horizontal hydraulic conductivity field was created using APT data and spatial interpolation 

between test locations using ordinary kriging. A raw variogram, which mathematically describes spatial 

continuity, was calculated using an exponential variogram model with a correlation length of 16,404 ft 

(5,000 m) and a sill, or asymptotic value of variance, of 0.3 of the log-transformed hydraulic 

conductivity values. Hydraulic conductivities were calculated by dividing measured transmissivity by 

the full thickness of the Biscayne aquifer at that location. The locations of transmissivity estimates from 

APTs in Miami-Dade County (Fish and Stewart, 1991; Renken and others, 2008) and the interpolated 

horizontal hydraulic conductivity field, which will be referred to hereafter as the “base-case” hydraulic 

conductivity field, are shown in figure 11. This base-case field was utilized in the transient model as 

well as in the deterministic steady-state model simulations. These deterministic steady-state simulations 

are referred to as the “base-case” simulations. 

The Biscayne aquifer was simulated as a convertible MODFLOW layer, for which unconfined 

flow conditions are assumed when the head is below the top of the model cell and confined conditions 

are assumed when the head exceeds the top of a model cell. For this one-layer Biscayne aquifer model, a 

specific storage coefficient of 1.0 ft-1 was assigned for the confined case in MODFLOW, such that when 

water levels rise above the top of the aquifer (land surface) a unit head change results in a unit change in 

storage volume per unit area. This means that the volume above land surface is 100 percent water by 

volume. The storage term for unconfined aquifers is known as specific yield, which is a measure of the 

volume of water stored or released from an aquifer under gravity-drainage of the pore space per unit 

surface area of aquifer per unit decline in head. Specific yield was set to 0.20, based on previous 

calibrated models of this area by Langevin (2001) and Merritt (1996, 1997). 

For the deterministic base-case simulations, effective porosity was set to 0.12 based on the 

average of effective porosity estimates from tracer tests (4 percent) and effective porosity values used in 

calibrated models of the area (20 percent) (Langevin, 2001; Merritt, 1996, 1997, Lohmann and others, 

2012). This effective porosity value is within the range of bulk porosity estimates from Cunningham 

and others (2004a, 2006b) and estimates for the NWWF from Renken and others (2008). The value of 

0.12 is also the same as the mean value for effective porosity of groundwater flow zone layers used in 

the Lake Belt area by Shoemaker and others (2008a).

Boundary Conditions

Flows into and out of the model were simulated using specified-head, specified-flux, or headdependent-flux boundary conditions in MODFLOW. Flow in the unsaturated zone was not represented 

in this model because (1) the volume of water in the unsaturated zone is small relative to the total 

volume of water in the Biscayne aquifer, and (2) based on analysis of rainfall data and groundwater 

response time, the delay between infiltration of rainfall and groundwater recharge is thought to be less 

than a day. A no-flow boundary condition was assigned to the model base as there is a low-permeability 

confining unit that limits leakage between the two aquifers (Fish and Stewart, 1991). Additonally, any 

leakage between the Biscayne aquifer and the underlying portion of the surficial aquifer is considered to 

be small relative to other water budget sources and sinks, such as evapotranspiration and rainfall, and 

therefore could be neglected. 
13

Biscayne Bay

The cells representing Biscayne Bay, Barnes Sound, and Card Sound are adjacent to actively 

simulated onshore cells and are simulated as prescribed-head boundaries using the time-varying 

specified head or “Constant-Head (CHD) Package” in MODFLOW-2005. Water levels in these cells 

were assigned using daily tailwater stage data from control structure S-123. This structure, located along 

the coast, is centered along Biscayne Bay in the model and has been used previously as a proxy for 

water levels in Biscayne Bay (Langevin, 2001). Although Biscayne Bay is saline and density variations 

can affect groundwater flow patterns near the coast, variable density effects were not taken into account 

in the model. It is reasonable to expect that coastal density variations would not affect model results in 

the NWWF and WWF areas because they are far from the coast and separated from it by the salinity 

control structures and other hydrologic stresses included in the model. Model cells in Biscayne Bay not 

adjacent to active onshore cells were inactivated and set as no-flow cells.

Groundwater Withdrawals for Public Supply

All supply well fields in Miami-Dade County were simulated using the MODFLOW-2005 well 

(WEL) Package. Pumping rates were based on measured, daily volumetric pumping rates, where 

available. In some cases, only monthly volumetric totals were available; therefore, daily groundwater 

withdrawal rates were distributed over the number of days in the month and per well. When only hourly 

withdrawal rates were available, hourly rates were averaged for each 24-hour period. All well-field

withdrawal data were provided by the Miami-Dade County Water and Sewer Department. The dataset is 

the same as that used for the 1996 to 2004 time period of Lohmann and others (2012) and withdrawals 

are estimated in the same manner as discussed in that report. 

Canals 

Canals were simulated as head-dependent fluxes using the River (RIV) Package in MODFLOW2005. The RIV Package requires canal stage levels for each stress period, the river-bottom altitude and a 

conductance value for each canal cell. Daily headwater and tailwater stages for primary canal structures 

were obtained from the SFWMD DBHYDRO database and used to define river stages for the RIV 

Package. Stage data collected for important secondary canals were also included in the model. The most 

important secondary canal included in the model was the Snapper Creek Extension Canal (fig. 8), which 

is used to maintain high heads in the aquifer during dry periods to prevent east to west flow and has a 

substantial effect on the heads and drawdown in the NWWF area. If stage data were not available for 

included secondary and tertiary canals, stage data from primary canals hydraulically connected to these 

canals were used or modified during calibration. 

Other variables required by the RIV Package were taken directly, or calculated, from a spatial 

dataset provided by the SFWMD (Giddings and others, 2006). The canal conductance value (CRIV) 

required in the RIV Package was calculated using the following equation: 

M

WLK CRIV×× = (1)

where 

K is the horizontal hydraulic conductivity of the cell, in feet per day; 
14

L is the length of the canal reach in the cell, in feet; 

W is the width of the canal in the cell, in feet; and 

M is the flow length, in feet. 

The length of each canal attributed to each model cell was determined by superimposing a spatial 

representation of the canal on the model grid. The RIV Package was created by geospatially intersecting 

the GIS canal coverage with the model grid to calculate the length (L) of each canal-reach polyline in 

each model cell. The width (W) used in CRIV calculations for each RIV cell was calculated using the 

average of the top and bottom widths in the GIS coverage. Unlike other canal parameters, which can be 

directly measured or estimated, the flow length parameter (M) is not explicitly defined by the RIV 

Package. The parameter is similar in magnitude to the grid-cell size and, as Swain and others (1996) 

state, “in purely mathematical terms,” it is an estimate of the distance between the center of a canal (in 

cross section) and the point in the aquifer where the head in the aquifer is equal to the average head in 

the model cell. Therefore, the calculation of conductance in this investigation comes from the 

MODFLOW definition presented for equation 1 (McDonald and Harbaugh, 1988), but is slightly 

adjusted such that the flow length parameter is equivalent to the thickness of the river bed b ' in Swain 

and others (1996, eq. A1). Although this calculation is based on a vertical-flow relation, it is assumed 

that most water exchange between the canals and the surrounding aquifer occurs along the canal walls 

where sedimentation would be less than that along the canal bottom. Sedimentation was not explicitly 

incorporated in the calculation of CRIV but is assumed to be included into the flow-length parameter. 

The assumption of flow occurring along the walls of the canal implies that the channel side dimension 

would be more important to use in the calculation of CRIV than the canal width. It is assumed, however, 

that the conductance value—and seepage as well—are proportional to the width of the canal. Therefore, 

the estimated canal width (the average of the top and bottom width) is used for W in the CRIV 

calculation, and the adjustment of this value to actually account for wall seepage instead of canalbottom seepage is assumed to be accounted for in the adjustment of the flow length parameter M. An M

value of 32.8 ft (10 m) was used in the model, based on the 164-ft (50 m) cell dimensions, and results 

from the model were not shown to be very sensitive to this parameter. MODFLOW-2005 was modified 

to calculate canal conductance internally from the horizontal hydraulic conductivity (K) of the cell and 

the assigned conductance value of each RIV cell (Appendix 2), which was specified in the RIV Package 

and is equal to L WM × . Canal-bottom altitudes, which are required in the RIV Package, were also 

taken from an attribute in the Lower East Coast Sub-Regional MODFLOW (LECsR) model coverage.

 The western extent of the active model boundary is located at the L-30 and the L-31N Canals. A 

head-dependent flux is applied at this boundary, using the RIV package. The boundary flux is controlled 

by the specified canal stage and the conductance of the RIV cells. This boundary condition controls all 

flow of water coming into the model from the west. In effect, this boundary condition simulates a 

composite flux of water recharging the aquifer from the canals, and flowing beneath the canals from the 

west. This boundary condition is similar to that used for the western model boundary in the 

investigation of Lohmann and others (2012). The model in Lohmann and others (2012) was a multiplelayer model, and a head-dependent flux boundary condition was applied to the layers beneath the canal. 

The top layer of the model of Lohmann and others (2012) actively simulated the surface-water network

in this area, and exchange between the canal and the underlying aquifer. The approach used for the onelayer model documented in this report does not allow the flux from canal recharge to be distinguished 

from flux beneath the canal from the west. Furthermore, it does not allow a determination of the entire 
15

contributing area for each well field, because particles can only track backward into the canal, and 

cannot track into possible source areas in the inactive area west of the canal. Although some estimates 

exist for canal underflow (Sonenshein, 2001), the magnitude and variability of the flux of water from 

the west under the L-30 and L-31N canals is poorly understood. This boundary condition brings water 

into the actively simulated area to meet calibration targets, without quantifying canal underflow from 

areas west of the model edge. Testing of this boundary by actively simulating an area 

directly west of the NWWF indicated that water levels in the NWWF are not affected by adjusting 

the model boundary, thus extending the model area would not affect drawdown in the well field 

(Appendix 3).

Recharge 

Recharge in the model was simulated as a specified flux using the Recharge (RCH) Package. No 

attempt was made to reduce rainfall rates by estimated runoff rates. Daily measured rainfall rates at 48 

stations in Miami-Dade County were averaged over the 1996-2004 period to determine the spatial 

distribution of rainfall fluxes (in inches per day). The average daily rainfall from 1996 to 2004 at each 

of the 48 sites was plotted and contoured. The contour plot indicated that rainfall could be classified into 

(1) an area of high rainfall of 65.0 in/yr in the northern half of the model domain, and (2) an area of low 

rainfall of 48.8 in/yr in the southern half of the model. Because the NWWF and WWF lies in the 

northern portion of the model, daily rates from 11 rainfall stations located near the NWWF and WWF 

within the northern rainfall zone were averaged (66.5 in/yr) and applied uniformly throughout the entire 

model domain. 

Evapotranspiration 

Evapotranspiration was modeled as a head-dependent flux using the Evapotranspiration (EVT) 

Package. Physical parameters required by the EVT Package include an evapotranspiration surface 

(SURF); the extinction depth, which corresponds to the depth below SURF where evapotranspiration no 

longer occurs; and maximum evapotranspiration rates. Physically, the extinction depth generally 

corresponds to the root depth of vegetation, which is a function of land use and land coverage. SURF 

was defined to be land-surface altitude (fig. 10), and the extinction depth was defined to be 3.28 ft (1 m) 

below SURF. Langevin (2001) used extinction depths ranging from 0.49 ft (0.15 m) for barren land to 

2.3 ft (0.70 m) for upland forests. 

Daily average maximum evapotranspiration rates used in the EVT Package were calculated and 

processed from 15-minute time-series data collected at site 8 in German (2000) for initial use in the 

model of Lohmann and others (2012) for the same model domain. Open-water stations located in the 

Everglades and Water Conservation Areas (WCAs) 1, 2, and 3 were not used in this calculation of daily 

average maximum evapotranspiration rates, as these rates at open-water stations are generally higher 

than those at stations with vegetation. Site 8 is dry part of the year as opposed to other sites listed in 

German (2000), which are over open water or consistently wet. 

Quarry Lakes 

In the NWWF and WWF area, rock-mining activities have created large man-made quarry 

lakes that form as groundwater seeps into newly excavated areas. Near the NWWF, lakes are as close as 

650 ft from the nearest supply well. These quarry lakes are estimated to reach depths of 50 to 80 ft into 

the Biscayne aquifer, which in some cases, is as deep as or deeper than the level of the production zone 
16

of the public-supply wells, indicating the potential for direct hydraulic connection between the lakes and 

the production zone. 

Lakes were simulated using high hydraulic conductivity cells. Each cell whose centroid was 

within a lake was assigned a horizontal hydraulic conductivity of 328,084 ft/d (100,000 m/d). As 

another option, the MODFLOW Lake (LAK) Package could have been used to simulate the lakes. An 

analysis was undertaken to ensure that the high hydraulic conductivity approach was comparable to 

explicit representation of lakes using the LAK Package (Appendix 4). This analysis indicates that the 

use of high hydraulic conductivity zones to represent lakes produces similar results compared to those 

obtained using the LAK Package. 

 To assign lake cells, a 1999 aerial photograph was used to determine the outlines of the lakes 

during the timeframe of the transient model (fig. 12). These outlines were overlain on the model grid, 

and cells that had their centroid within a lake were assigned as lake cells. Lake depths were not 

estimated for the transient model, although lake depths are estimated and taken into account for the 

steady-state scenario models, as discussed later. Groundwater levels in the study area tend to be higher 

than a 0.0 ft altitude (Appendix 5, fig. 13) and lake water levels were assumed to be similar to or lower 

than surrounding aquifer heads. Therefore, for lake cells, the top altitude and SURF parameter in the 

EVT Package were assigned a value of 0.0 ft altitude. This approach ensured that water would be 

removed at the maximum evapotranspiration rate from lake cells. Lakes are 100 percent water by 

volume, thus the storage property and effective porosity were both set to 1.0 at all lake cells. 

The 1999 lake shapefile was used to define lakes in the transient model because lake images 

were collected near the middle of the 1996-2004 simulation period. The steady-state scenario model 

runs utilized a newer lake shapefile developed from 2004 aerial photography (fig. 12), in order to 

simulate more recent conditions.

Initial Conditions 

Initial heads everywhere were assigned a value of 3.28 ft (1 m). Some initial errors are evident 

when comparing the observed water levels at monitoring wells to simulated heads, but these errors tend 

to dissipate quickly (within days) and have no effect on the steady-state particle tracking and drawdown 

results, which are the results of interest in this investigation.

Model Fit

Model calibration targets for this investigation included daily groundwater levels at 25 

monitoring wells in the study area and monthly net exchange for two surface-water basins in the study 

area. The model was calibrated by manual adjustment of evapotranspiration extinction depth, canal 

conductance flow length, secondary and tertiary canal stage (those without stage data), and hydraulic 

conductivity to qualitatively match observed groundwater levels and net exchange and to represent 

generally observed features of the flow system. A qualitatively higher standard of model fit was 

required within the study area compared to the surrounding area. Thus, in areas outside of the influence 

of withdrawals at the NWWF and WWF, the model was only required to capture general hydrologic 

trends of the Biscayne aquifer. 

Model calibration is typically required to yield realistic simulation results from a deterministic 

model, and thus the model was calibrated for transient conditions. The boundary conditions were 

modified to represent steady-state conditions to estimate the well-field capture zones and drawdown. 

Calibration to transient conditions ensures that the calibrated hydraulic parameters and model 
17

construction reasonably represent the hydrologic system, whereas the steady-state simulations yielded 

long-term hydrologic responses to input stresses that occur during the period of interest. 

Comparison of Simulated and Observed Water Levels 

Daily simulated and observed water levels were evaluated in 25 monitoring wells (fig. 14). 

These monitoring wells were selected based on their location relative to the area of interest and the 

availability of daily data for the transient simulation period. The water-level comparisons (fig. 13 and 

Appendix 5) show that, compared to historical observations, simulated water levels near the NWWF and 

WWF have similar seasonal fluctuation patterns and, in many cases, similar maxima and 

minima. Locations of groundwater observation data used to calibrate model and mean error calculated 

from observed and simulated daily water levels. Area shown is Base map 2 (fig. 1). NWWF is 

Northwest Well Field. WWF is West Well Field.

Three statistics used to describe the model fit are mean error (ME), mean absolute error (MAE) 

and the coefficient of determination, R2. The ME, or average difference between daily observed and 

simulated values, is an overall indicator of how well a model simulates observed conditions at a specific 

location. It should not be considered as the only indicator, however, because a small ME value does not 

necessarily indicate a good model fit (Anderson and Woessner, 1992). The ME is a standard indication 

of any model bias present at a location. The MAE, a good statistical indicator of the average quantitative 

misfit at a location, is calculated by taking the average of the absolute values of the differences between 

daily observed and simulated water levels at each monitoring well. The R2, or coefficient of 

determination, is a statistical approximation of the relative variance of the misfit. R2 is calculated 

as follows: 

∑

∑

=

=

−

−

−= n

i

i

n

i

ii

obs obs

obs sim

hh

hh

R

1

2

1

2

2

)(

)(

1 (2)

where 

obs hi is an observed water level at location i; 

isim h is a simulated water level at location i; 

obs h is the mean observed water level; and 

n is the total number of observations. 

The closer R2 is to unity, the better the fit is assumed to be. These statistics were calculated 

for all observations at each monitoring site and averaged for all wells to obtain a general value for the 

study area. 

Observed groundwater levels for all wells for the simulated period ranged from -5.08 to 8.038 ft. 

The ME values for individual wells ranged from -0.653 to 2.139 ft, with an average of 0.009 ft (fig. 14 

and table 1), which indicates simulated values are generally lower than observed values. The MAE for 

individual wells ranged from 0.087 to 2.186 ft (table 1), with an average MAE of 0.52 ft. The average 
18

MAE of 0.52 ft is approximately 4 percent of the overall range of observed values (highest observed 

head minus lowest observed head, in feet) and approximately 9 percent of the average observed range at 

each of the 25 individual wells. R2 values ranged from 0.324 to 0.967, with an average value of 0.753. 

The largest ME and MAE generally occurred at monitoring wells located near supply wells (G-551, 

G-3253 and G-3259A) (fig. 14). It can be difficult to match simulated and observed heads near pumping 

wells because of transient well effects and model grid resolution. Furthermore, errors introduced in 

cases where volumetric well-field pumping totals or rates were converted to individual well rates using 

maximum pumping rates may have introduced additional error. A combination of these types of 

measurement error could have possibly resulted in overestimated or underestimated simulated heads at 

individual monitoring wells. R2 values at these three wells and at G-973 (just east of the study area) 

were also low (table 1, fig. 13, fig. 14, and Appendix 5). 

Another useful statistic for determining model fit is the standard deviation of observed minus 

simulated groundwater heads, normalized to the range in observed heads (Kuniansky and others, 2004). 

A target value of 0.1 or less is preferred in order to indicate a good fit to the observed head data. The 

standard deviation of the observed minus simulated heads for all 78,443 values was 0.758; when divided 

by 13.12 ft (the range in observed values), this yields 0.058, which is below the target of 0.1. 

Sensitivity analyses of the model indicated that simulated heads are relatively insensitive to 

changes in hydraulic conductivity because of the high transmissivity of the Biscayne aquifer. Moreover, 

sensitivity analyses indicated observed heads can be matched relatively well even with changes to 

various model input parameters. This finding indicates that history-matching to groundwater heads at 

monitoring wells in the Biscayne aquifer, by itself, does not necessarily ensure the model is an accurate 

predictor of actual groundwater flow patterns. 

Comparison of Simulated and Observed Net Exchange between Canals and the Aquifer 

Sensitivity analyses indicated that the hydrologic system is strongly controlled by canal stage. 

High hydraulic conductivities in the model allow observed heads to be matched with relative ease, so 

the model was calibrated also to estimated monthly net exchange. Estimated net exchange observations, 

or net recharge or discharge rates for the primary canals within a basin, were calculated using the C-2, 

C-3, C-4, C-5, and C-6 surface-water basins delineated by Cooper and Lane (1987; fig. 15). Net 

exchange for the C-2, C-3, C-4, and C-5 basins (labeled with asterisks in figure 15) were grouped into a 

single basin because flow measurements were not available between the inflow structure for the C-4 

basin at the western edge of the urban areas and the outflow structures for the C-2, C-3, C-4, and C-5 

basins. Thus, the composite net exchange for the C-2, C-3, C-4, and C-5 basins is henceforth called C-4 

net exchange. 

To calculate net exchange observations, the following hydrologic components were ignored: 

storage in the canals; direct rainfall onto canals and evaporation out of the canals; and internal flows in 

primary, secondary, and tertiary canals. Estimated net exchange observations were determined using 

structure inflows and outflows at basin boundaries (Langevin, 2001; Lohmann and others, 2012). 

Measured structure flows were acquired from the SFWMD DBHYDRO database. Calculated net 

exchange for a surface-water basin was assumed to be the difference between the sum of structure 

inflows and the sum of structure outflows for that basin. 

The model simulates the daily volumetric exchange of water between the aquifer and the canal 

cells. The daily net exchange rate of water into or out of the composite C-4 and C-6 canals was 

extracted from the model results. Comparisons between observed and simulated net exchange show that 

the model simulates seasonal fluctuations that are similar to observed values for both basins (fig. 16). A 
19

positive exchange represents a net loss of water from the canal and net gain in the aquifer for the basin 

(recharge), whereas a negative exchange represents a net loss of water from the aquifer and net gain in 

the canals for the basin (discharge). Monthly volumetric rates for estimated observed net exchange were 

used to evaluate model performance to reduce noise associated with canal storage effects, direct 

recharge and evaporation, and internal canal flows within the basin. ME, MAE, and R2 statistics were 

also calculated for the net exchange comparisons (fig. 16). In general, simulated net exchange appears to 

match the general seasonal trend reasonably well, although the model does tend to underestimate canal 

losses for the C6 basin and overestimate canal losses for the C4 basin. 

Comparison of Dry and Wet Season Groundwater-Level Altitudes in Miami-Dade County 

A comparison between potentiometric surfaces generated from observed data and potentiometric 

surfaces generated from simulated water levels shows that, over most of the study area, differences in 

water levels are within 1 ft. Average groundwater levels were calculated from the transient model 

results for May (dry season) and November (wet season), 1998. For purposes of comparison, published 

water-table altitudes (Sonenshein and Koszalka, 1996) for May and November 1993 were interpolated 

to grid-cell centers. Although the year 1993 does not fall within the transient model period, the year 

1998 was shown to be most similar to the year 1993 based on a comparison of average annual canalstage levels at a control structure on the L-30 Canal. Simulated daily water levels in the aquifer during 

May and November 1998 were averaged for each grid cell, and water levels during May and November 

1993 from Sonenshein and Koszalka (1996) were subtracted to calculate difference maps (figs. 17 

and 18). Water-level discrepancies are generally greater in the north, east, and far south of the model 

domain and near supply well fields than in the area of interest. 

Estimation of Well-Field Capture Zones and Drawdown 

Once the model was calibrated to groundwater heads and net exchange, steady-state simulations 

representing a range of hydrologic conditions were executed to estimate travel-time capture zones for 

10, 30, 100, and 210 days, and 0.25-ft and 0.1-ft drawdown contours for the well fields. In addition, two 

scenarios representing (1) seepage mitigation controls and (2) lake expansion were simulated, and the 

resulting capture zones and drawdown contours were compared to the results of the system without 

these proposed modifications. 

Methods

Travel-time capture zone and drawdown were simulated under steady-state conditions. A 

deterministic approach used the base-case distribution of hydraulic parameters, and the Monte Carlo 

simulations employed statistical variations of these parameters to investigate parameter uncertainty. To 

represent a range of hydrologic conditions, a dry, average and wet day were selected for the analysis. 

Maximum design-capacity pumping rates for all county supply well fields near and within the study 

area were used in order to estimate the largest extent of the capture zones and drawdown contours. 

Steady-state conditions were represented by converting boundary conditions (input using EVT, RIV, 

CHD packages for the transient simulation) to specified flux for the selected dry, average and wet days. 

To represent steady-state recharge for the average and wet hydrologic conditions, daily rainfall rates for 

the selected average and wet-condition days were annualized. For the dry conditions, the steady-state 

rainfall rate used was zero. 
20

Lake Representation—Bulk Hydraulic Conductivity Calculation 

In the transient model, it was assumed that quarry lakes excavated the full thickness of the 

Biscayne aquifer. Lake depths in the steady-state scenario models were adjusted to improve the 

accuracy of lake volumes and simulated travel times and drawdown. 

Quarry lake depths for the NWWF were compiled by Miami-Dade County Department of 

Regulatory and Economic Resources (RER) and utilized to create lake-bottom spatial coverages. Lake 

depths were based on bathymetric maps and rock-mining permits (J. Baker, Miami-Dade Department of 

Regulatory and Economic Resources, unpub. data, 2011). A permitted lake depth of 60 ft below land 

surface was assumed for lakes that lacked bathymetry information or permit records. For lakes near the 

WWF, available bathymetry was utilized for the two southernmost lakes (Solo-Gabriele, 2000). RER 

provided a bathymetric survey for the northernmost lake, and the average depth of this lake was used for 

the smallest lake directly south of it. Although the lake spatial coverage used in the steady-state 

simulations was generated in 2004, bathymetric surveys were not necessarily collected in the same year. 

The available lake-bottom data were used with the general assumption that the lake depths in 2004 were 

the reported surveyed depths. 

Lake-bottom altitudes were calculated for each grid cell by intersecting the grid cell centroids 

with interpolated raster datasets developed from the bathymetric points. Where aerial lake data indicated 

expansion of a lake after 2004, only bathymetric data within the 2004 lake polygon were utilized. On a 

cell-by-cell basis, lake depths were subtracted from the estimated thickness of the Biscayne aquifer to 

yield the thickness of the Biscayne aquifer in each model cell beneath the lakes. 

To represent the lakes accurately in a one-layer model, a bulk hydraulic conductivity value was 

calculated for lake cells. The calculation used lake depth, an assumed 328,084-ft/d (100,000-m/d) lake 

hydraulic conductivity, the underlying aquifer thickness and the Biscayne aquifer hydraulic conductivity 

obtained from the corresponding transient model cell in the base-case simulation (fig. 11). The bulk 

hydraulic conductivity was calculated using: 

aquifer aquifer lake lake

bulk

total

K b Kb

K

b

× +× =

(3) 

where 

Kbulk is the bulk hydraulic conductivity of the cell used in the hydrologic scenario models, 

in feet per day; 

Kaquifer is the hydraulic conductivity of the aquifer beneath the lake, in feet per day; 

baquifer is the thickness of the aquifer beneath the lake, in feet; 

Klake is 328,084 ft/d (100,000 m/d); 

blake is the thickness of the lake in feet; and 

btotal is the total sum of baquifer and blake. 
21

Determination of Drawdown and Travel-Time Capture Zones 

Drawdown was calculated as the difference between heads from the steady-state simulation 

without pumping and heads from the steady-state simulation with maximum design-capacity pumping. 

Drawdown contours of 0.25 ft and 0.1 ft were evaluated.

Groundwater flow velocities were calculated by MODPATH 5.0 using cell-by-cell Darcy fluxes 

from MODFLOW-2005 divided by the specified effective porosity values used with MODPATH for 

each model cell. Particles were tracked backwards using MODPATH 5.0 until they discharged to 

boundaries or internal sources or sinks. The weak sink option was used in MODPATH, which allows 

particles to pass through cells with an internal boundary condition, such as a cell that has a well 

(specified flux) or canal (head-dependent flux). This allows some particles to pass beneath canals 

(representing underflow) or move past a low discharge well, rather than terminate at the cell with the 

applied boundary condition. Along the L-30 and L-31N Canals, which are at the western model 

boundary, particles terminate at the canal because the cells to the west are inactive. The canals provide 

a strong head control, thus the pathlines backtracking from the well were tested and shown to be 

unaffected even if the aquifer to the west is set to be actively simulated. Particles that terminate at this 

boundary, however, may have been able to move past this boundary, indicating that they may have 

moved from the natural wetlands towards the well field, if this area were actively simulated 

(Appendix 3). 

Particles were placed around individual public-supply wells in the NWWF and WWF areas in all 

steady-state simulations. A total of 600 particles were placed along the cell faces containing supply 

wells in the NWWF and a total of 120 particles were placed along the cell faces containing supply wells 

in the WWF. Particles were then tracked backwards in time until they reached internal sources or 

boundaries. Forward particle tracking may yield different results than backward tracking, particularly in 

complex three-dimensional simulations (Barlow, 1997). Groundwater flow in the Biscayne aquifer is 

predominantly horizontal and the model used in this study is two-dimensional. A forward-tracking 

analysis was used for the one-layer model with canal cells in the supplemental simulations 

discussed in Appendix 1. Capture-zone results from this forward tracking simulation were not markedly 

different than results from an analogous backward-tracking simulation. Additionally, with the weak sink 

option, particles did travel from areas beyond the test canals through the weaker canals (these were 

internal to the boundary of the actively simulated area) moving towards the well. 

For particles that travel through lakes, particle endpoints were calculated using a post-processing 

technique that applied residence-time adjustments (Chin and others, 2010) (Appendix 6) to account for 

transport in open water. Residence times of lakes were calculated by dividing the lake volume by the 

simulated steady-state lake outflow into the surrounding aquifer. Lake volumes were calculated using 

lake-bottom altitude data provided by RER. Lake volumes were unchanged between all simulations, 

whereas lake outflows were different in each simulation and calculated using ZONEBUDGET 

(Harbaugh, 1990). The cell-by-cell flow file from the maximum design-capacity pumping steady-state 

run was adjusted to calculate flows from just the lakes (as opposed to total flows between cells in the 

one-layer model) by multiplying lake-cell flows by the ratio of the product of the lake thickness (blake) 

and the lake hydraulic conductivity (Klake) to the total Biscayne aquifer transmissivity. This adjusted 

cell-by-cell flow file was used as the input into ZONEBUDGET. The particle locations throughout time 

in the steady-state models were assumed to follow the actual particle pathlines through the lakes, 

although the lakes actually have more complicated mixing dynamics. The use of a residence time is one 

way to implement lake mixing without explicitly simulating the mixing process because a residencetime approach assumes the lake is well mixed. 
22

Residence times were calculated and used to adjust the particle locations by determining where 

adjusted particles would be located at 10, 30, 100, and 210 days. Once a particle was tracked backwards 

to a lake, the particle location was adjusted to either the inflow side of the lake, if the residence time was 

such that the particle was located inside the lake at 10, 30, 100, or 210 days, or allowed to travel past the 

inflow of the lake and into the aquifer beyond, if the residence time was such that the particle was no 

longer inside the lake during the 10-, 30-, 100-, or 210-day period being tracked. For example, if a 

particle reached a lake within 10 days and this lake had a residence time of 90 days, the particle would 

need to be tracked backwards for an additional 110 days from the inflow side of the lake to obtain the 

210-day capture-zone location. For the 30-day capture zone, however, the particle location would be the 

inflow point of the lake (which corresponds to where the particle would ultimately exit the lake when 

tracked backwards) because the residence time (90 days) and tracking time (10 days) are collectively 

larger than the required capture-zone endpoint (30 days). This residence-time approach is conceptually 

more correct than assuming porous media flow and a porosity of 1.0, because the latter does not 

account for mixing. Further information and justification on the particle post-processing method is 

discussed in Appendix 6. Finally, ending particle locations for specified times were plotted on a map to 

yield a spatial estimate of capture zones for the individual well fields. 

Steady-state deterministic simulations under different hydrologic conditions were conducted 

using the base-case horizontal hydraulic conductivity field, with bulk hydraulic conductivities for lake 

cells, and a uniform value of 0.12 for the effective porosity of the Biscayne aquifer and 1.0 for lakes. 

These base-case deterministic simulations were compared with Monte Carlo simulations to evaluate the 

benefits of considering parameter uncertainty using the Monte Carlo analysis compared to the 

traditional deterministic approach.

Evaluation of Uncertainty Using Monte Carlo Analysis 

Hydraulic conductivity and effective porosity are two physical parameters that are important in 

controlling groundwater flow and transport velocity. Although these parameters are estimated at APT 

locations, a large degree of uncertainty typically exists for parameter values between APT locations, 

while a smaller measurement-error uncertainty is present at those locations. A Monte Carlo analysis was 

performed to estimate capture zones and drawdown while considering uncertainties in hydraulic 

conductivity and effective porosity. A typical Monte Carlo analysis evaluates the effect of parameter 

uncertainty on model results by simulating a large number (hundreds or thousands) of realizations for 

which parameters are varied according to an assumed distribution. For this investigation, horizontal 

hydraulic conductivity and effective porosity were both varied. Canal bed conductances are also 

assigned the value of the aquifer hydraulic conductivity of the cell in which the RIV boundary condition 

is applied, so they vary as well Preliminary analyses indicated that 10,000 stochastic runs would be 

required to accommodate sample parameter uncertainty. Therefore, 10,000 conditioned hydraulicconductivity fields were created, 10,000 RIV Packages were created based on these conductivity fields, 

and 10,000 effective porosity values were sampled from a distribution representative of field and 

literature values. Each of the 10,000 realizations was simulated for each selected hydrologic condition. 

Accounting for Uncertainty in Spatial Distribution of Hydraulic Conductivity

Aquifer performance test data (Fish and Stewart, 1991; Renken and others, 2008), which were 

used to create the base-case hydraulic conductivity field (fig. 11), were also used to condition the 10,000 

hydraulic conductivity realizations. Horizontal hydraulic conductivity fields were generated using the 
23

sequential Gaussian simulation (SGSIM) method in PEST’s FIELDGEN utility, which creates 

stochastic fields using a variogram model and conditioning data points (Deutsch and Journel, 1998; 

Doherty, 2010). Specifically, random variability multiplier arrays were created using the SGSIM 

method. The hydraulic conductivity field realization used in each model simulation was then 

calculated by multiplying the base-case hydraulic conductivity field by the stochastically generated 

multiplier array: 

jihjirjih ),(),(),( r BC = KmK . (4) 

where for each realization r, 

i is the model row number, 

j is the model column number, 

jih ),( r K is the horizontal hydraulic conductivity for cell ( ) i j , , 

m jir ),( is the multiplier for cell ( ) i j , , and 

h ij BC , ( ) K is the base-case horizontal hydraulic conductivity field interpolated from Fish and 

Stewart (1991) and Renken and others (2008) for cell ( ) i j , . 

The variogram used in the SGSIM controls the variability and correlation structure of the 

stochastic fields. For this study, an exponential variogram with a correlation length of 13,123 ft 

(4,000 m) was used to generate the multiplier arrays. This correlation length was selected by 

constructing a variogram through curve-fitting based on the percentage of Biscayne aquifer flow zone 

thicknesses to the total thickness of the aquifer as defined in Cunningham and others (2006a, b). At APT 

locations (conditioning data), the multiplier array was assigned an expected value of 1.0. This 

assignment results in a mean hydraulic conductivity at APT locations that equals the hydraulic 

conductivity derived from APTs. The SGSIM variogram included a log-based nugget of 0.01 to 

represent APT measurement error. For example, for an APT-measured hydraulic conductivity value of 

6,484 ft/d (1,976 m/d), the log-based nugget of 0.01 results in a 3-standard deviation hydraulic 

conductivity range between approximately 3,250 and 12,937 ft/d. The sill represents the total variance 

of the multiplier array for distances longer than the variogram range. A log-based sill of 0.1 was used 

that allows hydraulic conductivity to vary within the range of approximately 730 and 57,611 ft/d 

(3 standard deviations) about the mean of 6,484 ft/d. The maximum and minimum values for each of the 

10,000 realizations of hydraulic conductivity were sorted and the 2.5 percentile, median, and 

97.5 percentile of the maximums and minimums were determined to better understand the range of 

hydraulic conductivity (table 2). The maximum and minimum of each realization are extreme values at 

one individual grid cell and not representative of the majority of grid cell values for hydraulic 

conductivity, which will tend to be close to the mean value. Three realizations (fig. 19) are examples 

of the products of the multiplier array and the base-case hydraulic conductivity field. Lake 

bulk hydraulic conductivity has not been replaced where lakes are present in these three examples.
24

 All three realizations are unique, yet similar to each other in that areas of higher conductivity

 and lower conductivity are common to each of the three, because the conditioning 

 values constrain the estimated conductivity at APT locations. 

For cells with lakes, the hydraulic conductivity values from the base-case or 10,000 Monte Carlo 

realizations were not used directly. Instead, bulk hydraulic conductivities were calculated using the lake 

depth, lake hydraulic conductivity, and the aquifer hydraulic conductivity of the cell the lake is in for 

each realization. This was done in order to appropriately represent the high-conductivity lakes and their 

depths in each realization used in the model. 

Accounting for Uncertainty in Effective Porosity

The uncertainty associated with effective porosity was evaluated concurrently with uncertainty 

in horizontal hydraulic conductivity. Unlike hydraulic conductivity realizations that contained 

variability throughout the model domain, a uniform effective porosity field was used in each realization. 

Although effective porosity is expected to vary spatially, and correlation between effective porosity and 

hydraulic conductivity in flow zones may exist, few data exist to determine the spatial variability of 

effective porosity. Use of a non-uniform effective porosity would have introduced an additional element 

of uncertainty. Based on a range of effective porosity values reported in the literature for the Biscayne 

aquifer—4 percent (Renken and others, 2008) to 20 percent (Langevin, 2001; Merritt, 1996, 1997)—a 

representative intermediate effective porosity value of 12 percent was chosen. For the base-case 

simulations, this porosity value was used, except in the lake cells where the effective porosity was set 

to 1.0. 

Effective porosity values were sampled from a lognormal distribution created using the PEST 

utility RANDPAR (Doherty, 2010). The seed for each of the 10,000 instances of RANDPAR was 

determined using a random number generator. The log normal effective porosity distribution (fig. 20) 

had a mean of 0.12 (log-transformed -0.92082) and a log-transformed standard deviation of 0.16. Most 

effective porosity values used in the 10,000 stochastic realizations centered around the mean value of 

0.12, with values from 0.04 to 0.36 falling within 3 standard deviations about the mean. Values used for 

all realizations ranged from 0.04 to 0.50 (fig. 20).

Selection of Stress Period for Steady-State Particle-Tracking Simulation 

A variety of different hydrologic conditions were evaluated to account for all reasonable 

groundwater flow fields near the well fields. Dry conditions, defined as a period of no recharge, low 

aquifer heads and low canal stages, are expected to result in capture zones with the largest areal extent. 

Wet and average hydrologic conditions were also simulated to account for potential differences in flow 

directions and travel times. Stage levels for tailwater at the S-337 and G-119 control structures, which 

control surface-water levels near the NWWF and WWF, were used to guide the selection of 

representative daily hydrologic conditions during the 9-year transient model period (fig. 21). Control 

structure S-337 (tailwater) is located on the L-30 Canal (fig. 8) northwest of the NWWF, and control 

structure G-119 (tailwater) is located on the Tamiami Canal north of the WWF. Stage data from the 

G-119 structure were used to select wet, average, and dry days that would be used to create steady-state 

models to simulate these different conditions, and data from the S-337 structure were used to verify the 

selected days. 

Dry hydrologic conditions were evaluated using June 19, 2004, when the minimum stage value 

for the 9-year transient period was recorded at G-119 and a low stage was recorded at S-337 (fig. 21). 
25

Average hydrologic conditions were determined by selecting a day, December 19, 2002, when observed 

tailwater stages were close to the average calculated tailwater stage at G-119. Wet hydrologic conditions 

were evaluated using October 16, 1999, when the maximum stage value for the transient period was 

recorded at both structures. 

Simulated Time of Travel and Drawdown for Dry, Average, and Wet Hydrologic Conditions

Representative dry, average, and wet hydrologic steady-state conditions are simulated and 

results are compared herein to determine the largest combined area enclosed by capture zones and 

drawdown contours. Analyses included both deterministic and Monte Carlo simulations. 

Base-Case Deterministic Simulations 

For each of the three hydrologic conditions, a deterministic simulation was run. Simulated 

drawdowns were contained within the surrounding canals and comparable to one another for all three 

runs. Simulated backward-tracked particle positions for the base-case simulations indicate that the 

capture zones for the NWWF extend further to the west, and are generally more uniform and radial 

toward the west than toward the east (figs. 22, 23, and 24). Ambient flow direction in the aquifer is 

generally west to east, which increases the length of the path that a particle travels in a set amount of 

time, making the capture zone extent larger toward the west than the east. When a simulated particle 

reaches a lake, the residence-time adjustment effectively increases or decreases the areal extent of 

capture-zone contours when compared to particles that do not come into contact with lakes, depending 

upon the calculated residence time. This results in a generally more irregular shape of the capture zone 

east of the well field than west of the well field. In addition, the Snapper Creek Extension Canal, located 

east of the NWWF (fig. 8), maintains high water levels in the aquifer and acts as a strong sink for 

backward-tracked particles. The dry-conditions simulation (fig. 22) appears to yield more spatially 

extensive time-based capture zones when compared to average (fig. 23) and wet (fig. 24) conditions. 

This could be attributed to the strong influence on recharge by the canals during dry conditions. 

Particles reach the closest lake to the north of the NWWF within 10 days. 

For the WWF, travel-time results for the three conditions are generally similar, with particles 

reaching a lake closest to the WWF within 100 days (figs. 22, 23, and 24). Capture zones for wet 

hydrologic conditions extend slightly farther to the north of the WWF than capture zones for dry and 

average hydrologic conditions (fig. 24), whereas capture zones for dry and average hydrologic 

conditions were generally more radially symmetric around the well field. This indicates that the capture 

zones for the WWF are more sensitive to recharge rates and fluctuations of canal stages than those for 

the NWWF. 

In general, the 0.25-ft and 0.1-ft drawdown contours differ minimally between simulations for 

the three hydrologic conditions (figs. 22, 23, and 24). 

Monte Carlo Simulations 

The results from the base-case simulations represent deterministic outcomes for the capture 

zones and drawdown contours of the NWWF and WWF under different hydrologic conditions. These 

results only represent one possible realization of hydraulic conductivity, canal conductance, and 

effective porosity, and do not account for the effects of parameter uncertainty on model results. To 

address parameter uncertainty in model results, a set of 10,000 possible physical realizations using a 
26

reasonable range of hydraulic conductivity and effective porosity values was simulated for each of three 

hydrologic conditions. 

Particle positions after 10, 30, 100, and 210 days of travel for the 10,000 dry hydrologic 

condition realizations generally cover a larger extent than those simulated for the base-case, dry

conditions (figs. 22 and 25). The 210-day particle positions from the Monte Carlo simulations cover a 

larger areal extent than the deterministic base-case simulation particle endpoints because many of the 

10,000 realizations have hydraulic conductivity and effective porosity values that differ markedly from 

base-case values. Many 100-day and 210-day capture-zone particles from the Monte Carlo simulations 

reach the western model boundary (fig. 25), whereas the 100-day and 210-day capture-zone particles for 

the base-case deterministic run do not (fig. 22). The Snapper Creek Extension Canal east of the NWWF 

functions as a strong sink for the backward-tracked particles limiting their movement to the east of the 

NWWF in both the base-case and Monte Carlo simulations. The Snapper Creek Extension Canal is 

reached by backward-tracked particles within 100 days (fig. 25). Several lakes surrounding the NWWF 

were reached by backward-tracked particles within 10 days. Particles travel to lakes closest to the WWF 

within 10 and 30 days and many particles tracked backwards to the western boundary condition within 

100 days. 

A statistical technique developed by Varljen and Shafer (1991) was utilized to statistically 

evaluate results from the 10,000 individual simulations for each hydrologic condition. The approach 

uses particle endpoints to delineate median capture-zone contour lines and the 95-percent CI zone of 

uncertainty. Particle endpoints are grouped into radial classes based on a specified class size (fig. 26). 

The class size must be small to satisfy the assumption that radial differences between particles are small, 

but also large enough to produce smooth contour lines. 

Relative to the spatial center of withdrawals for each well field, particles were sequentially 

ordered in the counter-clockwise direction starting at the 0-degree baseline (fig. 26). Particles were then 

separated into particle classes using 12,000 and 2,400 particles per class for the NWWF and WWF, 

respectively. This partitioning resulted in a total of 500 radial-slice particle classes for each well field. 

Within each class, particle distances from the centroid of withdrawals were used to statistically estimate 

travel-time capture zones. Varljen and Shafer’s (1991) approach assumes that the particles are normally 

distributed in each direction; this assumption allows determination of the 95-percent CI zone of 

uncertainty using the calculated mean and standard deviation calculated for each particle class. 

However, because stochastic samples of effective porosity values are log-normally distributed and lakes 

noticeably change the areal extent of particles depending upon their residence time, the calculation of 

the 95-percent CI was based on the 2.5 and 97.5 percentiles of the cumulative distribution function of 

the particle endpoints in each class. The 50th percentile in a given particle class is equivalent to the 

median distance for that travel time from the center of the well field. The 95-percent CI about this 

median was approximated by discarding those endpoints whose distances were less than the 

2.5 percentile or greater than the 97.5 percentile within each radial class. 

Simulated particle positions for dry conditions surrounding the NWWF for the 210-day capture 

zone for 100 of the 10,000 Monte Carlo simulations show the dense distribution of particles along the 

inflow boundaries of the lakes because of the residence-time adjustment of the particles tracked 

backwards through lakes (fig. 26). West of the NWWF supply wells, a more normal particle density 

distribution is evident. 

For the steady-state dry, average, and wet hydrologic conditions, the 95-percent confidence 

interval, the median, and the base-case deterministic simulated particle endpoints for 210, 100, 30, and 

10 days are presented. In the absence of lakes, the width of the 95-percent CI band reflects ambient flow 

conditions and uncertainty in the distribution of hydraulic properties and fluxes. In the presence of 
27

lakes, the width of the 95-percent CI band also reflects the uncertainty of residence times of the lakes 

and particle adjustments in post processing. The extent of the 95-percent CI contours differ depending 

on the hydrologic conditions. 

Simulated median capture zones and 95-percent CI for dry hydrologic conditions were 

determined for 210, 100, 30, and 10 days (figs. 27, 28, 29, and 30). The 210-day median capture zone 

for the Monte Carlo simulations are similar to but cover a larger spatial extent than the capture-zone 

particle positions for the base-case simulation, showing a greater extent toward the west than toward the 

east, (fig. 27). The width of the 95-percent CI is generally narrower north, east and southeast of the well 

field, where lake density is high, than west of the well field, because of (1) the high values of bulk 

hydraulic conductivities of lake cells, (2) an effective porosity of 1.0 for lakes in each realization, and 

(3) residence-time adjustments of the backward-tracked particles. In the area with lakes around the 

NWWF, the 95-percent CI is irregular for the 210-day and 100-day capture zones because individual 

particles may traverse more than one lake depending upon the lake residence times. The width of the 95-

percent CI is approximately 3 mi west of the NWWF away from the lakes, indicating greater uncertainty 

in particle transport rates. As with the deterministic results, results from the Monte Carlo simulations 

indicate that particles reach lakes closest to the NWWF within 10 days (fig. 30). 

Simulated median capture zones and 95-percent CIs for the average hydrologic conditions 

(figs. 31, 32, 33 and 34) are similar to those for the dry hydrologic conditions (figs. 27, 28, 29 and 30) 

although less spatially extensive. Simulated median capture zones and 95-percent CIs for wet 

hydrologic conditions (figs. 35, 36, 37, 38) are also similar to those for average conditions for the 

NWWF, but are substantially different in the cases of 210 days and 100 days for the WWF (figs. 31, 32, 

35, and 36). The 95-percent CI bands for the WWF are more radially symmetric (except in the presence 

of lakes) for dry and average hydrologic conditions than for wet conditions, and extend farther to the 

north of the WWF for wet conditions than for dry and average conditions. For the NWWF, the outer 

95-percent CI contour for the 210-day capture zone for dry hydrologic conditions almost fully 

encompasses the outer 95-percent CI contour for the 210-day capture zone for wet hydrologic 

conditions (fig. 39). For the WWF, the outer 95-percent CI contour for the 210-day capture zone for dry 

conditions does not fully encompass the outer 95-percent CI contour for the 210-day capture zone for 

wet conditions, which extends farther to the north (fig. 39). This finding indicates that the WWF has 

different recharge source areas depending on hydrologic conditions. Therefore, consideration of dry, 

average and wet hydrologic conditions is important when defining the maximum areal extent of capture 

zones for the WWF. 

The 30-day capture zone for the WWF extends further toward the southwest for dry conditions 

than for average or wet conditions (figs. 29, 33, 37). The upper bound and median of the 95-percent CIs 

of the 30-day capture zone intersect a lake to the southwest of the WWF in this case, indicating faster

transport of an introduced-chemical constituent from this lake to the WWF during the dry season than 

during the wet season. 

The longer-travel-time capture zones, 210-day and 100-day, indicate the area surrounding the 

well fields that ultimately contribute groundwater to, and within which introduced-chemical constituents 

may reach well fields, whereas a 10-day travel-time capture zone, indicates an area within which 

introduced-chemical constituents may reach well fields in a short amount of time. Because of the 

residence-time approach to particle transport through lakes, which is based on residence-time mixing, 

any lake that intersects a capture zone can contribute to the well field within that capture-zone travel 

time. The nearest lake north of the NWWF is approximately 650 ft from the closest supply well. The 

median and upper bound of the 95-percent CIs of the 10-day capture zone intersect this lake, for dry, 

average, and wet conditions. The results indicate that, because of the proximity of the lake to the 
28

NWWF, an introduced-chemical constituent could reach the NWWF within 10 days. Mixing processes 

in lakes and interactions between the lake and the surrounding aquifer introduce substantial uncertainty 

as to whether an introduced-chemical constituent in a lake will reach a well field within 10 days if the 

lake is within a 10-day capture zone. The possibility exists, nonetheless, and is worthy of consideration 

in management of wells that are close to lakes.

A composite 210-day capture zone for each well field was estimated by delineating the entire 

area in which particles were captured by each well field for both the dry and wet conditions. The upper 

bounds of the 95-percent CIs were overlaid on one another and the union of these two contour lines was 

formed (fig. 40). The composite capture zone shows the largest possible estimate (within the 95-percent

CI zone of uncertainty) for 210 days of travel for both well fields based on this analysis. The composite

210-day capture zone differs substantially in shape and extent of the 210-day capture zones used by the 

county to establish Wellhead Protection Areas (figs. 3 and 40). There are substantial differences in the 

models used to generate the capture zones: (1) the model described herein accounts for the presence and 

transport of particles through lakes; (2) the Monte Carlo analysis used to determine a 95-percent CI

accounts for uncertainty in hydraulic property distribution; and (3) the pumping rate for the WWF has 

been updated from 140 Mgal/d to 25 Mgal/d herein. 

Drawdowns of 0.1 ft and 0.25 ft for the NWWF and WWF were also estimated using the Monte 

Carlo simulations. The median and 95-percent CI of the drawdown in each cell were first calculated for 

each of the hydrologic conditions (dry, average, and wet conditions). Then the median and the upper 

and lower bounds of the 95-percent CI were contoured for the 0.1-ft and 0.25-ft drawdowns for each of 

the hydrologic conditions (figs. 41, 42, and 43). There are minimal differences in the median and 

95-percent CIs of the drawdowns for the different hydrologic conditions, likely because the the canals 

provide an unlimited source of water and the hydraulic conductivity is high, so pumping simply draws 

more water into the aquifer from the canals, maintaining the groundwater levels. Drawdown contours 

around the WWF coalesce with drawdown contours around other major well fields. For purposes of this 

analysis, a maximum design-capacity pumping rate was also used for well fields outside of the NWWF 

and WWF. Deterministic base-case simulated drawdowns (figs. 22, 23, and 24) coincide with the 

median drawdowns obtained from the Monte Carlo simulations and, therefore, are not shown. 

Canals, both within the model domain and at the model boundaries, act as sources of water 

and particle discharge locations, as can be seen in drawdown and capture zone results. Except in a few 

cases, particles did not pass through internal canals, even when the weak sink option was used in 

MODPATH 5.0, which indicates internal canals act mainly as strong sources of water. Particles did, in 

some instances, pass through secondary internal canals when they were tracked to their final discharge 

locations beyond the 210-day extent shown. 

 
29

Two Proposed Scenarios: Potential Effects of Seepage Mitigation and 

Lake Expansion on Capture Zones and Drawdown 

Two Monte Carlo simulations of three proposed anthropogenic changes to the natural 

environment were performed to evaluate their potential effects on well-field capture zones and 

drawdown. The first scenario investigates the effects of a recharge canal near the NWWF and restriction 

of flow from the western boundary west of the WWF. The second scenario evaluates the effects of 

proposed lake expansion on the capture zones and drawdown of both well fields. 

Proposed Scenario 1: Investigation of Effects of Seepage Mitigation

Construction of a 7-mi-long slurry wall located west of the WWF has been proposed as a means 

to reduce the seepage of groundwater from ENP beneath the L-31N Canal. The proposed slurry wall 

would be constructed west of and parallel to the L-31N Canal between the Tamiami Canal and SFWMD 

primary canal structure G-211 (Lake Belt Mitigation Committee, 2009). A 1,000 ft section of this wall 

has been constructed to test the feasibility of this project. At the time of this writing, plans are in place 

to finish the 7-mi slurry wall section and possibly complete the wall to a greater subsurface depth than 

the test slurry wall. The existing test slurry wall section is 28 inches (in.) wide and 18 ft deep and does 

not extend to the full depth of the Biscayne aquifer. The location of the proposed slurry wall between 

the L-31N Levee and the L-31N Canal is designed to reduce lateral seepage of groundwater beneath the 

canal from the Everglades toward the eastern part of the county, but not to inhibit the seepage from the 

canal into the aquifer. The L-31N Canal and undeveloped and protected areas west of the canal could 

potentially serve as sources of water east towards the WWF. 

The model cannot explicitly evaluate the effects of a slurry wall west of the L-31N Canal, as the 

canal coincides with the western model boundary. The head-dependent flux boundary condition 

imposed at the L-31N Canal represents composite sources, leakage from the canal into the aquifer and 

underflow beneath the canal from the west. Sensitivity of the model results to the flux from the western 

model boundary was tested because the slurry wall could not be explicitly simulated. 

The model was modified, using the MODFLOW-2005 Horizontal Flow Barrier (HFB) Package 

to cut-off flow from the western boundary, to evaluate the effect on capture zones and drawdown near

the WWF. Using the HFB package, interfaces between the L-31N RIV cells and the cells east of these 

cells were assigned a barrier with hydraulic conductivity of zero to eliminate flow from the west 

(fig. 44). The effects on the capture zones and drawdown of the WWF were evaluated using a Monte 

Carlo analysis for the steady-state dry conditions. 

A proposed recharge canal west of the NWWF along the Dade-Broward Levee was also 

included in this scenario (fig. 44). This canal is proposed to be connected to a northern adjoining 

interceptor canal and will recharge 60 Mgal/d from an external source of water. The recharge canal was 

represented in the steady-state seepage mitigation simulation with the WEL Package. The 219 cells 

along the length of this canal each contained a well that injected at a rate of 0.27 Mgal/d. The 

MODPATH 5.0 option IFACE for particle tracking for these cells was set to 6 to ensure that particles 

could track backwards to the top of the cell to properly represent seepage from the recharge canal.

In the NWWF, the 210-day 95-percent CI and median capture zones differ from the dry 

conditions simulation that did not implement seepage mitigation (fig. 45). The presence of the recharge 

canal moves the extent of the simulated 210-day capture zone toward the east, away from the model 

boundary, relative to the 210-day capture zone in the absence of the recharge canal. Also shown in this 

figure are particle endpoints at the end of the steady-state Monte Carlo simulations (where particles 
30

eventually track backwards to water source locations sometime after 210 days). The endpoints 

represent the source areas within the simulated area for the NWWF under this seepage mitigation 

scenario. In the absence of the recharge canal, particles tracked backwards to the existing canals. With 

the introduction of seepage mitigation measures, some particles track backwards to the proposed 

recharge canal, although the majority of the particles track backwards to the western boundary. As this 

western boundary represents a composite flux of water from the L-30 Canal and west of the canal, 

particles that track backwards to this boundary may represent source-area water west of this canal. 

Results for the 100-, 30-, and 10-day capture zones and drawdowns for the seepage mitigation scenario 

(not shown) were not substantially different from those for the dry conditions without the slurry wall 

and recharge canal (figs. 28, 29, and 30). 

Results of this simulation also indicate that the imposed flow barrier affects the extent of 

simulated capture zones around the WWF (figs. 46 and 47). The 95-percent CI upper bounds for the 

210- and 100-day capture zones (figs. 46 and 47) for the seepage mitigation scenario expanded eastward 

and toward the north and south, and retracted from the west. Reduction of flow from the western 

boundary causes the WWF to capture more flow from all other directions, thus expanding the capture 

zone extent in those directions. 

Proposed Scenario 2: Investigation of Effects of Lake Expansion

Substantial expansion of rock-mining activities in the Lake Belt area has been proposed that 

would increase the surface area of quarry lakes by an additional 25 mi2 (fig. 48) (a 150-percent increase 

from 2004 conditions). This proposed expansion may have a measurable effect on particle travel times, 

the areal extent of capture zones, and drawdown for both the NWWF and WWF. The effect of this lake 

expansion was simulated stochastically assuming representative steady-state dry (June 19, 2004) 

hydrologic conditions. 

The 210-day median and 95-percent CI upper bound for these simulations show a marked 

decrease in the areal extent for the travel-time capture zones of the NWWF (fig. 49), relative to the 

results of the dry conditions simulations. This is attributed largely to the increase in volume of these 

lakes in the subsurface near the NWWF. In the dry conditions simulations, backward-tracked particles 

moving to the west of the NWWF may have traversed no lakes or one lake; in the proposed lakeexpansion simulation, nearly all NWWF backward-tracked particles reach lakes in all directions from 

the well field. The bounds of the 95-percent CI for the NWWF are difficult to distinguish from the 

210-day capture zone median because once a particle enters a lake, it moves to the inflow location for 

the duration of the residence time. The extent of the simulated 0.25-ft drawdown contour was also 

reduced relative to simulations of the dry conditions without additional lakes, although not as 

substantially as the reduction of the capture zone extent. 

Around the WWF, the 95-percent CI upper bound for the 210-day capture zone for steady-state 

dry conditions for the lake-expansion scenario is similar to results without additional lakes, except for 

the area west and northwest of the WWF (fig. 50). Expansion of the largest lake nearly to the L-31N 

Canal shifted the capture zone west of its location without lake expansion (fig. 50). The 0.1-ft 

drawdown decreased in extent near the largest lake northwest of the WWF, but otherwise remained 

unchanged throughout the well-field area when compared to the results from the dry conditions 

simulations (fig. 50). 
31

Model Limitations 

The groundwater flow model and Monte Carlo simulations developed for this investigation were 

designed to estimate well-field flow regimes and travel times and were used to delineate capture zones 

and drawdown associated with the NWWF and the WWF. The model and methodology were designed 

to address specific objectives and are not necessarily designed to address other problems or answer 

other questions in the study area, without consideration of the technical assumptions made for this 

particular investigation. 

Model errors can mainly be attributed to limitations of the conceptual model and inaccuracies in 

representation of the aquifer input parameters and boundary conditions. Aquifer properties vary 

throughout Miami-Dade County and limited data (hydraulic conductivity, canal conductances, storage 

properties, and effective porosity) are available to characterize hydraulic parameters in each model cell. 

Even measured parameters are associated with varying degrees of uncertainty from measurement error, 

which can be considerable in some cases. Transmissivity estimates derived from aquifer performance 

tests can range over orders of magnitude depending on the type of test, equipment, aquifer 

heterogeneity, method of analysis, and conditions during the test. Hydraulic conductivity derived from 

transmissivity estimates requires an assumption of the thickness of the contributing zone, which is not 

known in all locations where transmissivity data are derived from field tests. The conceptual model 

developed in this investigation was based on available estimates for these input parameters, with a 

built-in level of measurement error in some instances (implemented by use of a nugget in the variogram 

for hydraulic conductivity realizations) to stochastically estimate well-field capture zones and 

drawdown and quantify the uncertainty. Capture zone estimates are sensitive to hydraulic conductivity 

and effective porosity values and are therefore sensitive to any potential errors in estimation of these 

values that are incorporated in the 10,000 realizations used in the Monte Carlo analysis. 

MODFLOW-2005 is designed for simulating flow in aquifers under laminar, non-turbulent flow 

regimes and is based on the assumptions of Darcy’s law. In fact, the karstic Biscayne aquifer is known 

to be highly heterogeneous and anisotropic, with the presence of high-porosity, high permeability, flow 

zones near the NWWF and WWF (Shoemaker and others, 2008a). A limitation of the model is that it 

does not account for the possible existence of turbulent, non-Darcian flow, especially near the well 

fields during changes in pumping conditions. With turbulent flow, the linear relation between specific 

discharge and the hydraulic gradient in Darcy’s law may not apply at high velocity. Turbulence was not 

considered because steady-state models were used in this analysis.

The spatial distribution and variability of hydraulic properties generally contribute substantial 

uncertainty in hydrologic models. For this model, hydraulic conductivity values are derived from 

limited number of APTs at specific locations. These values vary over several orders of magnitude and 

may differ over short distances. Hydraulic conductivity in areas lacking data is uncertain. Furthermore, 

estimated values of hydraulic conductivity are subjected to error in measured field parameters and 

assumptions made about the flow system when calculating transmissivity from field test data. Derived 

conductivity values usually represent a composite value for the entire aquifer thickness, although 

observations of stratiform flow zones indicate vertical heterogeneity in conductivity. Effective porosity 

cannot be measured directly, and may be considerably different than measured porosity, and there are 

few transport data with which to calibrate effective porosity. Effective porosity may also be spatially 

variable and scale dependent. Uncertainty in hydraulic conductivity (including values used as canal 

conductances) and effective porosity in the model were addressed using the Monte Carlo approach. The 

Monte Carlo simulation results, although accounting for a reasonable range of variation of the lateral 

distribution of hydraulic conductivity and effective porosity values, may also be sensitive to the nugget 
32

and sill values specified in the variogram model used to generate the 10,000 individual fields of 

hydraulic conductivity. 

Although multiple flow zones are known to exist at some locations in the Biscayne aquifer, a 

one-layer model was used in this investigation because (1) the distribution and connectivity of 

individual flow zones across the model area is unknown, (2) few water-level data for individual flow 

zones exist in the study area with which to constrain model calibration, (3) to attempt to simulate 

individual flow zones would introduce substantially more uncertainty to the model, and (4) a 

simplification of the system made the stochastic approach more tractable. Development and calibration 

of a three-dimensional groundwater flow model that represents multiple flow zones is not possible based 

on present understanding of flow zone continuity and hydraulic conductivity. Observed datasets are not 

available to calibrate individual flow zones, and this precludes formal quantification of uncertainty in 

capture zones. The flow zones have not been mapped over the entire study area; there are only a few 

monitoring wells that measure water levels of specific flow zones; and few APT data exist in the study 

area to constrain vertical and horizontal hydraulic conductivity estimates of layers representing 

individual flow zones. The effect of the presence of flow zones on the extent of travel-time capture 

zones was implicitly included in the analysis by using a large range of effective porosity values. 

Although analyses indicated the travel-time capture zones and drawdown for a one-layer model differ 

somewhat from those for a multi-layer model (Appendix 1), one-layer models with a composite 

hydraulic conductivity that represent multiple lithologies have been used to simulate flow for the 

surficial aquifer system of Florida at a regional-scale (Sepulveda and Kuniansky, 2010; Sepulveda and 

others, 2012). While simulating multiple layers may improve precision of particle tracking, the decision 

was made to simplify representation of the Biscayne aquifer with a one-layer model to limit uncertainty, 

making the stochastic approach more tractable. 

The model developed in this study may not be appropriate to address questions regarding 

ultimate sources of groundwater to the well fields from west of the model boundary. Many of the Monte 

Carlo simulated particle tracks for longer travel times terminate in the western boundary of the model 

(L-30 and L-31N Canals). These hydrologic features separate the urban part of Miami-Dade County 

from the Everglades National Park and Water Conservation Area 3B, yet studies indicate that underflow 

may occur beneath the canals from the western wetland areas (Cunningham and others, 2006; 

Sonenshein, 2001). Field tests and vertical-head data are required to better characterize the temporal and 

spatial distribution of underflow. A primary objective of this study was to determine capture zones 

within the urban part of the study area. Therefore, detailed simulation of contributing areas within the 

Everglades and the Water Conservation Areas was not within the scope of this study. The headdependent boundary conditions at the western model boundary are designed to represent all flow 

coming from the canals and under the canals from the Everglades and Water Conservation Areas. 

Particles that track backwards to this western model boundary indicate that the source area could be 

either the canal or west of the canal. 

Simplifications were made to boundary conditions that affect recharge and discharge. Areal 

recharge is a function of precipitation and evapotranspiration. Precipitation and evapotranspiration 

rates were applied uniformly within the model domain, although they are known to be spatially variable. 

Canal stages and conductances of canal beds also affect recharge and discharge to the system. Errors 

in canal stage values are possible, particularly for smaller secondary canals where stage data were not 

directly available, or in ungaged reaches. Canal bed conductances, while calibrated, are largely 

unknown, and contribute uncertainty to the calibrated recharge. Explicit simulation of runoff and 

dynamic surface-water flow could improve the distribution of canal stages and thus the estimates 

of recharge, but were deemed beyond the scope of this study. Results of these analyses should be 
33

considered in light of the considerable uncertainty in the controls on recharge and discharge in 

the system. 

Mixing within lakes was not explicitly simulated in the model. Instead, a lake-residence time 

was incorporated, using an approach modified from Chin and others (2010), to estimate particle 

transport through lakes for the Monte Carlo analysis (Appendix 6). This adjustment was based on the 

assumption of lake mixing occurring over the period of the residence time, which may not be a valid 

assumption as the residence time is a measure of the average amount of time that water is retained in a 

lake. Through the Monte Carlo analysis, a range of residence times was calculated for each lake, based 

on the flow rate from the lake to the aquifer, which changed between each realization. Furthermore, the 

pathlines of backward-tracked particles through lakes are different between each realization. Therefore, 

although a simple residence-time approach to lake transport was utilized, the Monte Carlo analysis 

considers the uncertainty of particle pathlines and lake-aquifer exchanges, and therefore, ensures that 

the greatest possible capture area is estimated. 

A further assumption made in the modeling process is that groundwater recharge from rainfall 

occurs instantaneously. The unsaturated zone is not explicitly simulated because delay of infiltrating 

rainfall in the unsaturated zone is expected to be less than a day and is not expected to substantially 

affect particle transport compared to other water-budget components. 

The no-flow boundary condition at the bottom of the model assumes no water is exchanged 

between the lower Biscayne aquifer and the underlying gray limestone aquifer as there is a 

low-permeability confining unit that limits leakage between the two aquifers (Fish and Stewart, 1991). 

This no-flow boundary is also used in previous models (Langevin, 2001; Lohmann and others, 2012). 

There is likely some heterogeneity in the properties of the confining unit and leakage may occur but 

leakage is considered to be small relative to other sources/sinks in the model. 

Particle tracking evaluations were made for steady-state conditions using MODPATH, which 

accounts for advective transport, but does not consider dispersive transport. Advective particle tracking 

is commonly used to determine well-field capture zones and was used to develop the existing capture 

zones for the study area. Results of tracer tests near the NWWF indicate that dispersive transport is 

estimated to be small compared to advective transport at a scale of about 100 m (Renken and others, 

2008; Shapiro and others, 2008). These results may not apply at the scale of this investigation, 

nonetheless the assumption was made that dispersion contributed minimally to the reach of individual 

pathlines. MODPATH can simulate particle tracking in both steady-state and transient flow fields, but 

results from transient particle tracking were assumed to fall within the results from the steady-state

Monte Carlo analysis. Backward particle tracking simulations were analyzed with the Monte Carlo 

approach; forward particle tracking simulations were not analyzed, except as discussed in Appendix 1 to 

verify backward-tracked particle pathlines. It is possible that forward tracking my provide different 

results than backward tracking, but the results in Appendix 1 indicate that, in this case, the results would 

be essentially the same.

Simulated results of how the system would respond to proposed modifications should be 

interpreted with caution because the model was not designed specifically to evaluate those conditions. 

Groundwater recharge from a proposed recharge canal was modeled using the MODFLOW Well 

(WEL) Package, allowing direct input of a uniform rate of water entering the aquifer along the entire 

length of the proposed canal. In reality, recharge from the canal may not occur at this maximum 

specified rate and may not be controllable, and also may not be equally distributed along the length of 

the canal, as it is in the model. 

The proposed slurry wall could not be explicitly simulated because (1) the location of the 

proposed structure is beyond the model boundary, and (2) the one-layer model cannot directly account 
34

for a flow barrier that is installed only partially through the aquifer. Application of a horizontal flow 

barrier on the east side of the western head-dependent flux boundary effectively cut off flow from the 

west, and capture-zone extents and drawdowns were shown to be sensitive to this change in the flow 

field. 

Sensitivity analyses are used to determine the impact that parameter uncertainty may have on 

model predictions. A quantitative predictive sensitivity analysis to evaluate to what extent individual 

parameters and stresses affect the capture zones was not formally conducted. Instead, the Monte Carlo 

analysis incorporates uncertainty in the hydraulic properties, and was used in this investigation to 

investigate parameter uncertainty effects on capture zones and drawdown for steady-state wet, average 

and dry hydrologic conditions. 

Summary and Conclusions 

The U.S. Geological Survey investigated groundwater flow paths in the urban areas near the 

Northwest (NWWF) and West (WWF) Well Fields in Miami-Dade County, Florida, to evaluate traveltime capture zones and drawdowns. A one-layer groundwater flow model of the Biscayne aquifer was 

developed for the eastern part of Miami-Dade County. The model documented herein was modified 

from a previously existing integrated surface-water/variable-density groundwater flow model (Lohmann 

and others, 2012). The model accounts for daily changes in interactions between the aquifer and 

surface-water features, recharge and evapotranspiration, changing tidal stages in Biscayne Bay, as well 

as groundwater pumping. The model was calibrated to transient daily conditions during January 1996 to 

December 2004 by adjustment of evapotranspiration extinction depth, ungaged secondary and tertiary 

canal stages, canal conductance, and hydraulic conductivity. Calibration targets were observed 

groundwater levels and net exchanges between the canals and aquifer in primary canal basins. 

The model was modified to represent steady-state dry, average and wet hydrologic conditions to 

determine travel-time capture zones and drawdown for the NWWF and WWF. Canal stage data from 

two control structures near the NWWF and WWF were used to select one day each representing dry, 

average, and wet hydrologic conditions in the study area. Model input for these 3 days were selected 

from the transient model, and adjustments were made to annualize the rainfall rates for each of these 

hydrologic conditions. To simulate transport more accurately through the lakes, bulk hydraulic 

conductivities were calculated from information on lake depths and aquifer thicknesses and applied to 

model cells with lakes. Using the annualized rainfall rates and bulk hydraulic conductivities for lakes, 

three steady-state simulations representing dry, average, and wet hydrologic conditions were run with 

and without well-field pumping. 

To account for uncertainty in the spatial distributions of hydraulic conductivity, effective 

porosity, and canal conductances, 10,000 stochastic realizations of these model input parameters were 

created based on published data from field tests and documented models. For each of the three 

hydrologic conditions, the 10,000 realizations were simulated, and the 10,000 results were used to 

calculate drawdown and to simulate particle pathways. The median and 95-percent confidence intervals 

(CI) of the 0.25-foot (ft) and 0.1-ft drawdowns were contoured for each of the three hydrologic 

conditions. MODPATH 5.0 backward particle tracking was applied for each simulation, and particle 

endpoints through lakes were further adjusted to account for the lake residence time, to estimate the 

median and 95-percent CIs of the 10-, 30-, 100-, and 210- day capture zones for each of the three 

hydrologic conditions. Each particle path through a lake was adjusted to account for lake residence time 

without explicitly simulating mixing in the lake. 
35

Two scenarios were simulated to investigate the effects on capture zones and drawdown of 

(1) two seepage mitigation structures, and (2) future proposed expansion of man-made quarry lakes. The 

mitigation scenario simulated the addition of a recharge canal near the NWWF and a flow barrier west 

of the WWF, effectively cutting off flow from the western boundary. The lake-expansion scenario 

simulated increasing the spatial extent of the lakes by 25 square miles (mi2). 

The results of this investigation indicate several considerations for management of the 

groundwater resources near the NWWF and WWF. 

• Because the ambient groundwater-flow direction is generally from west to east near the well fields, 

the capture zones are more extensive toward the west than east of the well fields for all travel times 

and hydrologic conditions evaluated, particularly for the NWWF. The shape of the capture zones in 

areas with lakes is more irregular than in areas without lakes, because of the residence-time 

adjustment in lakes.

• In general, the deterministic and median and 95-percent CIs of the 0.25-ft and 0.1-ft drawdown 

contours differ minimally between simulations for the three hydrologic conditions, and the median 

CIs for the 0.25-ft and 0.1-ft drawdown contours are similar to the simulated drawdown contours for 

the deterministic base-case simulations. 

• The capture zones for the NWWF show the same basic shape for dry, average, and wet hydrologic 

conditions, with slightly greater extent toward the west for drier conditions than for average and wet 

conditions. 

• The simulated capture zones for the WWF indicate that recharge areas for this well differ slightly for 

dry and wet hydrologic conditions. The WWF capture zones, in general, are radially symmetric 

around the well field during dry hydrologic conditions, except where lakes are present west of the 

well field, but are slightly skewed toward the north during wet hydrologic conditions. The 30-day 

capture zone for the WWF extends into a lake southwest of the well field during dry conditions, 

indicating the possibility of faster transport of an introduced-chemical constituent from this lake to 

the WWF during the dry season than during the wet season. 

• The 10-day travel-time capture zone indicates an area within which introduced-chemical 

constituents may reach well fields in a short amount of time. One lake, just north of the NWWF, is 

included within the 10-day capture zone for all three hydrologic scenario runs. Mixing processes in 

lakes are likely more complicated than represented herein and therefore, there is substantial 

uncertainty whether an introduced-chemical constituent in a lake will reach a well field within 10 

days if the lake is within a 10-day capture zone, yet the possibility exists. 

• In general, backward-tracked particles ultimately terminate at canals, even when the weak-sink 

particle-tracking option in MODPATH 5.0 is used, indicating that most canals in the study area act 

as strong sources of water. Westward-tracking particles can originate no further than the western 

model boundary, which is defined by the L-31N and L-30 Canals. Because the western boundary 

condition represents a composite flux from the canal and from groundwater entering the model 

domain beneath the canal from the west, particles that track backwards to this boundary could 

represent water originating from sources west of the boundary, from the canals themselves, or a 

combination of western sources and canals. 
36

• The composite 210-day capture zones differ substantially in shape and extent from the 210-day 

capture zones used by the county to establish Wellhead Protection Areas (WHPAs) because (1) the 

model described in this report has accounted for the presence and transport of particles through 

lakes, (2) a stochastic approach was used for these evaluations to determine 95-percent confidence 

intervals that accounted for uncertainty in hydraulic property distribution, and (3) the pumping rate 

used for the WWF in this analysis was substantially lower than that used to establish the original 

WHPAs. 

• Alterations to the hydrologic system may affect the distribution of capture zones and drawdown 

contours. Simulation results show that (1) the recharge canal captures some particles that would 

have otherwise tracked back to the western boundary, and reduces the extent of the capture zones; 

(2) reduction of flow from the western boundary causes the WWF to capture more flow from all 

other directions, thus expanding the capture zone and drawdown extent in those directions; and (3) 

the addition of lakes resulted in a general decrease in the areal extent of capture zones and a 

decrease in the extent of the 0.25-ft drawdown contour for the NWWF, expanded the westward 

extent of the capture zone around the WWF, and decreased the extent of the drawdown contour of 

the WWF. 

• The model is limited by assumptions about the flow system, simplifications made in the conceptual 

model, and knowledge of the hydraulic properties of the units simulated. These include, among 

others: the distribution of hydraulic properties; the simulation of the Biscayne aquifer as a single 

composite layer; the proximity of the study area to the western model boundary; and the transport 

behavior in quarry lakes. While model results should be considered in light of these limitations, the 

Monte Carlo analysis incorporates uncertainty into the model predictions of capture zones and 

drawdown. This allows for calculation of confidence intervals of results, producing a more 

informative evaluation of the response of the hydrologic system to varying stresses and hydrologic 

conditions than an evaluation based solely on deterministic simulation results.
37

Selected References 

Anderson, M.P., Hunt, R.J., Krohelski, J.T., and Chung, K., 2002, Using High Hydraulic Conductivity 

Nodes to Simulate Seepage Lakes: Ground Water, v. 40, no. 2, p. 117–122. 

Anderson, M.P., and Woessner, W.W., 1992, Applied groundwater modeling, simulation of flow and 

advective transport: San Diego, Calif., Academic Press, Inc., 381 p. 

Bair, E.S., Safreed, C.M., and Stasny, E.A., 1991, A Monte Carlo-based approach for determining 

traveltime-related capture zones of wells using convex hulls as confidence regions: Ground Water, v. 

29, no. 6, p. 849–855. 

Barlow, P.M. 1997, Particle-tracking analysis of contributing areas of public-supply wells in simple and 

complex flow systems, Cape Cod, Massachusetts: U.S. Geological Survey Water-Supply Paper 2434, 

66 p. 

Camp Dresser & McKee Inc. (CDM), 1985, Groundwater flow model for the Northwest Well Field, 

Dade County, Florida: Consultant report prepared for Miami-Dade Department of Environmental 

Resources Management.

CH2MHill, 2001, Risk assessment and groundwater modeling of the Miami-Dade Northwest Wellfield: 

Report prepared for the Miami-Dade Department of Environmental Resources Management (DERM).

Chin, D.A., Iudicello, J.J., Kajder, K.C., Kelly, P.M., Porzilli, D.V., and Guha, H., 2010, Lake effect in 

wellhead protection: Journal of Water Resources Planning and Management, ASCE, May/June 2010. 

Chin, D.A., and Patterson, R.D., 2005, Quantification of hydrologic processes and assessment of 

rainfall-runoff models in Miami-Dade County, Florida: U.S. Geological Survey Scientific 

Investigations Report 2004–5191. 

Code of Miami-Dade County, Florida, 2009, Ordinance No. 09–27, accessed June 26, 2012, 

at http://library.municode.com/index.aspx?clientId=10620. 

Cooper, R.P., and Lane, J., 1987, An atlas of eastern Dade County surface water management basins: 

South Florida Water Management District Technical Memorandum. 

Cunningham, K.J., Carlson, J.L., Wingard, G.L., Robinson, Edward, and Wacker, M.A., 2004a, 

Characterization of aquifer heterogeneity using cyclostratigraphy and geophysical methods in the 

upper part of the karstic Biscayne aquifer, southeastern Florida: U.S. Geological Survey WaterResources Investigations Report 03–4208, 46 p. 

Cunningham, K.J., Renken, R.A., Wacker, M.A., Zygnerski, M.R., Robinson, E., Shapiro, A.M., and 

Wingard, G.L., 2006a, Application of carbonate cyclostratigraphy and borehole geophysics to 

delineate porosity and preferential low in the karst limestone of the Biscayne aquifer, SE Florida, in

Harmon, R.S., and Wicks, C., eds., Perspectives on Karst Geomorphology, Hydrology, and 

Geochemistry: Geological Society of America, v. 404, p. 191–208. 

Cunningham, K.J., and Sukop, M.C., 2011, Multiple technologies applied to characterization of the 

porosity and permeability of the Biscayne aquifer, Florida: U.S. Geological Survey Open-File Report 

2011–1037, 8 p. 

Cunningham, K.J., Sukop, M.C., Huang, H., Alvarez, P.F., Curran, H.A., Renken, R.A., and Dixon, J.F., 

2009, Prominence of ichnologically influenced macroporosity in the karst Biscayne aquifer: 

Stratiform "super-K" zones: Geological Society of America Bulletin, v. 121, no. 1/2, p. 164–180. 
38

Cunningham, K.J., Wacker, M.A., Robinson, E., Dixon, J.F., and Wingard G.L., 2006b, A 

cyclostratigraphic and borehole geophysical approach to development of a three-dimensional 

conceptual hydrogeologic model of the karst Biscayne aquifer, southeastern Florida., U.S. Geological 

Survey Scientific Investigations Report, 2005–5235, 69 p. 

Cunningham, K.J., Wacker, M.A., Robinson, Edward, Gefvert, C.J., and Krupa, S.L., 2004b, 

Hydrogeology and ground-water flow at Levee-31N, Miami-Dade County, Florida, July 2003 to May 

2004: U.S. Geological Survey Scientific Investigations Map I–2846, 1 sheet. 

Department of Environmental Resources Management, 1992, Modeling for the West Wellfield updated 

protection area: Miami-Dade County Department of Environmental Resources Management 

Technical Memorandum. 

Deutsch, C.V., and Journel, A.G., 1998, GSLIB: Geostatistical software library: Oxford University 

Press, 384 p. 

Doherty, J., 2010, PEST: Model independent paramater estimation, user documentation: Brisbane, 

Queensland, Australia, Watermark Numerical Computing. 

Federal Register Notice, 1979, Biscayne aquifer; notice of determination: Notice FRL 1264–3, v. 44, no. 

198, p. 58797–58798. 

Fish, J.E., and Stewart, M.T., 1991, Hydrogeology of the surficial aquifer system, Dade County, 

Florida: U.S. Geological Survey Water-Resources Investigations Report 90–4108. 

Florea, L.J., Cunningham, K.J., and Altobelli, S., 2009, NMR imaging of fluid exchange between 

macropores and matrix in eogenetic karst: Groundwater, v. 47, p. 382–390. 

Franke, O.L., Reilly, T.E., Pollock, D.W., and LaBaugh, J.W., 1998, Estimating areas contributing 

recharge to wells, lessons from previous studies: U.S. Geological Survey Circular 1174, 14 p. 

Franzetti, S., and Guadagnini, A., 1996, Probabilistic estimation of well catchments in heterogenous 

aquifers: Journal of Hydrology, v. 174, p. 149–171. 

German, E.R., 2000, Regional Evaluation of Evapotranspiration in the Everglades: U.S. Geological 

Survey Water-Resources Investigations Report 00–4217, 48 p. 

Giddings, J.B., Kuebler, L.L., Restrepo, J.I., Rodberg, K.A., Montoya, A.M., and Radin, H.A., 2006, 

Draft report of Lower East Coast subRegional (LECsR) MODFLOW Model Documentation, South 

Florida Water Management District, 294 p. 

Guha, H., 2008, A stochastic modeling approach to address hydrogeologic uncertainties in modeling 

wellhead protection boundaries in karst aquifers: Journal of the American Water Resources 

Association, v. 44, no. 3. 

Harbaugh, A.W., 1990, A computer program for calculating subregional water budgets using results 

from the U.S. Geological Survey modular three-dimensional ground-water flow model: U.S. 

Geological Survey Open-File Report 90–392, 46 p. 

Harbaugh, A.W., 2005, Modflow–2005, The U.S. Geological Survey modular ground-water model—

The Ground-Water Flow Process, U.S. Geological Survey Techniques and Methods 6–A16. 

Harvey, R.W., Metge, D.W., Shapiro, A.M., Renken, R.A., Osborn, C.L., Ryan, J.N., 

 Cunningham, K.J., and Landkamer, L.L., 2008, Pathogen and chemical transport 
39

 in the karst limestone of the Biscayne aquifer: 3. Use of microspheres to estimate the 

 transport potential of Cryptosporidium parvum oocysts: Water Resources Research, v. 

 44, W08431, doi:10.1029/2007WR006060. 

Hunt, R.J., Haitjema, H.M., Krohelski, J.T., and Feinstein, D.T., 2003, Simulating ground-water-lake 

interactions: Approaches and insights: Ground Water, v. 41, no. 2, p. 227–237. 

Kerans and Tinker, 1997, Sequence stratigraphy and characterization of carbonate reservoirs: SEPM 

Short course notes, Vol. 40, 130 p. 

Kuniansky, E.L., Gómez-Gómez, Fernando, and Torres-González, Sigfredo, 2004, Effects of aquifer 

development and changes in irrigation practices on ground-water availability in the Santa Isabel area, 

Puerto Rico: U.S. Geological Survey Water-Resources Investigations Report 03–4303, 56 p. 

Lake Belt Mitigation Committee, 2009, Lake Belt Mitigation Committee meeting summary, April 17, 

2009: Accessed March 9, 2010, 

at http://www.sfwmd.gov/portal/page/portal/pg_grp_sfwmd_watersupply/subtabs%20-

%20miamidadeco%20%20-%20mitigation/tab1610177/meeting%20summary4-17-09.pdf

Langevin, C.D., 2001, Simulation of ground-water discharge to Biscayne Bay, Southeastern Florida: 

U.S. Geological Survey Water-Resources Investigations Report 00–4251. 

Lee, T.M., 1996, Hydrogeologic controls on the groundwater interactions with an acidic lake in karst 

terrain, Lake Barco, Florida: Water Resources Research, v. 32, no. 4, p. 831–844. 

Lohmann, M.A., Swain, E.D., Wang, J.D., and Dixon, J., 2012, Evaluation of Effects of Changes in 

Canal Management and Precipitation Patterns on Salinity in Biscayne Bay, Florida, using an 

Integrated Surface-Water/Groundwater Model: U.S. Geological Survey Scientific Investigations 

Report 2012-5099, 94 p. 

Marella, R.L., 2009, Water withdrawals, use, and trends in Florida, 2005: U.S. Geological Survey 

Scientific Investigations Report 2009–5125, 49 p. 

McDonald, M.G., and Harbaugh, A.W., 1988, A modular three-dimensional finite-difference groundwater flow model: U.S. Geological Survey Techniques of Water Resources Investigations, book 6, 

chap. A1, 586 p. 

Merritt, M.L., 1997, Simulation of the water-table altitude in the Biscayne aquifer, Southern Dade 

County, Florida, water years 1945–89: U.S. Geological Survey Water-Supply Paper 2458, 148 p. 

Merritt, M.L., 1996, Numerical simulation of a plume of brackish water in the Biscayne aquifer 

originating from a flowing artesian well, Dade County, Florida: U.S. Geological Survey WaterSupply Paper 2464, 74 p. 

Merritt, M.L. and Konikow, L.F., 2000, Documentation of a computer program to simulate lake-aquifer 

interaction using the MODFLOW ground-water flow model and the MOC3D solute-transport model: 

U.S. Geological Survey Water-Resources Investigations Report 00–4167, 146 p. 

Parker, G.G., 1951, Geologic and hydrologic factors in the perennial yield of the Biscayne aquifer: 

American Water Works Association Journal, v. 8, p. 817–835. 

Pollock, D.W., 1994, User's Guide for MODPATH/MODPATH-PLOT, Version 3: A particle tracking 

post-processing package for MODFLOW, the U.S. Geological Survey finite-difference ground-water 

flow model: U.S. Geological Survey Open-File Report 94–464, 6 ch. 
40

Prickett, T.A., and Lonnquist, C.G., 1971, Selected digital computer techniques for groundwater 

resource evaluation: Illinois State Water Survey Bulletin 55, Champaign, Ill. 

Reese, R.S. and Cunningham, K.J., 2000, Hydrogeology of the gray limestone aquifer in southern 

Florida: U.S. Geological Survey Water-Resources Investigations Report 99–4213, 244 p. 

Reese, R.S., and Richardson, Emily, 2008, Synthesis of the hydrogeologic framework of the Floridan 

aquifer system and delineation of a major Avon Park permeable zone in central and southern Florida: 

U.S. Geological Survey Scientific Investigations Report 2007–5207, 60 p., 4 pls., plus apps. 

Reilly, T.E. and Pollock, D.W., 1993, Factors affecting areas contributing recharge to wells in shallow 

aquifers: U.S. Geological Survey Water Supply Paper 2412, 21 p. 

Renken, R.A., Cunningham, K.J., Shapiro, A.M., Harvey, R.W., Zygnerski, M.R., Metge, D.W., and 

Wacker, M.A., 2008, Pathogen and chemical transport in the karst limestone of the Biscayne aquifer: 

1. Revised conceptualization of groundwater flow: Water Resources Research, v. 44, W08429, doi: 

10.1029/2007WR006058 

Renken, R.A., Shapiro, A.M., Cunningham, K.J., Harvey, R.W., Metge, D.W., Zygnerski, M.R., 

Osborn, C.L., Wacker, M.A., and Ryan, J.N., 2005, Assessing the vulnerability of a municipal well 

field to contamination in a karst aquifer: Environmental and Engineering Geoscience, v. 11, no. 4, p. 

319-331. 

Sepúlveda, Nicasio and Kuniansky, Eve L., 2010, Effects of model layer simplification using composite 

hydraulic properties: Hydrogeology Journal, v. 18, p. 405-416. 

Sepúlveda, Nicasio, Tiedeman, C.R., O’Reilly, A.M., Davis, J.B., and Burger, Patrick, 2012, 

Groundwater flow and water budget in the surficial and Floridan aquifer systems in east-central 

Florida: U.S. Geological Survey Scientific Investigations Report 2012–5161, 214 p. 

Shapiro, A.M., Renken, R.A., Harvey, R.W., Zygnerski, M.R., and Metge, D.W., 2008, Pathogen and 

chemical transport in thekarst limestone of the Biscayne aquifer: 2. Chemical retention from diffusion 

and slow advection: Water Resources Research, v. 44, W08430, doi:10.1029/2007WR006059. 

Shoemaker, W.B., Cunningham, K.J., Kuniansky, E.L., and Dixon, J.F., 2008a, Effects of 

 turbulence on hydraulic heads and parameter sensitivities in preferential groundwater 

flow layers: Water Resources Research, v. 44, W03501, doi:10.1029/2007WR006601. 

Shoemaker, W.B., Kuniansky, E.L., Birk, Steffen, Bauer, Sebastian, and Swain, E.D., 2008b, 

Documentation of a Conduit Flow Process (CFP) for MODFLOW-2005: U.S. Geological Survey 

Techniques and Methods 6-A24, 50 p. 

Solo-Gabriele, H., 2000, An isotopic study of two rock mining lakes: Second interim technical report 

for isotope study, Phase I, Submitted to Kendall Properties and Investments, Inc., University of Miami 

Departmant of Civil, Architectural, and Environmental Engineering, Coral Gables, Fla. 

Sonenshein, R.S., 2001, Methods to quantify seepage beneath Levee 30, Miami-Dade County, Florida: 

U.S. Geological Survey Water-Resources Investigations Report 01-4074, 36 p. 

Sonenshein, R.S., and Koszalka, E.J., 1996, Trends in water-table altitude (1984–93) and saltwater 

intrusion (1974–93) in the Biscayne aquifer, Dade County, Florida: U.S. Geological Survey Open-File 

Report 95–705, 2 sheets. 
41

Snyder, D.T., Wilkinson, J.M., and Orzol, L.L., 1998, Use of a ground-water flow model with particle 

tracking to evaluate ground-water vulnerability, Clark County, Washington: U.S. Geological Survey 

Water-Supply Paper 2488, 63 p. 

Starn, J.J., Bagtzoglou, A.C., and Robbins, G.A., 2010, Using atmospheric tracers to reduce uncertainty 

in groundwater recharge areas: Ground Water v. 48, no. 6, p. 858–868, doi:10.1111/j.1745–

6584.2010.00674.x. 

Starn, J. Jeffrey and Brown, Craig J., 2007. Simulations of Ground-Water Flow and Residence Time 

near Woodbury, Connecticut. Geological Survey (U.S.) Scientific Investigations Report 2007-5210, 

viii, 45 p. 

Starn, J.J., Stone, J.R., and Mullaney, J.R., 2000, Delineation and analysis of uncertainty of contributing 

areas to wells at the Southbury Training School, Southbury, Connecticut: U.S. Geological Survey 

Water-Resources Investigation Report 00–4158, 53 p. 

Swain, E.D., Howie, B., and Dixon, J., 1996, Description and field analysis of a coupled ground 

water/surface-water flow model (MODFLOW/BRANCH) with modifications for structures and 

wetlands in southern Dade County, Florida: U.S. Geological Survey Water-Resources Investigations 

Report 96–4118, 71 p. 

U.S. Environmental Protection Agency, 1987, Guidelines for delineation of wellhead protection areas: 

Office of Ground-Water Protection, U.S. Environmental Protection Agency report no. 440687010. 

Varljen, M.D., and Shafer, J.M., 1991, Assessment of uncertainty in time-related capture zones using 

conditional simulation of hydraulic conductivity: Ground Water, v. 29, no. 5, p. 737–748. 

Wilson, J.D. and Naff, R.L., 2004, The U.S. Geological Survey modular ground-water model—GMG 

linear equation solver package documentation: U.S. Geological Survey Open-File Report 2004–1261, 

47 p. 

43

Appendix 1. Comparative analysis of a one-layer and a multi-layer 

representation of a hydrologic system 

A conceptual exercise was developed to evaluate whether a one-layer groundwater flow model 

that uses bulk aquifer properties can accurately estimate capture zones and drawdown for a 

heterogeneous, anisotropic aquifer composed of two groundwater flow zones and a relatively less 

permeable zone in between. A single-layer steady-state MODFLOW groundwater flow model 

(conceptual Model A) was created with a well located in the center of the domain and constant-head 

boundaries on the left- and right-hand boundaries (fig. 1-1). The grid consisted of 100 rows and 100 

columns with 1,640-feet (ft) (500-meter (m)) grid spacing in the x- and y-directions. The top altitude of 

the model was assigned a uniform value of 0.0 ft (0.0 m) and the bottom altitude was assigned a uniform 

value of -69.0 ft (-21.0 m), similar to the surface and average bottom altitudes of the Biscayne aquifer 

(fig. 4). Constant-head cell boundaries were set to 0.0 IW. 

A second conceptual steady-state MODFLOW groundwater flow model (conceptual Model B) 

was developed containing three layers (fig. 1-1). A uniform horizontal grid discretization was used and 

specified constant-head cells on the left- and right-hand boundaries of the model in layer 1 were the 

same as those in Model A. The uppermost and lowermost layers of Model B represent continuous 

megaporosity flow zones mapped near the NWWF (Cunningham and others, 2006a, 2006b) and WWF. 

In this analysis, these flow zones were hypothesized to be well connected over the scale of the model 

domain. It was assumed that horizontal hydraulic conductivity in these two layers was considerably 

larger than in layer 2. This assumption corresponds well with a simplified conceptualization of the 

framework of the Biscayne aquifer near the NWWF and WWF (Cunningham and others, 2006; Renken 

and others, 2008) in which groundwater flow occurs primarily within the highly porous, touching-vug 

flow zones. The top and bottom altitudes of Model B were assigned the same values as those in 

conceptual Model A, with the top and bottom altitudes of the less permeable layer 2 set to -10 ft 

(-3.048 m) and -45 ft (13.716 m), respectively. 

Equivalent bulk hydraulic properties of horizontal hydraulic conductivity and effective porosity 

were used in the one-layer conceptual Model A, whereas uniform individual hydraulic properties were 

assigned to each layer in the multi-layer Model B. The vertical anisotropy in layers 1 and 3 of 

conceptual Model B was set at 1.0, with horizontal hydraulic conductivity being equivalent to vertical 

hydraulic conductivity. The vertical conductivity of layer 2 was one order of magnitude larger than the 

horizontal hydraulic conductivity of layer 2 and is consistent with air permeability data in Cunningham 

and others (2004a). The hydraulic conductivity was 32,808 feet per day (ft/d) (10,000 meters per day 

(m/d)), 328 ft/d (100 m/d), and 32,808 ft/d (10,000 m/d) in model layers 1, 2, and 3, respectively. An 

equivalent bulk hydraulic conductivity of 16,309 ft/d (4,971 m/d) was used in conceptual Model A and 

was calculated from the values of conductivities and thicknesses of the layers in conceptual Model B 

using a thickness-weighted average. The effective porosity for layers 1, 2, and 3 of conceptual Model B 

were set to 0.5, 0.176, and 0.5, respectively, to represent high-porosity flow zones in layers 1 and 3. The 

bulk equivalent effective porosity for conceptual Model A was calculated by multiplying the porosity of 

layers 1 and 3, which is assumed to dominate the flow, by the fraction of the thickness of the aquifer for 

which the groundwater flow zones compose, yielding an effective porosity of 0.246. 

Recharge and evapotranspiration were not applied in either model. In both models, a well was 

placed in a cell in the center of the model domain and assigned a specified flux in the MODFLOW-2005 

WEL package based on the maximum design-capacity pumping rate of the NWWF (225 million gallons 
44

per day (Mgal/d)). In the multi-layer conceptual Model B, the well was placed in the lowest highconductivity flow zone (layer 3). 

Particles were tracked backwards from the simulated well for 500 days, effectively delineating 

the 500-day travel-time capture zone. Particles were placed around the well in layer 3 of Model B, with 

particles distributed at evenly spaced vertical intervals throughout the layer. Particles in Model A were 

distributed at evenly spaced vertical intervals in the bottom third of the single model layer to remain 

consistent with the locations of particles in Model B. Results of backward-tracked particles are shown in 

figure 1-1. Particle pathlines and head contours are similar for both conceptual models indicating the

single-layer deterministic model of the Biscayne aquifer yields capture zones and drawdown 

comparable to those generated by a multi-layer model with explicit representation of flow zones. The 

parameter values used in Models A and B were representative of the Biscayne aquifer, however, 

considerable uncertainty exists for these parameters and deterministic results do not consider parameter

uncertainty.

To further evaluate the validity of the one-layer regional model, the stochastic approach used to 

address the effects of parameter uncertainty in the regional-model study was also applied to conceptual 

Models A and B. Horizontal hydraulic conductivity fields were generated using the sequential Gaussian 

simulation (SGSIM) method in the utility FIELDGEN (Doherty, 2010) to create stochastic fields using 

an exponential variogram model. The correlation length of 13,123 ft (4,000 m), the same as that used in 

the regional model, was used to create the perturbed multiplier fields, although in this case (Model A) 

hydraulic conductivity fields, and not multiplier fields, were created. A log-based sill of 0.1 and a mean 

value of 6,647 ft/d (2,026 m/d) were used to generate the stochastic hydraulic conductivity fields 

evaluated with conceptual Models A and B (these parameters are the same as those used in the regional

model). Unlike the regional model, however, the stochastic hydraulic conductivity fields in the models 

presented herein were not conditioned using APT point data, as these models are conceptual. The 

generated hydraulic conductivity realizations were used, along with the values of effective porosity from 

the regional model, to create 10,000 versions of Model A that were statistically similar to the fields 

evaluated with the regional model. 

For Model B, calculations were made to create horizontal hydraulic conductivity realizations for 

layers 1, 2, and 3 from the horizontal hydraulic conductivity realizations used in Model A. The bulk 

horizontal hydraulic conductivity used in Model A was assumed to be the thickness-weighted horizontal 

hydraulic conductivity of the individual layers of Model B. It was assumed that the horizontal hydraulic 

conductivity of the flow-zone layers 1 and 3 were equal (K1 = K3) and that the horizontal hydraulic 

conductivity of layer 2 was two orders of magnitude smaller than in layers 1 and 3 (100*K2=K1=K3). 

Model B layers 1 and 3 are isotropic and the vertical hydraulic conductivity of model layer 2 is an order 

of magnitude greater than the horizontal hydraulic conductivity.

Effective porosity values used in the regional model were applied with spatial uniformity in 

Model A runs. Based on field observations (Cunningham and other, 2004a, b, 2006a, b, 2009; Renken 

and others, 2008; Shapiro et al, 2008; Harvey and others, 2008), flow zones in the study area (layers 1 

and 3) may rapidly move groundwater and potential contaminants. Therefore, the effective porosity of 

the flow zones were calculated, using their thicknesses, from the effective porosity values used in the 

one-layer Model A runs. It was assumed that the effective porosities of layer 1 and layer 3 were equal

because they are both flow zones. The effective porosity of matrix layer 2 was set at 0.176. The 

following equation was solved for rB

ne ,3, (which is equivalent to rB

ne ,1, ) for each realization:
45

∑=

+ =

3,2,1

,

3,1, )( ,3,

,1,

l

lB

BBe

e b

bbn

n

rB

rA (1.1) 

where 

rlK

ne ,, is the effective porosity of layer l in realization r of Model K, and lK b , is the 

thickness of layer l in Model K, where K is either A or B. 

In the Monte Carlo simulations, Models A and B were further modified to include two canals, 

which were simulated using the MODFLOW River (RIV) Package. The locations of these two canals 

were the same in Models A and B, except in Model B the canals are located only in layer 1. The stages, 

conductances and bottom altitudes of the canals were modeled to resemble the L-30 Canal (to the west 

of the NWWF) and the canal directly to the north of the NWWF. The stage of the western canal was set 

at 3.12 ft (0.95 m) and the stage of the northern canal was set at 1.31 ft (0.40 m) based on data collected 

during dry conditions on June 19, 2004. The locations and shapes of the canals were meant to mimic the 

general location and shape of these canals with respect to the well fields, not the actual shape of the 

canals themselves. These canals in the regional model effectively enclose the well field (along with 

other connecting canals) whereas in the conceptual models simulated herein, the canals are implemented 

to investigate and compare the hydraulic effects of strong and weak canals, in general, in both a one 

layer and multi-layer context, on capture zones and drawdowns of the well field. The modified version 

of MODFLOW-2005, which uses the horizontal hydraulic conductivity of the cell to compute the RIV 

conductance, was used in Model A and B simulations; this allowed canal conductances to be varied 

stochastically as well. The weak sink option was used in MODPATH 5.0 to allow particles to track 

backwards through RIV cells that were not strong sources of water. 

Lakes were not added to Models A and B Monte Carlo simulations because the lakes range in 

depth from 60 to 80 ft. As a result, the lakes would fully penetrate all model layers in Models A and B 

and would not result in any noticeable difference in the results. 

Constant heads (fig. 1-2) were imposed at the left and right boundaries to allow an ambient flow 

field similar to that present in the regional model during dry conditions on June 19, 2004. The western 

boundary head was set at 2.95 ft (0.9 m) and the eastern boundary head was set at 0.98 ft (0.3 m). For 

Model B runs, the constant heads were set only in layer 1. 

The wells (fig. 1-2) were placed in a line of 6 cells and the 80 Mgal/d pumpage was distributed 

evenly between these cells to reflect the configurations of the actual well fields of the NWWF and 

WWF. The wells in the Model B runs were placed in layer 3. Particles were placed around the wells in 

layer 1 in Model A runs and around the wells in layer 3 in Model B runs. Particles were placed in 10 

different vertical locations throughout layer 3 in the Model B runs. To ensure particle endpoints in 

Models A and B would be comparable to one another, particles were vertically spaced throughout the 

bottom 34 percent of the Model A runs, because this is equivalent to the depth of the layer 3 production 

zone in Model B runs. The particles were tracked backwards for 210 days to obtain the 210-day capture 

zone. 
46

Initially, a control run was completed without well-field withdrawals to obtain a distribution of 

starting heads with which to calculate drawdown. Next a run was completed with well-field withdrawals 

of 80 Mgal/d. 

After 210 days of backward tracking from the wells, the particle endpoints were statistically 

analyzed using the technique adopted from Varljen and Shafer (1991), described in this report. The 

median and 95-percent confidence intervals (CIs) for the 210-day capture zone and for 1-ft drawdown 

were calculated. 

In some simulations for Model B, cells in layer 1 became dry. These simulations were excluded 

in both the Model A and Model B run results so that results could be compared between models in a 

consistent manner. Following the exclusion process, results from 7,662 runs were available for 

comparison. Results show that 1-ft drawdown and 210-day capture zones were comparable but slightly 

more extensive for 1-layer Model A than for 3-layer Model B (fig. 1-3 and fig. 1-4). Results from this 

investigation are consistent with findings of Sepulveda and Kuniansky (2010) in which it was 

determined that flow in the unconfined surficial aquifer system in Florida could be simulated reasonably 

well as a one-layer unit with calculated composite horizontal hydraulic conductivity based on a 

thickness-weighted arithmetic mean of multiple lithologic units. 

Reilly and Pollock (1993) define the contributing area of a well as the surface area that 

contributes recharge to the well. This differs from a capture zone, which is defined by the volume of 

aquifer water that enters the well within a specified amount of time. Forward and backward tracking 

will not always yield the same contributing area and forward tracking is preferred for contributing area 

determination, while backward tracking is often used for delineation of capture zone (Barlow, 1997).

Forward tracking with the simple model was accomplished with particles placed on the tops of cells for 

each cell in the area surrounding the wells and canals in the one-layer model (not shown). These 

forward-tracked particles all ended within the well cells and their pathlines were consistent with 

pathlines of backward-tracked particles from the wells. Particles to the west of the western canal and to 

the north of the northern canal moved past the canals to the wells. Thus, the test simulations indicate 

that with a weak sink option chosen with MODPATH, particles will move past boundary conditions 

such as canals towards the wells with the one-layer model. 
47

Appendix 2. Modification of MODFLOW-2005 Code 

Conductances of canals simulated using the MODFLOW-2005 River (RIV) Package were 

stochastic parameters in this investigation and were unique and variable along each canal for each of the 

10,000 realizations. MODFLOW-2005 was modified and recompiled to calculate canal conductance 

internally from the assigned horizontal hydraulic conductivity of the aquifer cell in which a canal is 

located (found in the Layer-Property Flow (LPF) Package in the model), and the assigned conductance 

value of each RIV cell, which was specified in the RIV Package. This conductance value was set for 

each realization and equal to L*W/M (see Equation 1). 

 Within the module GWF2RIV7.f, the subroutine GWF2RIV7FM was changed to calculate 

CRIV from HK (from the LPF Package) and RIVR (from the RIV Package):

 SUBROUTINE GWF2RIV7FM(IGRID)

C ******************************************************************

C ADD RIVER TERMS TO RHS AND HCOF

C ******************************************************************

C

C SPECIFICATIONS:

C ------------------------------------------------------------------

 USE GLOBAL, ONLY:IBOUND,HNEW,RHS,HCOF

 USE GWFRIVMODULE, ONLY:NRIVER,RIVR

C SEE LINZY

 USE GWFLPFMODULE, ONLY:HK

C ------------------------------------------------------------------

 CALL SGWF2RIV7PNT(IGRID)

C

C1------IF NRIVER<=0 THERE ARE NO RIVERS. RETURN.

 IF(NRIVER.LE.0)RETURN

C

C2------PROCESS EACH CELL IN THE RIVER LIST.

 DO 100 L=1,NRIVER

C

C3------GET COLUMN, ROW, AND LAYER OF CELL CONTAINING REACH.

 IL=RIVR(1,L)

 IR=RIVR(2,L)

 IC=RIVR(3,L)

C

C4------IF THE CELL IS EXTERNAL SKIP IT.

 IF(IBOUND(IC,IR,IL).LE.0)GO TO 100

C

C5------SINCE THE CELL IS INTERNAL GET THE RIVER DATA.

 HRIV=RIVR(4,L)

C COND CHANGED HERE TO MULTIPLY BY HK

!<CDL> CHANGE MADE BY LINZY TO COMPUTE CANAL CONDUCTANCE BASED ON HK

!O CRIV=RIVR(5,L)

 CRIV=RIVR(5,L)*HK(IC,IR,IL)

 ! OUTPUT CONDUCTANCE VALUE

 !WRITE(75,'(3I10,E10.3)') IL,IR,IC,CRIV

!</CDL>

 RBOT=RIVR(6,L)

 RRBOT=RBOT

C

C6------COMPARE AQUIFER HEAD TO BOTTOM OF STREAM BED.

 IF(HNEW(IC,IR,IL).LE.RRBOT)GO TO 96

C

C7------SINCE HEAD>BOTTOM ADD TERMS TO RHS AND HCOF.
48

 RHS(IC,IR,IL)=RHS(IC,IR,IL)-CRIV*HRIV

 HCOF(IC,IR,IL)=HCOF(IC,IR,IL)-CRIV

 GO TO 100

C

C8------SINCE HEAD<BOTTOM ADD TERM ONLY TO RHS.

 96 RHS(IC,IR,IL)=RHS(IC,IR,IL)-CRIV*(HRIV-RBOT)

 100 CONTINUE

 

C

C9------RETURN

 RETURN

 END

Also, the subroutine GWF2RIV7BD is changed to output the calculated CRIV to a file fort.75: 

 SUBROUTINE GWF2RIV7BD(KSTP,KPER,IGRID)

C ******************************************************************

C CALCULATE VOLUMETRIC BUDGET FOR RIVERS

C ******************************************************************

C

C SPECIFICATIONS:

C ------------------------------------------------------------------

 USE GLOBAL, ONLY:IOUT,NCOL,NROW,NLAY,IBOUND,HNEW,BUFF

 USE GWFBASMODULE,ONLY:MSUM,ICBCFL,IAUXSV,DELT,PERTIM,TOTIM,

 1 VBVL,VBNM

 USE GWFRIVMODULE,ONLY:NRIVER,IRIVCB,RIVR,NRIVVL,RIVAUX

C SEE LINZY

 USE GWFLPFMODULE, ONLY:HK

C

 DOUBLE PRECISION HHNEW,CHRIV,RRBOT,CCRIV,RATIN,RATOUT,RRATE

 CHARACTER*16 TEXT

 DATA TEXT /' RIVER LEAKAGE'/ 

C ------------------------------------------------------------------

 CALL SGWF2RIV7PNT(IGRID)

C

C1------INITIALIZE CELL-BY-CELL FLOW TERM FLAG (IBD) AND

C1------ACCUMULATORS (RATIN AND RATOUT).

 ZERO=0.

 RATIN=ZERO

 RATOUT=ZERO

 IBD=0

 IF(IRIVCB.LT.0 .AND. ICBCFL.NE.0) IBD=-1 

 IF(IRIVCB.GT.0) IBD=ICBCFL

 IBDLBL=0

C

C2------IF CELL-BY-CELL FLOWS WILL BE SAVED AS A LIST, WRITE HEADER.

 IF(IBD.EQ.2) THEN

 NAUX=NRIVVL-7 

 IF(IAUXSV.EQ.0) NAUX=0

 CALL UBDSV4(KSTP,KPER,TEXT,NAUX,RIVAUX,IRIVCB,NCOL,NROW,NLAY,

 1 NRIVER,IOUT,DELT,PERTIM,TOTIM,IBOUND)

 END IF

C

C3------CLEAR THE BUFFER.

 DO 50 IL=1,NLAY

 DO 50 IR=1,NROW

 DO 50 IC=1,NCOL

 BUFF(IC,IR,IL)=ZERO
49

50 CONTINUE

C

C4------IF NO REACHES, SKIP FLOW CALCULATIONS.

 IF(NRIVER.EQ.0)GO TO 200

C

C5------LOOP THROUGH EACH RIVER REACH CALCULATING FLOW.

 DO 100 L=1,NRIVER

C

C5A-----GET LAYER, ROW & COLUMN OF CELL CONTAINING REACH.

 IL=RIVR(1,L)

 IR=RIVR(2,L)

 IC=RIVR(3,L)

 RATE=ZERO

C

C5B-----IF CELL IS NO-FLOW OR CONSTANT-HEAD MOVE ON TO NEXT REACH.

 IF(IBOUND(IC,IR,IL).LE.0)GO TO 99

C

C5C-----GET RIVER PARAMETERS FROM RIVER LIST.

 HRIV=RIVR(4,L)

C COND CHANGED HERE TO MULTIPLY BY HK

!<CDL>

!O CRIV=RIVR(5,L)

 CRIV=RIVR(5,L)*HK(IC,IR,IL)

 WRITE(75,'(3I10,E10.3)') IL,IR,IC,CRIV

!</CDL>

 RBOT=RIVR(6,L)

 RRBOT=RBOT

 HHNEW=HNEW(IC,IR,IL)

C

C5D-----COMPARE HEAD IN AQUIFER TO BOTTOM OF RIVERBED.

 IF(HHNEW.GT.RRBOT) THEN

C

C5E-----AQUIFER HEAD > BOTTOM THEN RATE=CRIV*(HRIV-HNEW).

 CCRIV=CRIV

 CHRIV=CRIV*HRIV

 RRATE=CHRIV - CCRIV*HHNEW

 RATE=RRATE

C

C5F-----AQUIFER HEAD < BOTTOM THEN RATE=CRIV*(HRIV-RBOT).

 ELSE

 RATE=CRIV*(HRIV-RBOT)

 RRATE=RATE

 END IF

C

C5G-----PRINT THE INDIVIDUAL RATES IF REQUESTED(IRIVCB<0).

 IF(IBD.LT.0) THEN

 IF(IBDLBL.EQ.0) WRITE(IOUT,61) TEXT,KPER,KSTP

 61 FORMAT(1X,/1X,A,' PERIOD ',I4,' STEP ',I3)

 WRITE(IOUT,62) L,IL,IR,IC,RATE

 62 FORMAT(1X,'REACH ',I6,' LAYER ',I3,' ROW ',I5,' COL ',I5,

 1 ' RATE',1PG15.6)

 IBDLBL=1

 END IF

C

C5H------ADD RATE TO BUFFER.

 BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+RATE

C

C5I-----SEE IF FLOW IS INTO AQUIFER OR INTO RIVER.
50

 IF(RATE.LT.ZERO) THEN

C

C5J-----AQUIFER IS DISCHARGING TO RIVER SUBTRACT RATE FROM RATOUT.

 RATOUT=RATOUT-RRATE

 ELSE

C

C5K-----AQUIFER IS RECHARGED FROM RIVER; ADD RATE TO RATIN.

 RATIN=RATIN+RRATE

 END IF

C

C5L-----IF SAVING CELL-BY-CELL FLOWS IN A LIST, WRITE FLOW. ALSO

C5L-----COPY FLOW TO RIVR.

 99 IF(IBD.EQ.2) CALL UBDSVB(IRIVCB,NCOL,NROW,IC,IR,IL,RATE,

 1 RIVR(:,L),NRIVVL,NAUX,7,IBOUND,NLAY)

 RIVR(NRIVVL,L)=RATE

 100 CONTINUE

C

C6------IF CELL-BY-CELL FLOW WILL BE SAVED AS A 3-D ARRAY,

C6------CALL UBUDSV TO SAVE THEM.

 IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IRIVCB,BUFF,NCOL,NROW,

 1 NLAY,IOUT)

C

C7------MOVE RATES,VOLUMES & LABELS INTO ARRAYS FOR PRINTING.

 200 RIN=RATIN

 ROUT=RATOUT

 VBVL(3,MSUM)=RIN

 VBVL(4,MSUM)=ROUT

 VBVL(1,MSUM)=VBVL(1,MSUM)+RIN*DELT

 VBVL(2,MSUM)=VBVL(2,MSUM)+ROUT*DELT

 VBNM(MSUM)=TEXT

C

C8------INCREMENT BUDGET TERM COUNTER.

 MSUM=MSUM+1

C

C9------RETURN.

 RETURN

 END
51

Appendix 3. Test of Boundary Condition near the NWWF 

The western boundary canals are hydraulically connected to the Biscayne aquifer and control the 

stage near the well fields. The aquifer does extend west of the L-30 and L-31 Canals, however, and 

underflow has been documented. A simple test with the base-case, steady-state, average condition was 

made to check that the western canal does indeed control the head near the NWWF, and that simulated 

water levels near the well field would be identical if another boundary condition was applied. The test 

was accomplished by making the aquifer system active west of the L-30 Canal. Additionally, a general 

head boundary (GHB) was set at the western edge of the model using the average head at two gages for 

December 19, 2002, the day chosen to represent average hydrologic conditions for the steady-state 

model (Station “USGS 254754080344300 SHARK RIVER SLOUGH NO.1 IN CONS.3B NR 

COOPERTOWN” with daily stage of 7.74 ft and Station “USGS 255250080335001 SITE 71 IN 

CONSERVATION AREA 3B NR COOPERTOWN, FL” with daily stage of 8.07 ft, which averages to 

7.91 ft NGVD29). December is at the end of the wet season and much of the area west of L-30 Canal is 

inundated with water. Stage in the L-30 Canal is 4.79 ft (1.46 m), thus there is an eastward gradient in 

the Biscayne aquifer from the conservation area and the Everglades. Two tests were conducted with the 

general-head conductance set at 269 square feet per day (ft2

/d) (25 square meters per day (m2/d)) (cell 

area 26,910 square feet (ft2

) (2500 square meters (m2)) multiplied by a vertical hydraulic conductivity 

(Kv) of 0.0328 ft/d (0.01 m/d) and a 3.28 ft (1 m) thickness of sediment) and 26,910 ft2

/d (2500 m2/d) 

(cell area 26,910 ft2

 (2500 m2) multiplied by a Kv of 3.28 ft/d (1 m/d) and a 3.28 ft (1 m) thickness of 

sediment). Water-level contours remain the same near the NWWF (fig. 3-1). The change in conductance 

of the GHB only affects water levels between the levee and the western extent of the new active area. 

The test in Appendix 1 indicated that particles can move beneath a canal when the weak sink 

option is used in MODPATH. The boundary condition tests in this appendix indicate that the L-30 

Canal provides a strong hydraulic head boundary condition and pathlines from the canals to the well 

fields are the same if the area west of the L-30 Canal is actively simulated or not actively simulated. 

The area west of the L-30 Canal is not urban and the goal of the study was to estimate capture zones 

within urban areas for a maximum of 210 days of particle travel. In order to save on computation times, 

this area was made inactive and thus the simulations cannot be used to determine which parts of the 

natural wetlands are in the overall contributing area of the NWWF. 
52

Appendix 4. Comparative analysis of lake representation using high 

hydraulic conductivity zones and the LAK Package 

In groundwater models, lakes are commonly represented by cells with high hydraulic 

conductivity (Lee, 1996; Anderson and others, 2002; Hunt and others, 2003). The use of high hydraulic 

conductivity cells as an approach to represent lakes was compared to the explicit representation of lakes 

using the MODFLOW Lake (LAK) Package (Merritt and Konikow, 2000) to support the approach taken 

for this investigation. 

When the LAK Package is used, model cells delineated as lake cells are excluded from 

groundwater flow calculations and the relation between stage, volume, and surface area for individual 

lakes are determined directly from the model layer geometry. Individual lakes can be defined as single 

grid cells or contiguous groups of grid cells. Aquifer-lake exchanges are a function of the difference 

between the simulated lake stage and simulated head in active groundwater cells immediately adjacent 

(horizontally and vertically) to the defined lake perimeter. Lake conductances are calculated as the 

harmonic mean of lake sediment and aquifer hydraulic properties. The LAK Package includes 

functionality to simulate lakes with varying lateral extents, as well as coalescing lakes; available 

options, functionality, and numerical implementation of the LAK Package are detailed in Merritt and 

Konikow (2000). 

Some of the advantages of using the LAK Package are (1) the ability to simulate a single stage 

for a lake, (2) the ability to easily calculate lake water budgets, and (3) the elimination of model 

convergence issues related to use of high hydraulic conductivities to represent lakes. The most 

substantial limitations of the LAK Package for this investigation are (1) the difficulty of achieving 

steady-state lake stages when initial lake stages and groundwater levels differ substantially from steadystate conditions and lake volumes are large, and (2) the termination of MODPATH 5.0 (Pollock, 1994) 

particles in lake cells, which limits the ability to determine well-field capture zones. 

A one-layer steady-state MODFLOW model was developed using parameters representative of 

the Biscayne aquifer in northwestern Miami-Dade County to evaluate the appropriateness of the high 

hydraulic conductivity approach to represent lakes surrounding the NWWF (fig. 4-1). A total of 38 

separate quarry lakes were simulated. The model domain comprised 300 columns and rows, with a

constant horizontal grid spacing of 164 ft (50 m). Uniform upper and lower model altitudes of 9.84 ft (3 

m) and -88.6 ft (-27 m) NAVD 88, respectively, were specified and resulted in an aquifer thickness of 

approximately 100 ft. A hydraulic conductivity of 7,136 ft/d (2,175 m/d) was assigned to all active 

model cells in both models. For the high hydraulic conductivity simulation, a hydraulic conductivity of 

7.136 × 105 ft/d (2.175 × 105 m/d) was assigned to the lake cells. 

A constant head of 9.51 ft (2.90 m) and 6.92 ft (2.11 m) was assigned to the left and right 

boundaries, respectively, in both models. A no-flow boundary was assigned to the top and bottom 

boundaries in both models; recharge was not applied in either model. A maximum evapotranspiration 

rate of 32.5 in/yr was applied to active groundwater and lake cells in both models. An 

evapotranspiration extinction depth of 3.28 ft (1 m) was applied in all active groundwater cells. An 

evapotranspiration surface altitude of 9.84 ft (3 m) was applied in cells representing the aquifer and an 

evaporation surface altitude of -65.62 ft (-20 m) was applied to quarry lakes represented using the high 

hydraulic conductivity approach. Fifteen production wells, each having a constant individual 

groundwater withdrawal rate of 15 Mgal/d for a total well-field withdrawal rate of 225 Mgal/d, were 

simulated in both models. Canals were not simulated in either model. For the LAK Package model, a 

lake-bed leakance value of 2.175 × 105

 feet per day per foot ((ft/d)/ft) was assigned to each lake cell and 

is based on a sediment hydraulic conductivity of 7,136 ft/d (2,175 m/d) and a sediment thickness of 
53

0.328 ft (0.1 m). Backward-tracked particle pathlines from the 15 production wells were calculated 

using MODPATH 5.0 (Pollock, 1994) for both models. For the LAK Package simulation, particle 

locations on the upgradient side of lakes were determined for the simulations in which the lakes were 

represented with high hydraulic conductivity. These particle locations on the upgradient side of lakes 

were used to define additional starting particle locations for the LAK Package simulations, and the 

resulting pathlines were directly compared with particle pathlines for the simulations representing lakes 

with high hydraulic conductivity throughout the model domain. 

The geometric multigrid (GMG) solver (Wilson and Naff, 2004) was used in both models. 

Residual (l2-norm) and head convergence criteria of 1,766 cubic feet per day (ft3/d) (50 cubic meters 

per day (m3/d)) and 3.28 × 10-6 ft (1 × 10-6 m) were used to ensure that lake stages and volumes were

not changing in the LAK Package model. A head convergence criterion greater than 3.28 × 10-6 ft (1 × 

10-6 m) resulted in lake levels that would change if final groundwater heads and lake levels from one 

simulation were used as initial conditions for a second simulation. In all cases, the final l2

-norm of the 

residual and maximum head change for both models was less than 272 ft3/d (7.7 m3/d) and 3.28 × 10-6 ft 

(1 × 10-6 m), respectively. 

Simulated steady-state groundwater heads and particle pathlines for the high hydraulic 

conductivity model and LAK Package model are shown in figure 4-1. Quantitative comparison of 

simulated groundwater inflows and outflows and evapotranspiration from the 38 defined quarry lakes

are summarized in table 4-1. All lake water-budget flow terms simulated using the high hydraulic 

conductivity approach and the LAK Package differed by less than ±2.0 percent. Comparison of 

groundwater heads, particle pathlines, and quantitative lake water-budget terms indicates the high 

hydraulic conductivity approach for representing quarry lakes adjacent to the NWWF yields similar 

results to the approach in which lakes are explicitly represented using the LAK Package. Thus, the 

approach using high hydraulic conductivity cells for lakes, which is easier and introduces fewer 

parameters, was used for this investigation. 
54

Appendix 5. Additional Groundwater Head Calibration Plots 

Figure 5-1. Daily observed and model-simulated water levels for monitoring wells in the study area not 

shown in figure 13.

Figure 5-1 is shown on page .

Appendix 6. Comparison with an Alternative Particle-Tracking Analysis 

Approach 

Chin and others (2010) proposed a method to estimate well field travel-time capture zones in 

areas characterized by lakes. The method uses an advective-dispersive backward particle-tracking 

algorithm that accommodates lake residence times to determine travel-time contours. Lake residence 
55

times of particles that traverse lakes are added to the time required for the particle to reach the lake from 

the pumping well in order to delineate the backward-tracked travel-time contour. If the sum of the lake 

residence time and travel time from the pumping well to the lake is greater than the travel-time contour 

of interest, then the entire lake should be included in the capture zone. Conversely, if the sum of the lake 

residence time and travel time from the pumping well to the lake is less than the capture zone of interest, 

the difference between the capture zone of interest and the sum is calculated, and particles are tracked 

backwards from the inflow side of the lake for this additional length of time. 

This method is compared to an alternative method of tracking particles backwards through lakes 

with no residence-time adjustment (fig. 6-1). The particle tracks backwards to the lake in 24 days (fig. 

6-1A). The particle continues to move through the lake and reaches the other side, where it exits, enters 

a second lake at 70 days, and exits that lake at 80 days to continue backward tracking. In contrast, figure 

6-1B shows results with the residence-time post-processing adjustment to the particle times. This 

particle tracks backwards to the first lake, the residence time (10 days) is added to the travel time from 

the pumping well to the lake (24 days), and the particle location is adjusted to the exit point of the lake 

(actually the entry point because this is backward tracking), for a total travel time thus far of 34 days. 

By comparison, the particle in figure 6-1A does not exit the lake until 57 days of travel. The particle in 

figure 6-1B continues to track backwards for 13 more days, encountering the second lake, with a 

residence time of 5 days, where its location is adjusted to the exit point of the lake where the travel time

is 52 days. By comparison, the particle location in figure 6-1A had a travel time of 80 days. 

Chin and others (2010) illustrated this residence-time backward particle-tracking approach with 

an advective-dispersive code for particle tracking for the NWWF using a previously published multilayer calibrated model in which residence times for lakes near the pumping wells were estimated by 

dividing the volume of the lake by the steady-state outflow from the lake into groundwater (Chin and 

others, 2010). This calculation of residence time is based on the assumption that the lake is well-mixed. 

To compare results of the model used in the current study with the results presented by Chin and 

others (2010), a modified version of the deterministic base-case simulation representing dry hydrologic 

conditions (June 19, 2004) was created. Parameters and details were changed in the dry hydrologic 

conditions steady-state model to better match details in the model of Chin and others (2010). The 

maximum pumping rate for the NWWF was changed to 158 Mgal/d. Lakes were specified by Chin and 

others (2010) using a hydraulic conductivity of 1.00 × 106

 ft/d (305,000 m/d), a value greater than that 

used in base-case and Monte Carlo simulations evaluated in the current study. Both hydraulic 

conductivities were tested and calculated residence times were not found to be sensitive to this 

parameter, likely because both values for lake conductivity are very large. The lakes highlighted in Chin 

and others (2010), namely, 22, 23, 28, 30, 32 and 33, were simulated for steady-state dry hydrologic 

conditions using the lake coverage of Chin instead of the 2004 lake coverage (fig. 6-2 and fig. 12).

Residence times were estimated using lake outflow rates calculated from steady-state groundwater 

water-budget terms for individual lakes by multiplying the percentage of lake transmissivity to total cell 

transmissivity by the value for that cell in the cell-by-cell flow file. This scaling created an adjusted cellby-cell flow file. This adjusted cell-by-cell flow file, which is scaled to better represent flows from lakes 

into the groundwater as opposed to the total flow between cells, was input into ZONEBUDGET 

(Harbaugh, 1990) to calculate lake outflows into the groundwater for each lake. Lake volumes were 

calculated by multiplying lake cell areas and the depths of the lakes in each cell, which were estimated 

from bathymetric data and rock-mining permits (J. Baker, Miami-Dade Department of Permitting, 

Environment, and Regulatory Affairs, unpub. data, 2011). Where lake-bottom altitude data were 

unavailable, a depth of 60 ft was assumed based on the majority of rock-mining permits for the Lake 

Belt region.
56

A quantitative comparison between results of this investigation and Chin and others (2010) was 

made for the lakes shown in figure 6-2. Residence times for the model used in this investigation are 

mostly smaller than those calculated by Chin and others (2010) (table 6-1), likely because the lake 

volumes are generally smaller than those used by Chin and others (2010) (lake-bottom altitudes were 

not reported in Chin and others, 2010). Furthermore, flow rates from the lake into the aquifer were 

generally larger in this investigation (except for lakes 30 and 33) than those in Chin and others (2010),

likely because dry conditions represented in this investigation differ from the conditions used in the 

model of Chin and others (2010). Smaller residence times indicate that particles will remain in the lakes

in this investigation for a shorter duration than was calculated by Chin and others (2010). The one 

exception is lake 30, located southwest of the NWWF, for which the residence time calculated by the 

current model is approximately twice the value calculated by Chin and others (2010). Smaller residence 

times imply more particles will potentially need to be backward-tracked from the inflow side of the lake 

and potentially for a longer amount of time than those particles in the model of Chin and others (2010). 

This finding could lead to, in general, more spatially extensive capture zone estimates using the model 

in this investigation. The residence-time post-processing particle adjustment used in this investigation 

addresses the need to track particles backward beyond lake inflows to estimate capture zones, or to 

include the entire lake within the capture zone, depending on the duration of the capture zone and lake 

residence time, as recommended by Chin and others (2010). 

Model grid extent

Base map 1

Base map 2

Base map 3

Base map 4

Base map 5

Base map 6

Study area

0 5 10 MILES EXPLANATION

0 5 10 KILOMETERS

Water 

Conservation

41 Area 3B

27 75

195

1

826

826

874

836

997

1

997

COLLIER COUNTY

BROWARD COUNTY

MIAMI-DADE

COUNTY

MONROE COUNTY

EVERGLADES

NATIONAL

PARK

Homestead

Hialeah

Miami

25°45'

25°30'

25°15'

26°00'

80°45' 80°30' 80°15'

Figure 1. Locations of all base maps used throughout this report. Base maps are referred to in each 

map figure caption.

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

North 

Miami

Cutler

Ridge

L-67A Canal

L-67C Canal

Mowry Canal

Tamiami Canal

L-29 Canal

Florida Bay

Biscayne Bay

ATLANTIC OCEAN

L-67 Canal

C-111 Canal

Card

Sound

Barnes

Sound

L-30 Canal

FLORIDA

ATLANTIC OCEAN

GULF OF MEXICO

Lake

Okeechobee

Base

map 1

58
0 5 10 MILES

0 5 10 KILOMETERS

997

1 27 75

195

1

826

826

874

836

997

EXPLANATION

Line of section

Model grid extent

Study area

Lake Belt area

Production well

Biscayne National 

 Park boundary

A A'

A A'

B' B

BROWARD COUNTY

MIAMI-DADE

COUNTY

EVERGLADES

NATIONAL

PARK

25°45'

25°30'

25°15'

26°00'

80°45' 80°30' 80°15'

Figure 2. Model grid extent, study area, hydrogeologic sections A–A’ and B–B’ , location of production 

well fields and location of canals, southeastern Florida. Area shown is Base map 1 (fig. 1). NWWF is 

Northwest Well Field. WWF is West Well Field.

WWF

Water Conservation

Area 3B

NWWF

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

Homestead

Hialeah

Miami

North 

Miami

Cutler

Ridge

L-67A Canal

L-67C Canal

Mowry Canal

Tamiami Canal

L-29 Canal

Biscayne Bay

ATLANTIC OCEAN

L-67 Canal

C-111 Canal

Card

Sound

Barnes

Sound

L-30 Canal

59
EXPLANATION

Study area

Lake Belt area

Lakes

Production well

Wellhead Protection Areas for 

 production well fields

210-day travel-time capture zone

100-day travel-time capture zone

30-day travel-time capture zone

10-day travel-time capture zone

Outer protection boundary

0 2 4 MILES

0 2 4 KILOMETERS

Figure 3. Study area and existing Wellhead Protection Areas for production well fields in Miami-Dade County 

(Miami-Dade County Department of Regulatory and Economic Resources, 

http://www.miamidade.gov/derm/library/water/wellfield_protection_areas_poster.pdf, accessed 

February 8, 2012). Area shown is Base map 2 (fig. 1). NWWF is Northwest Well Field. WWF is West Well Field.

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EVERGLADES

NATIONAL

PARK

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-30 Canal Pennsuco Wetlands

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

60
Figure 4. Hydrogeologic cross section with subset section A–A’. Modified from Reese and 

Cunningham (2000). Location of subset cross section shown in fig. 2.

VERTICAL SCALE GREATLY EXAGGERATED

Pamlico Sand

Pinecrest Sand

Fort Thompson

Formation

Miami Limestone

Key Largo Limestone

Tamiami

Formation

Semiconfining

Unit

CSU

Water-table Aquifer

ICU Peace River Formation

Ochopee Limestone

Ochopee

Limestone

SA

Gray Limestone

Aquifer

A A'

Biscayne

Aquifer

Biscayne

Aquifer

Oolitic limestone Hydrostratigraphic boundary

Lithostratigraphic boundary

Confining to semiconfining unit

Intermediate confining unit

Sand aquifer

Formation

Reef rock and marine carbonate

Quartz sand or sandstone

Silt

Terrigenous mudstone

Pelecypod rudstone and floatstone;

 minor calcrete beds and laminations;

 locally planorbella floatstone

Pelecypod rudstone and floatstone;

 minor pelecypod-rich quartz sand and

 sandstone

CSU

ICU

SA

Fm

EXPLANATION

NGVD 29 NGVD 29

Altitude, in feetAltitude, in feet

–100 –100

–50 –50

–200 –200

–150 –150

Miami limestone

WEST EAST

Unnamed

Formation

Thin intercalations

of Holocene Lake Flirt Marl

and peat or peat only partly

cover Miami Limestone

Anastasia

Fm

Anastasia

Fm

0 5 10 MILES

0 5 10 KILOMETERS

61
West East

? DIST

No

DBI

Data

G-3773 G-3772

No

DBI

Data

VLS 2h 3a 3b 2g3 2e2 2d 2c2 4 5 2g22g1 2b 2a

G-3257A

No

DBI

Data

G-3816

No

DBI

Data

G-3817

No

DBI

Data

G-3258A

No

DBI

Data

VLS 2h 3aHFC HFC 5 HFC 4

3b 2g3 2e2 2d 5 4

BISCAYNE AQUIFER HU

FMQU

SEMICONFINING UNIT

UPPER CLASTIC UNIT OF TAMIAMI FORMATION FORT THOMPSON FORMATION MIAMI LIMESTONE

HFC 2

Q5 Q4

HFC 3

Q3 Q2 Q1?

2c2 2b

S-3163

2g2 2g1 2a

S-3164

No

DBI

Data

No

DBI

Data

NGVD 29 NGVD 29

4

–4

–8

–12

–16

–20

–24

–28

4

–4

–8

–12

–16

–20

–24

–28

 ?′ ANCE BETWEEN WELLS, IN METERS

9.99 9.99 9.99 9.99 9.99 9.99 9.99

5,012.7 1,869.6 32.75 263.52 6.00 2.19 1.31 62.44 34.61 2,680.4

NGVD 29 NGVD 29

Figure 5. Hydrogeologic cross section B–B’ across central part of the study area, including the 

Northwest Well Field (modified from Renken and others, 2008). Location of cross section shown 

in fig. 2.

EXPLANATION

Distance between wells, in feet

31,754

G-3795

G-3770

S-3163

G-3257A G-3258A G-3773 G-3816 G-3817 G-3772

S-3164

G-3792

SEMICONFINING UNIT BISCAYNE AQUIFER

Miami Limestone FORT THOMPSON FORMATION

UPPER CLASTIC UNIT OF 

TAMIAMI FORMATION

WEST EAST

B B'

No

DBI

Data

No

DBI

Data

No

DBI

Data

No

DBI

Data

No

DBI

Data

No

DBI

Data

No

DBI

Data

No

DBI

Data

No

DBI

Data

No

DBI

Data

No

DBI

Data

No

DBI

Data

No

DBI

Data

No

DBI

Data

HU

FM

QU

Production well— Brown fill 

 corresponds to cased section

 of well and no fill corresponds 

 to open-hole interval

Cased section of borehole

Digital borehole image

Megaporosity zone

Geologic formation

Hydrogeologic unit

Q-unit (Perkins, 1977)

Cycle boundary—Dashed, 

 where uncertain

Matrix porosity zones

DBI

HU

FM

QU

Q5

Q4

Q3

Q2

Q1

?

–10

–20

–30

–40

–50

–60

–70

–80

–90

10

–10

–20

–30

–40

–50

–60

–70

–80

–90

10

Altitude, in feet Altitude, in feet 

62
Mowry Canal

Tamiami Canal

L-29 Canal

Florida Bay

Biscayne Bay

ATLANTIC OCEAN

L-67 Canal

C-111 Canal

Card

Sound

Barnes

Sound

L-30 Canal

L-67A Canal

L-67C Canal

Water 

Conservation

Area 3B

APPROXIMATE

WESTERN

LIMIT

OF

BISCAYNE

AQUIFER

41

75

27

195

826

874

836

1

997

0 5 10 MILES

0 5 10 KILOMETERS

EXPLANATION

-20

Structure contour—Shows

 altitude of base of Biscayne

 aquifer, relative to NGVD 29. 

 Contour interval 10 feet. 

 Hachures indicate 

 depression 

Study area

-20

-30

-30

-20

-40

-40

-50

-50

-60

-60

-70

-70

-70

-80

-80

-80

-80

-170-180

-160 -150 -140

-130 -120

-90

-90

-90

-90

-100

-100

-110

-110

-100

COLLIER COUNTY

BROWARD COUNTY

MIAMI-DADE

COUNTY

MONROE COUNTY

EVERGLADES

NATIONAL

PARK

Homestead

Hialeah

Miami

25°45'

25°30'

25°15'

26°00'

80°45' 80°30' 80°15'

Figure 6. Altitude of the base of the Biscayne aquifer and its western extent throughout most of 

Miami-Dade County (Fish and Stewart, 1991).

North 

Miami

Cutler

Ridge

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

63
Figure 7. Daily municipal supply well withdrawals from the Biscayne aquifer in Miami-Dade County 

for January 1, 1996, through December 31, 2004. NWWF is Northwest Well Field. WWF is West Well Field.

0

50

100

150

200

250

300

350

400

450

1996 1997 1998 1999 2000 2001 2002 2003 2004

Year

Municipal groundwater withdrawals, 

in millions of gallons per day

EXPLANATION

County totals

NWWF/WWF totals

64
#

#

G119

S337

0 2 4 MILES

0 2 4 KILOMETERS

EXPLANATION

G119# Surface-water stage 

 recorder location

 and name

Study area

Lakes

Production well

Figure 8. Primary and secondary canal system in study area. Area shown is Base map 2 (fig. 1). NWWF is 

Northwest Well Field. WWF is West Well Field.

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

BROWARD COUNTY

MIAMI-DADE COUNTY

Water 

Conservation

Area 3B

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EVERGLADES

NATIONAL

PARK

L-31 N Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-30 Canal

L-29 Canal

Pennsuco Wetlands

65
EXPLANATION

Study area

Active cells (RCH, EVT, RIV)

Active cells (RCH, EVT)

Active cells (RCH, EVT, WEL)

Inactive cells

CHD cells

Model grid extent

Canals, not simulated

0 5 10 MILES

0 5 10 KILOMETERS

BROWARD 

COUNTY

MIAMI-DADE

COUNTY

EVERGLADES

NATIONAL

PARK

25°45'

25°30'

25°15'

26°00'

80°45' 80°30' 80°15'

Figure 9. Location of active and inactive model cells, and type of MODFLOW boundary condition. 

[RCH = Recharge Package, EVT = Evapotranspiration Package, RIV = River Package, 

CHD = Time-Variant Specified-Head Package, and WEL = Well Package]. Area shown is Base 

map 1 (fig. 1). NWWF is Northwest Well Field. WWF is West Well Field.

WWF

NWWF

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

Water 

Conservation

Area 3B

C-9 Canal

L-31 N Canal

L-67A Canal

L-67C Canal

Mowry Canal

Tamiami Canal

L-29 Canal

Biscayne Bay

ATLANTIC OCEAN

L-67 Canal

C-111 Canal

Card

Sound

Barnes

Sound

L-30 Canal

66
EXPLANATION

Model grid extent

Land surface altitude, 

 in feet NAVD 88

Less than 0

0 to 1.99

2 to 3.99

4 to 5.99

6 to 7.99

8 to 9.99

10 to 11.99

12 to 13.99

14 to 15.99

0 5 10 MILES

0 5 10 KILOMETERS L-31 N Canal

C-9 Canal

Mowry Canal

Biscayne Bay

ATLANTIC OCEAN

C-111 Canal

Card

Sound

Barnes

Sound

L-30 Canal

Figure 10. Altitude of the top of the Biscayne aquifer model. Area shown is Base map 1 (fig. 1).

80°40' 80°30' 80°20' 80°10'

26°00'

25°50'

25°40'

25°30'

25°20'

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

BROWARD 

COUNTY

MIAMI-DADE

COUNTY

EVERGLADES

NATIONAL

PARK

Water 

Conservation

Area 3B

L-67A Canal

L-67C Canal

Tamiami Canal

L-29 Canal L-67 Canal

67
Tamiami Canal

Levee 29

Levee 67

L-67A Canal

L-67C Canal

EXPLANATION

0 5 10 MILES

0 5 10 KILOMETERS

EVERGLADES

NATIONAL

PARK

Figure 11. Hydraulic conductivity field for the base-case simulation and locations of transmissivity 

estimates derived from aquifer performance tests. Area shown is Base map 1 (fig. 1).

WWF

NWWF

80°40' 80°30' 80°20' 80°10'

26°00'

25°50'

25°40'

25°30'

25°20'

BROWARD 

COUNTY

MIAMI-DADE

COUNTY

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

Water 

Conservation

Area 3B

C-9 Canal

L-31 N Canal

L-67A Canal

L-67C Canal

Mowry Canal

Tamiami Canal

L-29 Canal

Biscayne Bay

ATLANTIC OCEAN

L-67 Canal

C-111 Canal

Card

Sound

Barnes

Sound

L-30 Canal

Model grid extent

Production well

Simulated canals

Canals, not simulated

Transmissivity estimate

K, in feet per day

25,250

20,250

15,250

10,250

5,250

250

>30,000

68
0 2 4 MILES

0 2 4 KILOMETERS

EXPLANATION

Study area

Lakes (2004)

Production well

Lakes (1999)

Simulated canals

Canals, not simulated

L-31 N Canal

Snapper Creek Extension Canal

L-30 Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

Figure 12. Locations of lakes in 1999 (transient simulation) and 2004 (steady-state simulations). Area shown 

is Base map 2 (fig. 1). NWWF is Northwest Well Field. WWF is West Well Field.

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EVERGLADES

NATIONAL

PARK

69
1996 1997 1998 1999 2000 2001 2002 2003 2004

–3

–2

–1

0

1

2

3

4

5

6

7

Date

G-3253

1996 1997 1998 1999 2000 2001 2002 2003 2004

0

1

2

3

4

5

6

Date

G-3259A

1996 1997 1998 1999 2000 2001 2002 2003 2004

–8

–6

–4

–2

0

2

4

6

Date

G-551

1996 1997 1998 1999 2000 2001 2002 2003 2004

1

2

3

4

5

6

7

8

Date

G-855

1996 1997 1998 1999 2000 2001 2002 2003 2004

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

5.5

6.0

6.5

Date

G-975

1996 1997 1998 1999 2000 2001 2002 2003 2004

0

1

2

3

4

5

6

Date

G-976

Figure 13. Daily observed and model-simulated water levels for selected monitoring wells in the study area.

Head (NAVD 88), in feet Head (NAVD 88), in feet Head (NAVD 88), in feetHead (NAVD 88), in feet Head (NAVD 88), in feet Head (NAVD 88), in feet

EXPLANATION

Simulated groundwater heads

Observed groundwater heads

70
Krome

0.012

G-976

0.141

G-551

2.139

G-970

-0.029

G-1488

0.222

G-1637

0.439

G-973

-0.065

G-975

-0.139

G-3253

0.526

G-3574

0.605

G-3559

0.115

G-3551

0.026

G-855

-0.316 G-3561

-0.293

G-3555

-0.569

G-3560

-0.392

G-3259A

0.205

G-3572

-0.653

G-3439

-0.486

G-3558

-0.614

G-3556

-0.256

G-3554

-0.007

G-3552

-0.338

G-3264A

-0.038

G-3553

-0.017

0 2 4 MILES

0 2 4 KILOMETERS

EXPLANATION

Lakes (2004)

Production well

G-3561

-0.293

Groundwater well name 

 and mean error value 

 in feet

Simulated canals

Canals, not simulated

Figure 14. Locations of groundwater observation data used to calibrate model and mean error calculated 

from observed and simulated daily water levels. Area shown is Base map 2 (fig. 1). NWWF is Northwest 

Well Field. WWF is West Well Field.

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal

Snapper Creek Extension Canal

L-30 Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EVERGLADES

NATIONAL

PARK

71
0 2 4 MILES

0 2 4 KILOMETERS

EXPLANATION

Basins

Study area

Production well

Simulated canals

Canals, not simulated

Figure 15. Selected surface-water basins (shaded) in simulated area, locations of production wells, 

and canals. Area shown is Base map 3 (fig. 1). NWWF is Northwest Well Field. WWF is West Well Field.

BROWARD COUNTY

MIAMI-DADE COUNTY

80°30' 80°20'

26°00'

25°50'

25°40'

NWWF

WWF

C-4*

C-6

C-2*

C-3*

C-5*

Water 

Conservation

Area 3B

L-31 N Canal

L-30 Canal

L-67A Canal

L-67C Canal

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data EVERGLADES NATIONAL PARK

72
Figure 16. Monthly model-simulated (red) and observed (blue) net exchange for composite C-4 and C-6 basins, and 

mean error (ME), mean absolute error (MAE) and the coefficients of determination (R2

) for simulated versus observed 

net exchange for each basin.

–16

–12

–8

–4

0

4

8

12

16

Date

–16

–12

–8

–4

0

4

8

12

16

Date

C4 basin net exchange, 

in inches per month

C6 basin net exchange, 

in inches per month

Jan 96 Jul 96 Jan 97 Jul 97 Jan 98 Jul 98 Jan 99 Jul 99 Jan 00 Jul 00 Jan 01 Jul 01 Jan 02 Jul 02 Jan 03 Jul 03 Jan 04 Jul 04 Jan 05

Jan 96 Jul 96 Jan 97 Jul 97 Jan 98 Jul 98 Jan 99 Jul 99 Jan 00 Jul 00 Jan 01 Jul 01 Jan 02 Jul 02 Jan 03 Jul 03 Jan 04 Jul 04 Jan 05

ME = 0.62 inches per month

MAE = 2.21 inches per month

R2

 = 0.24

ME = 2.37 inches per month

MAE = 2.38 inches per month

R2

 = 0.78

Model-simulated

Observed

Canal recharge

to aquifer

Aquifer discharge 

to canal

Canal recharge

to aquifer

Aquifer discharge 

to canal

Model-simulated

Observed

73
Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

Water 

Conservation

Area 3B

C-9 Canal

L-67A Canal

L-67C Canal

Mowry Canal

Tamiami Canal

L-29 Canal

Biscayne Bay

ATLANTIC OCEAN

L-67 Canal

C-111 Canal

Card

Sound

Barnes

Sound

L-30 Canal

Figure 17. Average water-level differences for the Biscayne aquifer between May 1993 (observed) 

and May 1998 (simulated). Area shown is Base map 1 (fig. 1). White shaded area within model 

boundary comprises inactive areas of model. NWWF is Northwest Well Field. WWF is West 

Well Field.

EVERGLADES

NATIONAL

PARK

Study

area

BROWARD 

COUNTY

MIAMI-DADE

COUNTY

80°40' 80°30' 80°20' 80°10'

26°00'

25°50'

25°40'

25°30'

25°20'

0 5 10 MILES

0 5 10 KILOMETERS

WWF

NWWF

EXPLANATION

Water-level differences 

 (May), in feet

Less than –4

–3.9 to –3

–2.9 to –2

–1.9 to –1

–0.9 to 0

0.1 to 1

1.1 to 2

2.1 to 3

3.1 to 4

Model grid extent

Production well

Simulated canals

Canals, not simulated

Study area

74
0 5 10 MILES

0 5 10 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

Water 

Conservation

Area 3B

Mowry Canal

C-9 Canal

L-67A Canal

L-67C Canal

Mowry Canal

Tamiami Canal

L-29 Canal

Biscayne Bay

ATLANTIC OCEAN

L-67 Canal

C-111 Canal

Card

Sound

Barnes

Sound

L-30 Canal

Figure 18. Average water-level differences for the Biscayne aquifer between November 1993 

(observed) and November 1998 (simulated). Area shown is Base map 1 (fig. 1). White shaded 

area within model boundary comprises inactive areas of model. NWWF is Northwest Well Field. 

WWF is West Well Field.

80°40' 80°30' 80°20' 80°10'

26°00'

25°50'

25°40'

25°30'

25°20'

EVERGLADES

NATIONAL

PARK

Study

area

BROWARD 

COUNTY

MIAMI-DADE

COUNTY

WWF

NWWF

EXPLANATION

Water-level differences 

 (November), in feet

Less than –4

–3.9 to –3

–2.9 to –2

–1.9 to –1

–0.9 to 0

0.1 to 1

1.1 to 2

2.1 to 3

3.1 to 4

Model grid extent

Production well

Simulated canals

Canals, not simulated

Study area

75
EXPLANATION

0

0

5 MILES

5 KILOMETERS

2.5

2.5

Simulated canals

Canals, not simulated

K, in feet per day

Production well

Transmissivity data point

>30,000

25,250

20,250

15,250

10,250

5,250

0

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

Figure 19. Hydraulic conductivity field used in transient and deterministic base-case steady-state simulations in the study 

area for the (A) base-case simulation, and (B, C and D) for three selected additional realizations. These realizations are 

products of the multiplier array and the base-case hydraulic conductivity field, and do not have lake cells replaced with 

calculated bulk hydraulic conductivities.

80°30' 80°25' 80°30' 80°25' 80°30' 80°25'

25°55'

25°50'

25°45'

80°30' 80°25'

25°55'

25°50'

25°45'

25°55'

25°50'

25°45'

25°55'

25°50'

25°45'

A

B C D L-30 Canal L-30 Canal L-30 Canal

L-30 Canal

76
Figure 20. Log-normal probability distribution of effective porosity 

used for sampling of random values for Monte Carlo analysis. 

Values range from 0.04 to 0.50.

0

0.02

0.04

0.06

0.08

0.10

0.12

0.14

0.16

0.18

0.20

0.01 0.10 1.00

Probability of value

Effective porosity value

77
EXPLANATION

G-119 (tailwater)

S-337 (tailwater)

Selected day

Figure 21. Tailwater stages at S-337 and G-119 with selected days representing dry, average, 

and wet hydrologic conditions identified with black circles.

0

1

2

3

4

5

6

7

8

1996 1997 1998 1999 2000 2001 2002 2003 2004

Canal stage, in feet NAVD 88

Date

78
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

10 

30 

100 

210 

0.1-foot drawdown

0.25-foot drawdown

Study area

Lakes (2004) Production well

Particle endpoints, 

 days of travel

Simulated canals

Canals, not simulated

Figure 22. Simulated 0.25-foot and 0.1-foot drawdown contours, and particle endpoints for travel times 

of 10, 30, 100, and 210 days for the base-case simulation of dry conditions (June 19, 2004). Area shown 

is Base map 2 (fig. 1). NWWF is Northwest Well Field. WWF is West Well Field. 

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY

Snapper Creek Extension Canal

L-31 N Canal

L-30 Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

79
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

Simulated canals

Canals, not simulated

0.25-foot drawdown

Production well

10 

30 

100 

210 

0.1-foot drawdown

Study area

Lakes (2004)

Particle endpoints, 

 days of travel

Figure 23. Simulated 0.25-foot and 0.1-foot drawdown contours, and particle endpoints for travel times 

of 10, 30, 100, and 210 days for the base-case simulation of average conditions (December 19, 2002). 

Area shown is Base map 2 (fig. 1). NWWF is Northwest Well Field. WWF is West Well Field.

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY

Snapper Creek Extension Canal

L-31 N Canal

L-30 Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

80
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

10 

30 

100 

210 

0.1-foot drawdown

0.25-foot drawdown

Study area

Lakes (2004) Production well

Particle endpoints, 

 days of travel

Simulated canals

Canals, not simulated

Figure 24. Simulated 0.25-foot and 0.1-foot drawdown contours, and particle endpoints for travel times 

of 10, 30, 100, and 210 days for the base-case simulation of wet conditions (October 16, 1999). Area 

shown is Base map 2 (fig. 1). NWWF is Northwest Well Field. WWF is West Well Field.

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY

Snapper Creek Extension Canal

L-31 N Canal

L-30 Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

81
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

10 

30 

100 

210 

Study area

Lakes (2004)

Production well

Particle endpoints, 

 days of travel

Simulated canals

Canals, not simulated

Figure 25. Simulated particle endpoint locations for travel times of 10, 30, 100, and 210 days for dry conditions 

for 10,000 realizations. Area shown is Base map 2 (fig. 1). NWWF is Northwest Well Field. WWF is West Well Field.

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY

Snapper Creek Extension Canal

L-31 N Canal

L-30 Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

82
Snapper Creek Extension Canal

0°

0 1

1

2 MILES

0 2 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

Particle endpoint

Production well

Lakes (2004)

Simulated canals

Canals, not simulated

Figure 26. Simulated particle endpoint locations for 210 days of travel for dry conditions for 100 realizations 

around the Northwest Well Field (NWWF), and an example of the radial class division described by Varljen 

and Shafer (1991). Area shown is Base map 4 (fig. 1).

80°30' 80°25'

25°55'

25°50'

Tamiami Canal

Water Conservation

Area 3B

L-30 Canal

83
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 27. Median and 95-percent confidence intervals for the 210-day capture zones for the Northwest 

Well Field (NWWF) and West Well Field (WWF) for steady-state, dry conditions, and particle endpoints 

for 210 days of travel for the base-case simulation of dry conditions. Area shown is Base map 2 (fig. 1).

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY

Snapper Creek Extension Canal

L-31 N Canal

L-30 Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

84
L-31 N Canal

L-30 Canal

0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 28. Median and 95-percent confidence intervals for the 100-day capture zones for the Northwest 

Well Field (NWWF) and West Well Field (WWF) for steady-state, dry conditions, and particle endpoints 

for 100 days of travel for the base-case simulation of dry conditions. Area shown is Base map 2 (fig. 1).

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY

Snapper Creek Extension Canal

L-31 N Canal

L-30 Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

85
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 29. Median and 95-percent confidence intervals for the 30-day capture zones for the Northwest 

Well Field (NWWF) and West Well Field (WWF) for steady-state dry conditions, and particle endpoints 

for 30 days of travel for the base-case simulation of dry conditions. Area shown is Base map 2 (fig. 1). 

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

86
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 30. Median and 95-percent confidence intervals for the 10-day capture zones for the Northwest 

Well Field (NWWF) and West Well Field (WWF) for steady-state, dry conditions, and particle endpoints 

for 10 days of travel for the base-case simulation of dry conditions. Area shown is Base map 2 (fig. 1).

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

87
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 31. Median and 95-percent confidence intervals for the 210-day capture zones for the Northwest 

Well Field (NWWF) and West Well Field (WWF) for steady-state, average conditions, and particle 

endpoints for 210 days of travel for the base-case simulation of average conditions. Area shown is 

Base map 2 (fig. 1). 

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

88
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 32. Median and 95-percent confidence intervals for the 100-day capture zones for the Northwest 

Well Field (NWWF) and West Well Field (WWF) for steady-state, average conditions, and particle 

endpoints for 100 days of travel for the base-case simulation of average conditions. Area shown is 

Base map 2 (fig. 1). 

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY

EXPLANATION

L-31 N Canal

L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

89
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 33. Median and 95-percent confidence intervals for the 30-day capture zones for the Northwest 

Well Field (NWWF) and West Well Field (WWF) for steady-state, average conditions, and particle 

endpoints for 30 days of travel for the base-case simulation of average conditions. Area shown is 

Base map 2 (fig. 1).

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY

EXPLANATION

L-31 N Canal

L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

90
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 34. Median and 95-percent confidence intervals for the 10-day capture zones for the Northwest 

Well Field (NWWF) and West Well Field (WWF) for steady-state, average conditions, and particle 

endpoints for 10 days of travel for the base-case simulation of average conditions. Area shown is 

Base map 2 (fig. 1).

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

EXPLANATION

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

91
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 35. Median and 95-percent confidence intervals for the 210-day capture zones for the Northwest 

Well Field (NWWF) and West Well Field (WWF) for steady-state, wet conditions, and particle endpoints 

for 210 days of travel for the base-case simulation of wet conditions. Area shown is Base map 2 (fig. 1).

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

92
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 36. Median and 95-percent confidence intervals for the 100-day capture zones for the Northwest 

Well Field (NWWF) and West Well Field (WWF) for steady-state, wet conditions, and particle endpoints 

for 100 days of travel for the base-case simulation of wet conditions. Area shown is Base map 2 (fig. 1).

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

93
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 37. Median and 95-percent confidence intervals for the 30-day capture zones for the Northwest 

Well Field (NWWF) and West Well Field (WWF) for steady-state, wet conditions, and particle endpoints 

for 30 days of travel for the base-case simulation of wet conditions. Area shown is Base map 2 (fig. 1). 

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

94
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

95-percent CI

Median

Base-case deterministic 

 endpoint

Lakes (2004)

Production well

Simulated canals

Canals, not simulated

Figure 38. Median and 95-percent confidence intervals for the 10-day capture zones for the Northwest Well 

Field (NWWF) and West Well Field (WWF) for steady-state, wet conditions, and particle endpoints for 10 days 

of travel for the base-case simulation of wet conditions. Area shown is Base map 2 (fig. 1). 

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

95
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

Lakes (2004)

Production well

Dry 210-day upper bound 

 of 95-percent CI

Wet 210-day upper bound 

 of 95-percent CI

Simulated canals

Canals, not simulated

Figure 39. Upper bounds of the 95-percent confidence intervals (CI) for the 210-day capture zones for the 

Northwest Well Field (NWWF) and West Well Field (WWF) for steady-state, wet and dry conditions. Area 

shown is Base map 2 (fig. 1). 

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

96
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

Figure 40. Composite areal extent of the upper bounds of the 95-percent confidence intervals (CI) for 

the 210-day capture zones for the Northwest Well Field (NWWF) and West Well Field (WWF) for 

steady-state, wet and dry conditions. Also shown are the current 210-day Wellhead Protection Areas 

for the NWWF and the WWF (fig. 3). Area shown is Base map 2 (fig. 1). 

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

Lakes (2004)

Production well

210-day travel-time capture zone 

 Wellhead Protection Areas 

 (see figure 3)

Simulated canals

Canals, not simulated

NWWF composite 210-day upper 

 bound of 95-percent CI

WWF composite 210-day upper 

 bound of 95-percent CI

97
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

Lakes (2004)

Production well

0.1-foot 95-percent CI

0.1-foot median

0.25-foot 95-percent CI

0.25-foot median

Simulated canals

Canals, not simulated

Figure 41. Median and 95-percent confidence intervals (CI) for 0.1-foot and 0.25-foot drawdowns for the 

Northwest Well Field (NWWF) and West Well Field (WWF) for steady-state, dry conditions. Area shown is 

Base map 2 (fig. 1).

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

98
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data Lakes (2004)

Production well

0.1-foot 95-percent CI

0.1-foot median

0.25-foot 95-percent CI

0.25-foot median

Simulated canals

Canals, not simulated

Figure 42. Median and 95-percent confidence intervals (CI) for 0.1-foot and 0.25-foot drawdowns for the 

Northwest Well Field (NWWF) and West Well Field (WWF) for steady-state, average conditions. Area 

shown is Base map 2 (fig. 1). 

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY

EXPLANATION

Snapper Creek Extension Canal

L-31 N Canal

L-30 Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

99
0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

Lakes (2004)

Production well

0.1-foot 95-percent CI

0.1-foot median

0.25-foot 95-percent CI

0.25-foot median

Simulated canals

Canals, not simulated

Figure 43. Median and 95-percent confidence intervals (CI) for 0.1-foot and 0.25-foot drawdowns for the 

Northwest Well Field (NWWF) and West Well Field (WWF) for steady-state, wet conditions. Area shown is 

Base map 2 (fig. 1).

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY

Snapper Creek Extension Canal

L-31 N Canal

L-30 Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

100
Dade-Broward Levee

#

0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

Control structure and 

 name

Proposed recharge canal

Proposed slurry wall

HFBs

Lakes (2004)

Production well

G211#

Simulated canals

Canals, not simulated

Levee

Figure 44. Locations of the recharge canal west of the Northwest Well Field (NWWF), the proposed 

slurry wall between the L-31N Levee and the L-31N Canal, and the Horizontal Flow Barriers (HFBs) 

west of the West Well Field (WWF) for Scenario 1. Area shown is Base map 2 (fig. 1). 

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Levee L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

EVERGLADES

NATIONAL

PARK

G211

101
DDDDDDDDDDD

D

DDDDDDD

DDD

D

D

DDDDD

D

DDDDDDDD

D

DDDDDDDDDDDDD

D

D

D

DDDDD

D

D

DDDDDD

D

DDDD

D

DDDDDDDD

D

D

D

D

D

D

DDDDD DDDDDD DDDDDD DDDDDDD DDD DDDDDDD

DD

DDDD

D

D

D

D

D

DD

D

DDD DDD DDD DD D DDDDD DD DD DDDDDDDDDDDD DDDDDD DDD DD D DD D

D

D

D

D

D

D

D

D

D

D

D

D

DDD

DD

D

D

D

D

DD

D

D

D

D

DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD

D

DD

DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD

D

DDDDDD

DDDDDDDD

DDD

D

DDDDDD

D

DDD

DDDDD

D

D

D

D

D

D

DDD

D

D

D

DD

D

D

D

D

D

D

DDDDDD

D

DD

D DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD D

D

D

D

D

D

D

D

D

D

D

DD

D

D

D

DD

D

D

D

DD

DD

DDD

D

DDDDDDDDDDDDDDDDDDDDD

DD

D

DDDDDD

DDD

DDDDDDDDD

0 1 2 MILES

0 1 2 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

Figure 45. Median and 95-percent confidence intervals for the 210-day capture zones for the 

Northwest Well Field (NWWF) for steady-state, dry conditions for seepage mitigation scenario 

and without seepage mitigation, and final particle endpoints for both conditions beyond 210 days. 

Area shown is Base map 4 (fig. 1).

80°30' 80°25'

25°55'

25°50'

NWWF

Simulated canals

Canals, not simulated

Lakes (2004)

Production well

95-percent CI for 210-day capture zone 

 for seepage mitigation scenario

Median 210-day capture zone for 

 seepage mitigation scenario

95-percent CI for 210-day capture zone 

 without seepage mitigation 

Median 210-day capture zone without 

 seepage mitigation 

Seepage mitigation particle endpoint

Proposed recharge canal

D Dry condition particle endpoint

Tamiami Canal

Water Conservation

Area 3B

Snapper Creek Extension Canal

L-30 Canal

102
Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

0 1

1

2 MILES

0 2 KILOMETERS

Upper bound of 95-percent Cl for 210-day capture 

 zone for seepage mitigation scenario

Upper bound of 95-percent Cl for 210-day capture 

 zone without seepage mitigation 

Upper bound of the 95-percent Cl for 0.1-foot 

 drawdown for seepage mitigation scenario

Upper bound of the 95-percent Cl for 0.1-foot 

 drawdown without seepage mitigation 

Lakes (2004)

Production well

G211 Control structure and name

Simulated canals

Canals, not simulated

HFBs

G211

Figure 46. Upper bound of the 95-percent confidence intervals (CI) for the 210-day capture zones and 

upper bound of the 95-percent CI for the 0.1-foot drawdown for the West Well Field (WWF) for 

steady-state, dry conditions for seepage mitigation scenario and without seepage mitigation. Area 

shown is Base map 5 (fig. 1).

80°30' 80°25'

25°45'

25°40'

L-31 N Canal

Tamiami Canal

WWF

C-1 W Canal

EVERGLADES

NATIONAL

PARK

103
#

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

0 1

1

2 MILES

0 2 KILOMETERS

G211

EXPLANATION

Lakes (2004)

Production well

G211 Control structure and name

Upper bound of 95-percent Cl for 100-day capture 

 zone for seepage mitigation scenario

Upper bound of 95-percent Cl for 100-day capture 

 zone without seepage mitigation 

Simulated canals

Canals, not simulated

HFBs

Figure 47. Upper bound of the 95-percent confidence intervals (CI) for the 100-day capture zones for the 

West Well Field (WWF) for steady-state, dry conditions for seepage mitigation scenario and without 

seepage mitigation. Area shown is Base map 5 (fig. 1). 

80°30' 80°25'

25°45'

25°40'

L-31 N Canal

Tamiami Canal

WWF

C-1 W Canal

EVERGLADES

NATIONAL

PARK

104
5

56

79

70

6

64

67

63

7

93

24

16

26

81

31

44

62

84

2

83

80

28

71

69

1

68

23

30 48

82

92

61

94

85

16

13

87

34

75

19

51

32

13

52

29

12 76

40

41

3

91

15

59

66

60

14

65

33

22

95 53

17

43

46 88 45 52

86

90

54

42

27

18

47

20

72

11

10

59

0 2 4 MILES

0 2 4 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

Lakes (proposed expansion)

Production well

Figure 48. Locations of existing lakes, projected expansion of existing lakes and new proposed lakes for 

lake-expansion scenario. Area shown is Base map 2 (fig. 1). NWWF is Northwest Well Field. WWF is 

West Well Field. 

80°35' 80°30' 80°25' 80°20'

25°55'

25°50'

25°45'

25°40'

NWWF

WWF

Water 

Conservation

Area 3B

BROWARD COUNTY

MIAMI-DADE COUNTY L-31 N Canal L-30 Canal

Snapper Creek Extension Canal

Tamiami Canal

L-67A Canal

L-67C Canal

L-29 Canal

Simulated canals

Canals, not simulated

EVERGLADES

NATIONAL

PARK

105
0 1 2 MILES

0 1 2 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

EXPLANATION

Production well

Median 210-day capture zone for lake 

 expansion scenario

Upper bound of 95-percent for 210-day capture 

 zone for lake expansion scenario Upper bound of the 95-percent Cl for 0.25-foot 

 drawdown for lake expansion scenario

Upper Bound of 95-percent Cl for 210-day capture 

 zone without lake expansion

Median 210-day capture zone without 

 lake expansion

Upper bound of the 95-percent Cl for 0.25-foot 

 drawdown without lake expansion

Lakes (proposed expansion)

Simulated canals

Canals, not simulated

Figure 49. Median and upper bound of the 95-percent confidence intervals (CI) for the 210-day capture 

zones and upper bound of the 95-percent CI for the 0.25-foot drawdown for the Northwest Well Field 

(NWWF) for steady-state dry conditions for the lake expansion scenario and without lake expansion. 

Area shown is Base map 4 (fig. 1). 

80°30' 80°25'

25°55'

25°50'

NWWF

Tamiami Canal

Water Conservation

Area 3B

Snapper Creek Extension Canal

L-30 Canal

106
0 1

1

2 MILES

0 2 KILOMETERS

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

Upper bound of 95-percent Cl for 210-day capture 

 zone for lake expansion scenario

Upper bound of 95-percent Cl for 210-day capture

 zone without lake expansion

Production well

Upper bound of 95-percent Cl for 0.1-foot 

 drawdown for lake expansion scenario

Upper bound of the 95-percent Cl for 0.1-foot 

 drawdown without lake expansion

Lakes (proposed expansion)

EXPLANATION

Simulated canals

Canals, not simulated

Figure 50. Upper bound of the 95-percent confidence intervals (CI) for the 210-day capture zones and 

upper bound of the 95-percent CI for the 0.1-foot drawdown for the West Well Field (WWF) for 

steady-state, dry conditions for the lake expansion scenario and without lake expansion. Area shown 

is Base map 5 (fig. 1).

80°30' 80°25'

25°45'

25°40'

WWF

L-31 N Canal

Tamiami Canal

C-1 W Canal

EVERGLADES

NATIONAL

PARK

107
500-day particle pathline

Head contour

Constant-head cells

Simulated head, in feet

0

-4.5

-9.0

EXPLANATION 0 5 10 MILES

0 5 10 KILOMETERS

-0.5

-0.5

-0.5

-0.5

-1.5

-1.5

-2.0

-2.0

-2.5

-2.5

-1.0

-1.0

-1.0

-1.0

Figure 1–1. Simulated 500-day particle pathlines and head contours for one-layer conceptual Model A (left) and three-layer 

conceptual Model B (layer 2) (right).

108
0 2.5 5 10 Miles

0 5 10 MILES

0 5 10 KILOMETERS

5

5

5

5

5

5

EXPLANATION

5 Pumping well

River cells

Constant-head cells

Model grid

Figure 1–2. Model grid, locations of constant-head cells, river cells and pumping wells for 

Model A and Model B.

0 2 4 MILES

0 2 4 KILOMETERS

109
EXPLANATION

0 5 10 MILES

0 5 10 KILOMETERS

5

5

5

5

5

5

River cells

Constant-head cells

95-percent Cl for 1-foot 

 drawdown for Model A

Median 1-foot drawdown for 

 Model A

5 Pumping well

95-percent Cl for 1-foot 

 drawdown for Model B

Median 1-foot drawdown for 

 Model B

Figure 1–3. Median and 95-percent confidence intervals (CI) for 

1-foot drawdown for Model A and Model B.

110
5

5

5

5

5

5

Figure 1–4. Median and upper bound of 95-percent confidence interval (CI) for 210-day capture zones for Model A and Model B.

0 2 4 MILES

0 2 4 KILOMETERS

EXPLANATION

Upper bound of 95-percent Cl of 210-day 

 capture zone for Model A

Median 210-day capture zone for 

 Model A

Upper bound of 95-percent Cl of 210-day 

 capture zone for Model B

Median 210-day capture zone for 

 Model B

Model grid

5 Pumping well

River cells

111
A

B

C

Figure 3–1. Simulated water levels for base-case, steady-state average conditions 

(black area represents inactive model cells) (A) with model boundary as described 

in the main body of the report, (B) with model boundary extended west of the western 

canals and general head boundary conductance of 269 square feet per day (ft2

/d) 

(25 square meters per day (m2

/d)), and (C) with model boundary extended west of the 

western canals and general head boundary conductance of 26,910 ft2

/d (2,500 m2/d).

EXPLANATION

Inactive model cells

Area of standing surface 

 water in Everglades

Simulated canals

Water-level contour, 

 contour interval 0.5 foot

112
0 5 MILES

0 5 KILOMETERS

2.5

2.5

2.2

2.4

2.6 2.8

0.8

0.8

-0.8

-0.2

-0.6

-0.4

-1.0

0

2.0

1.2

1.2

1.8

1.6

0.6

0.6

1.6

1.8

0.4

0.2

2.0

1.0

1.0

1.4

1.4

1.4

2.2

2.4

2.6 2.8

0.8

0.8

-0.8

-0.2

-0.6

-0.4

-1.0

2.0

1.2

1.2

1.8

1.6

0.6

0.6

1.6

1.8

0.4

0.2

2.0

1.0

1.0

1.4

1.4

EXPLANATION

80°28' 80°26' 80°24' 80°22' 80°20'

25°54'

25°52'

25°50'

25°48'

25°54'

25°52'

25°50'

25°48'

80°28' 80°26' 80°24' 80°22' 80°20'

A

B

Figure 4–1. Simulated steady-state groundwater head and backward-tracked 

particle pathlines for (A) the model representing lakes with high hydraulic 

conductivity and (B) the model using the LAK Package to represent lakes. Lake 

numbers correspond to data in Table 4-1.

L-30

C-6

C-7 EXT

C-6

L-30

C-6

C-7 EXT

C-6

0 2 4 MILES

0 2 4 KILOMETERS

5

Production well

Constant-head cells

Lake and number

Head contour

Particle pathline

5

6

24

8

12

3

9

32

37

7

10

1

18

11

28

25

30

22

36

29

20

27

33

24

21

38

34

17

19

31

23

26

15

35

13 14

16

0

6 5

24

8

12

3 4

4

9

32

37

7

10

1

18

11

28

25

30

22

36

29

20

27

33

24

21

38

34

17

19

31

23

26

15

35

1314

16

Base modified from U.S. Geological Survey 

1:2,000,000-scale digital data

113
G-1488 G-1637

G-3264A G-3439

G-3551

Head (NAVD 88), in feet Head (NAVD 88), in feet Head (NAVD 88), in feet

Head (NAVD 88), in feet Head (NAVD 88), in feet Head (NAVD 88), in feet

G-3552

Figure 5–1. Daily observed and model-simulated water levels for monitoring wells in the study area not shown 

in figure 13.

EXPLANATION

Simulated groundwater heads

Observed groundwater heads

Date Date

Date Date

Date

1996 1997 1998 1999 2000 2001 2002 2003 2004 1996 1997 1998 1999 2000 2001 2002 2003 2004

1996 1997 1998 1999 2000 2001 2002 2003 2004 1996 1997 1998 1999 2000 2001 2002 2003 2004

1996 1997 1998 1999 2000 2001 2002 2003 2004 1996 1997 1998 1999 2000 2001 2002 2003 2004

Date

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

5.5

6.0

6.5

1.5

1.5

1.0

1.0

0.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

5.5

3.0

3.0

2.5

2.5

2.0

2.0

3.5

3.5

4.0

4.0

4.5

4.5

5.0

5.0

5.5

6.0

6.5

7.0

1

2

3

4

5

6

7

1

2

3

4

5

6

7

114
G-3553 G-3554

G-3555 G-3556

G-3558 G-3559

Date Date

Date Date

Date Date

3

2

1

4

5

6

7

3

2

1

4

5

6

7

3

2

1996 1997 1998 1999 2000 2001 2002 2003 2004 1996 1997 1998 1999 2000 2001 2002 2003 2004

1996 1997 1998 1999 2000 2001 2002 2003 2004 1996 1997 1998 1999 2000 2001 2002 2003 2004

1996 1997 1998 1999 2000 2001 2002 2003 2004 1996 1997 1998 1999 2000 2001 2002 2003 2004

1

4

5

6

7

3

2

1

4

5

6

7

1.5

2.5

3.5

4.5

5.5

6.5

2.0

3.0

4.0

5.0

6.0

1

2

3

4

5

6

7

8

Figure 5–1. Daily observed and model-simulated water levels for monitoring wells in the study area not shown 

in figure 13.—Continued.

Head (NAVD 88), in feet Head (NAVD 88), in feet Head (NAVD 88), in feetHead (NAVD 88), in feet Head (NAVD 88), in feet Head (NAVD 88), in feet

EXPLANATION

Simulated groundwater heads

Observed groundwater heads

115
G-3560 G-3561

G-3572 G-3574

G-970 G

KROME

G-973

Date Date

Date Date

Date

Date

1996 1997 1998 1999 2000 2001 2002 2003 2004 1996 1997 1998 1999 2000 2001 2002 2003 2004

1996 1997 1998 1999 2000 2001 2002 2003 2004 1996 1997 1998 1999 2000 2001 2002 2003 2004

1996 1997 1998 1999 2000 2001 2002 2003 2004

1996 1997 1998 1999 2000 2001 2002 2003 2004

1996 1997 1998 1999 2000 2001 2002 2003 2004

Date

3

2

1

4

5

6

7

3

2

1

4

5

6

7

2.5

3.5

3.5

4.5

4.5

5.5

6.5

7.0

2.0

2.5

2.0

1.5

1.0

0.5

0

3.0

3.0

4.0

4.0 3.5

4.5

2.5

2.0

1.5

1.0

0.5

0

3.0

4.0

5.0

5.0

6.0

1

2

3

4

5

6

7

8

9

1.5

2.5

3.5

4.5

5.5

6.5

2.0

3.0

4.0

5.0

6.0

Head (NAVD 88), in feet Head (NAVD 88), in feet Head (NAVD 88), in feet Head (NAVD 88), in feetHead (NAVD 88), in feet Head (NAVD 88), in feet Head (NAVD 88), in feet

EXPLANATION

Simulated groundwater heads

Observed groundwater heads

Figure 5–1. Daily observed and model-simulated water levels for monitoring wells in the study area not shown 

in figure 13.—Continued.

116
EXPLANATION

Lake 1

80

70

60 57

50

40

30

24 24

24

24

24

34 37

47

52

20 20

10 10

Residence time = 10

Residence time = 5

Lake 1

Lake 2 Lake 2

Figure 6–1. Backward-tracked particle pathlines from a pumping well 

through lakes (A) without residence-time post-processing, (B) with 

residence-time post-processing.

Production well

Lake

MODPATH backward-tracked 

 particle pathline

MODPATH particle travel time

Adjusted MODPATH particle 

 travel time

20

24

A B

117
Figure 6–2. Map showing locations and numbering of lakes referred to in paper by Chin and 

others (2010). Lake numbers correspond to data in Table 6–1. Area shown is Base map 6 (fig. 1).

22 23

28

33

30 32

80°27' 80°24'

25°48'

25°51'

NWWF

EXPLANATION

Lake and number

Production well

Simulated canals

33

0 1 MILE

0 0.5 1 KILOMETER

Base modified from U.S. Geological Survey 0.5

1:2,000,000-scale digital data

Snapper Creek Extension Canal

118

Table 1. Monitoring well names, mean errors of head, mean absolute errors of head, and the coefficients of determination of

simulated versus observed values of head.

[ME is plotted in figure 14. ME, mean error; MAE, mean absolute error; R2, coefficient of determination]

Well name

ME

(feet)

MAE

(feet) R2 Well name

ME

(feet)

MAE

(feet) R2

G-1488 0.222 0.265 0.915 G-3559 0.115 0.126 0.967

G-1637 0.439 0.464 0.800 G-3560 -0.392 0.429 0.821

G-3253 0.526 1.417 0.443 G-3561 -0.293 0.333 0.823

G-3259A 0.205 0.911 0.324 G-3572 -0.653 0.665 0.781

G-3264A -0.038 0.372 0.689 G-3574 0.605 0.612 0.790

G-3439 -0.486 0.503 0.883 G-551 2.139 2.186 0.519

G-3551 0.026 0.220 0.908 G-855 -0.316 0.364 0.832

G-3552 -0.338 0.430 0.739 G-970_G -0.029 0.149 0.834

G-3553 -0.017 0.384 0.668 G-973 -0.065 0.377 0.335

G-3554 -0.007 0.269 0.839 G-975 -0.139 0.330 0.761

G-3555 -0.569 0.581 0.831 G-976 0.141 0.555 0.704

G-3556 -0.256 0.343 0.800 KROME 0.012 0.087 0.954

G-3558 -0.614 0.619 0.859

120
Table 2. Statistics of minimum value and maximum value of hydraulic

conductivity for all realizations, prior to calculation of bulk hydraulic

conductivities.

[Values shown are in feet per day]

Descriptor Minimum value for all

realizations

Maximum value for all

realizations

2.5-percentile value 60.5 104,205

Median value 109 147,672

97.5-percentile value 159 256,076

121
Lake

High hydraulic conductivity lakes LAK Package lakes

Groundwater 

inflow

(ft

3

/d)

Evaporation

(ft

3

/d)

Groundwater 

outflow 

(ft

3

/d)

Groundwater 

inflow

(ft3/d)

Evaporation

(ft

3

/d)

Groundwater 

outflow 

(ft

3

/d)

1 504,399 33,514 470,886 510,968 33,514 477,631

2 1,222,029 79,423 1,142,606 1,267,902 79,423 1,188,798

3 964,196 61,448 902,714 1,062,512 61,448 1,001,347

4 1,183,465 68,828 1,114,672 1,237,108 68,828 1,168,704

5 3,389,502 171,982 3,217,661 3,525,322 171,982 3,353,940

Table 4–1. Simulated groundwater inflow, groundwater outflow, and evaporation lake water-budget terms for lakes simulated using 

the high-hydraulic-conductivity approach and the LAK Package. Lakes are shown in figure 4-1.

[ft

3

/d, cubic foot per day]

6 3,104,089 126,709 2,977,591 3,116,696 126,709 2,990,552

7 1,712,655 43,084 1,669,571 1,733,526 43,084 1,690,619

8 3,921,941 64,661 3,857,951 3,993,736 64,661 3,929,463

9 2,513,239 44,885 2,468,213 2,561,373 44,885 2,516,770

10 4,512,508 84,791 4,427,047 4,582,078 84,791 4,497,676

11 5,996,819 57,669 5,939,009 6,055,406 57,669 5,998,196

12 4,255,205 123,495 4,131,675 4,395,970 123,495 4,272,722

13 917,016 3,178 913,838 1,001,453 3,178 998,381

14 684,186 3,002 681,185 695,664 3,002 692,768

15 756,122 5,580 750,543 730,484 5,580 724,975

16 548,896 2,578 546,283 563,693 2,578 561,185

17 1,265,784 12,360 1,253,423 1,283,829 12,360 1,271,575

18 3,970,675 65,262 3,905,414 3,970,075 65,262 3,905,096

19 1,292,340 10,771 1,281,569 1,341,251 10,771 1,330,621

20 2,244,706 31,112 2,213,594 2,303,682 31,112 2,272,746

21 2,534,710 21,542 2,513,133 2,552,615 21,542 2,531,249

22 3,238,355 34,326 3,204,064 3,240,650 34,326 3,206,642

23 1,793,067 9,994 1,783,108 1,806,345 9,994 1,796,492

24 1,940,082 28,322 1,911,759 1,944,920 28,322 1,916,809

25 2,003,401 39,093 1,964,308 2,063,542 39,093 2,024,660

26 1,174,707 8,370 1,166,337 1,188,550 8,370 1,180,287

27 1,328,114 26,733 1,301,381 1,349,515 26,733 1,322,993

28 776,464 40,506 735,922 773,497 40,506 733,238

29 1,093,625 32,313 1,061,312 1,097,121 32,313 1,065,020

30 1,670,808 38,493 1,632,315 1,700,754 38,493 1,662,579

31 425,365 10,559 414,771 441,469 10,559 431,051

32 2,065,025 89,381 1,975,679 2,102,918 89,381 2,013,854

33 424,094 23,131 400,927 419,397 23,131 396,407

34 1,531,986 14,550 1,517,436 1,546,465 14,550 1,532,021

35 610,308 5,403 604,940 613,239 5,403 607,907

36 928,811 33,514 895,262 967,975 33,514 934,567

122
37 1,687,264 89,381 1,597,883 1,745,604 89,381 1,656,646

38 678,607 16,351 662,221 696,864 16,351 680,549

 Sum 70,864,529 1,656,470 69,208,059 72,184,168 1,656,470 70,536,738

Table 4–1. Simulated groundwater inflow, groundwater outflow, and evaporation lake water-budget terms for lakes simulated using 

the high-hydraulic-conductivity approach and the LAK Package. Lakes are shown in figure 4-1.—Continued

[ft

3

/d, cubic foot per day]

Lake

High hydraulic conductivity lakes LAK Package lakes

Groundwater 

inflow

(ft

3

/d)

Evaporation

(ft

3

/d)

Groundwater 

outflow 

(ft

3

/d)

Groundwater 

inflow

(ft3/d)

Evaporation

(ft

3

/d)

Groundwater 

outflow 

(ft

3

/d)

123
Lake 

number

Residence time (days)

This study Chin and others (2010)

Table 6-1. Residence times used in model documented in this report and those calculated

by Chin and others (2010). Lakes are shown in figure 6-2.

22 239 314

23 325 1,290

28 115 259

30 165 88

32 73 87

33 99 159

124


Brakefield and others—Estimation of Capture Zones and Drawdown at the Northwest and West Well Fields—Open-File Report 2013–1086