Prepared in cooperation with the U.S. Navy

Groundwater Withdrawals and Regional Flow Paths at

and near Willow Grove and Warminster, Pennsylvania

—Data Compilation and Preliminary Simulations for

Conditions in 1999, 2010, 2013, 2016, and 2017

By Daniel J. Goode and Lisa A. Senior 

Open-File Report 2019–1137 

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

DAVID BERNHARDT, Secretary

U.S. Geological Survey

James F. Reilly II, Director

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

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

https://www.usgs.gov/ or call 1–888–ASK–USGS (1–888–275–8747). 

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

visit https://store.usgs.gov/.

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.

Suggested citation:

Goode, D.J., and Senior, L.A., 2020, Groundwater withdrawals and regional flow paths at and near

Willow Grove and Warminster, Pennsylvania—Data compilation and preliminary simulations for

conditions in 1999, 2010, 2013, 2016, and 2017: U.S. Geological Survey Open-File Report 2019–1137, 

127 p., https://doi.org/10.3133/ofr20191137.

Associated data for this publication:

Goode, D.J., and Senior, L.A., 2020, MODFLOW 6 and MODPATH 7 model data sets used to evaluate 

groundwater flow in the vicinity of Horsham and Warminster, Bucks and Montgomery Counties,

Pennsylvania—Preliminary simulations for conditions in 1999, 2010, 2013, 2016, and 2017: U.S.

Geological Survey data release, https://doi.org/10.5066/P9K36P5S. 

Senior, L.A., Ludlow, R.A., and Zarr, L.F, 2020, Groundwater levels, groundwater withdrawals, and

point-source discharges to streams in the vicinity of Willow Grove and Warminster, Bucks and

Montgomery Counties, Pennsylvania, for selected years during 1999–2017: U.S. Geological Survey data 

release, https://doi.org/10.5066/P9ZGEI67. 

ISSN 2331-1258 (online) 
iii

Acknowledgments 

The support of U.S. Navy personnel and their contractors, Resolution Consultants, Tetra Tech

Incorporated, and Battelle Memorial Institute is gratefully acknowledged. Data and other technical support, 

as well as permission to access wells from the Horsham Water and Sewer Authority and the Warminster 

Township Municipal Authority is gratefully acknowledged. Additionally, water-level data were provided by 

the Air National Guard and their contractor Leidos. William McLay, Montgomery County, Pennsylvania, and 

Robert H. Keough, Bucks County, Pennsylvania, provided stream geographic information system (GIS) 

data. 

Numerous U.S. Geological Survey colleagues provided assistance to the project as follows: 

Michael Fienen and Christian Langevin for Python instructions and scripts for Pest++ calibration; Richard 

Winston for ModelMuse development and support; Linda Zarr and Russ Ludlow for database compilation; 

Matthew Conlon for data processing; Matthew Gyves and Scott Hoffman for GIS support; Dennis Risser for 

groundwater modeling and data analysis guidance; and Richard Yager and Paul Misut for colleague 

technical reviews.
iv

Contents

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

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

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

Previous Investigations .............................................................................................................................. 5

Purpose and Scope.................................................................................................................................... 5

Description of Study Area.............................................................................................................................. 6

Geologic Setting......................................................................................................................................... 6

Hydrologic Setting .................................................................................................................................... 10

Hydrologic Conditions and Water Use, 1999–2017 ..................................................................................... 11

Precipitation ............................................................................................................................................. 12

Streamflow ............................................................................................................................................... 14

Groundwater Levels ................................................................................................................................. 22

Groundwater Withdrawals ........................................................................................................................ 27

Discharges to Streams............................................................................................................................. 32

Simulation of Regional Groundwater Flow................................................................................................... 33

Conceptual Model of the Groundwater-Flow System............................................................................... 33

Model Approach and Methods ................................................................................................................. 34

Model Area and Spatial Discretization...................................................................................................... 34

Model Hydrogeologic Zones..................................................................................................................... 39

Aquifer Inflows and Outflows.................................................................................................................... 40

Recharge ............................................................................................................................................. 40

Groundwater Withdrawals.................................................................................................................... 41

Flow to and from Streams.................................................................................................................... 44

Bulk Hydraulic Conductivity and Anisotropy ............................................................................................. 46

Model Calibration ..................................................................................................................................... 47

Effects of Recharge and Withdrawals on Groundwater-Flow Paths ......................................................... 52

1999 Conditions................................................................................................................................... 53

2017 Conditions................................................................................................................................... 58

Groundwater-Flow Paths and PFAS Contamination in Wells................................................................... 60

Stream Loss to Groundwater ................................................................................................................... 66

Limitations and Suggestions for Improvements ........................................................................................... 68

Summary ..................................................................................................................................................... 69

References Cited......................................................................................................................................... 71

Appendix 1. Annual Base Flow as Determined from Measured Streamflow at Selected Gages and 

Estimated for Missing Streamflow Records During 2010–2015............................................................ 79

Appendix 2. Model Calibration Results ........................................................................................................ 91

Appendix 3. Simulated Water Levels and Groundwater-Flow Paths.......................................................... 101
v 

Figures

1. Locations of former Willow Grove Naval Air Station Joint Reserve Base, former Naval Air

Warfare Center Warminster, and active Horsham Air Guard Station, land-surface altitudes,

streams, and model grid boundary, Bucks and Montgomery Counties, Pennsylvania. .............. 3

2. Physiographic provinces of Pennsylvania and physiographic-province sections in the study

area in southeastern Pennsylvania. ........................................................................................... 7

3. Geology of the study area and major faults, streams, selected U.S. Geological Survey stream

gages and observation wells, and location of model grid, Bucks and Montgomery Counties,

Pennsylvania.............................................................................................................................. 8

4. Long-term daily mean depth to water in U.S. Geological Survey observation well BK-1020,

completed in the Stockton Formation, Warminster Township, Bucks County, Pennsylvania. .. 11

5. Mean annual total precipitation for Bucks and Montgomery Counties, Pennsylvania,

1990–2017. .............................................................................................................................. 12

6. Location of selected U.S. Geological Survey stream gages and long-term observation wells,

wells with water levels used in model calibration, selected point-source discharges to streams

above and near gages where data were used in model calibration, and model grid. ............... 15

7. Data collected at Little Neshaminy Creek at Valley Road and Pennypack Creek gages, where

base flow was determined using the HYSEP hydrograph separation method from the U.S.

Geological Survey Groundwater Toolbox: A, Daily base flow and total streamflow,

1999–2018; B, Monthly base-flow and runoff components of streamflow, 2014–17; C, monthly

base flow and sum of reported monthly discharges from sewage treatment plants above the

gage, 2014–17; and D, base flow at Pennypack Creek at Horsham, Pennypack Creek tributary

at Hatboro, and Pennypack Creek at Willow Grove with reported discharges from National

Pollution Discharge Elimination System sources above Hatboro gage and Upper Moreland

Hatboro STP above Willow Grove gage, 2010–16................................................................... 18

8. Groundwater levels in well BK-1020 and five Montgomery County network wells, as well as

groundwater-level data used for model calibration in the time periods selected to represent

calendar years A, 1999, and B, 2010, 2013, and 2016 conditions. .......................................... 23

9. A, Distribution of well depths and B, relation between well depths and reported pumping rates

for withdrawal wells in 1999 and 2010 in the model area near Willow Grove and Warminster,

Pennsylvania............................................................................................................................ 28

10. Total annual average reported well and quarry withdrawals in the model area near Willow

Grove and Warminster, Pennsylvania, in 1995, 1999, 2010, 2013, 2014, 2016, and 2017...... 30

11. ModelMuse screenshot of full model grid with model streams near Willow Grove and

Warminster, Pennsylvania. ...................................................................................................... 35

12. ModelMuse screenshot of detail of model grid cells and withdrawal well locations close to the

bases near Willow Grove and Warminster, Pennsylvania........................................................ 36

13. ModelMuse screenshot of details near Willow Grove of model grid in plan and cross section

showing seven model layers with illustrative regional groundwater-flow paths from possible

sources to discharge cells........................................................................................................ 37

14. ModelMuse screenshot of model hydrogeologic zones: Upper Stockton, Middle and Lower

Stockton, Lockatong, Brunswick, Carbonate, and Other.......................................................... 39
vi 

15. ModelMuse screenshot of model recharge rates in 2010 and 2013 in inches per year, zoned by

geologic unit. Recharge rates for other periods were reduced to 90 percent of the values

shown for 1999 and 70 percent for 2016 and 2017.................................................................. 40

16. ModelMuse screenshot of model groundwater withdrawal rates visualized as well-withdrawal

footprints for simulation of conditions in 1999. ......................................................................... 41

17. ModelMuse screenshot of model groundwater withdrawal rates visualized as well-withdrawal

footprints for simulation of conditions in 2017. ......................................................................... 42

18. ModelMuse screenshot of changes in groundwater withdrawal rates in the Willow Grove

and Warminster area visualized by the difference between well-withdrawal footprints for

1999 and 2017. ........................................................................................................................ 43

19. ModelMuse screenshot of locations of model point discharges to streams.............................. 45

20. ModelMuse screenshot showing model down dip angles to the northwest for hydraulic

conductivity tensor in model layers 3 to 7, which represent unweathered fractured rock......... 46

21. Simulated and measured water levels in wells for optimum parameters identified by calibration

for each model period, 1999, 2010, 2013, and 2016................................................................ 49

22. Histogram of water-level residuals for all calibration data for all model periods. ...................... 50

23. Simulated and measured streamflow for optimum parameters identified by calibration for all

model periods. ......................................................................................................................... 51

24. ModelMuse screenshot of simulated potentiometric surface in model layer 1 for

1999 conditions. ....................................................................................................................... 53

25. ModelMuse screenshot of simulated areas that contribute recharge to discharging wells and

streams for 1999 conditions. .................................................................................................... 54

26. ModelMuse screenshot of simulated regional groundwater-flow paths from locations

of possible per- and polyfluoroalkyl substances contamination in groundwater for

1999 conditions. ....................................................................................................................... 55

27. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 1999

conditions near Willow Grove. .................................................................................................. 56

28. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 1999

conditions near Warminster. .................................................................................................... 57

29. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2017

conditions near Willow Grove................................................................................................... 58

30. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2017

conditions near Warminster. .................................................................................................... 59

31. Private-well sampling areas and levels of perfluorooctane sulfonate and perfluorooctanoic

acid in well-water samples relative to drinking-water health advisories for Horsham Air

Guard Station and the former Naval Air Station Joint Reserve Base Willow Grove as of

September, 2019...................................................................................................................... 61

32. Private-well sampling areas and levels of perfluorooctane sulfonate and perfluorooctanoic acid

in well-water samples relative to drinking-water health advisories for the former Naval Air

Warfare Center Warminster as of January, 2019............ ....................................................... 62

33. Composite simulated regional groundwater-flow paths from possible per- and

polyfluoroalkyl substances contamination sources at Willow Grove and Warminster for

conditions in 1999, 2010, 2013, 2016, and 2017 and maximum levels of perfluorooctane
vii

sulfonate and perfluorooctanoic acid in 2014–2018 well-water samples relative to drinkingwater health advisories. ........................................................................................................... 63

34. Detail at Willow Grove of composite simulated regional groundwater-flow paths from possible

per- and polyfluoroalkyl substances contamination sources for conditions in 1999, 2010, 2013,

2016, and 2017 and maximum levels of perfluorooctane sulfonate and perfluorooctanoic acid

in 2014–2018 well-water samples relative to drinking-water health advisories. ....................... 64

35. Detail at Warminster of composite simulated regional groundwater-flow paths from possible

per- and polyfluoroalkyl substances contamination sources for conditions in 1999, 2010, 2013,

2016, and 2017 and maximum levels of perfluorooctane sulfonate and perfluorooctanoic acid

in 2014–2018 well-water samples relative to drinking-water health advisories. ....................... 65

36. Detail near Willow Grove, Park Creek, and Little Neshaminy Creek of ModelMuse screenshot

of simulated potentiometric surface in model layer 1, regional flow paths from possible per- and

polyfluoroalkyl substances contamination in groundwater, and locations of stream loss for

conditions in 1999. ................................................................................................................... 66

37. ModelMuse screenshot detail near Willow Grove, Park Creek, and Little Neshaminy Creek of

the simulated potentiometric surface in model layer 1, regional flow paths from possible perand polyfluoroalkyl substances contamination in groundwater, and locations of stream loss for

conditions in 2017. ................................................................................................................... 67
viii

Tables

1. Annual total county-wide precipitation for Bucks and Montgomery Counties, Pennsylvania,

1990–2017. . ............................................................................................................................... 13

2. Period of record and type of data recorded for selected U.S. Geological Survey stream gages

and observation wells in the study area in and near Willow Grove and Warminster,

Pennsylvania, 1990–2017......... ................................................................................................ 16

3. Annual average base flow at selected stream gages near Willow Grove and Warminster,

Pennsylvania, 1999 and 2010–2017. .......................................................................................... 21

4. Summary statistics for wells with water levels used in model calibration scenarios for

1999, 2010, 2013, and 2016 conditions, including land-surface altitude, depth to top and bottom

of well opening, date of water-level measurement, and depth to water near Willow Grove and

Warminster, Pennsylvania........................................................................................................... 25

5. Summary statistics for pumping wells with withdrawals used in model scenarios for

1999, 2010, 2013, 2016 and 2017 conditions, including land-surface altitude and depth to top

and bottom of well opening, near Willow Grove and Warminster, Pennsylvania. ........................ 29

6. Annual average withdrawals for wells and quarries in model area for selected years, 1995–2017,

near Willow Grove and Warminster, Pennsylvania...................................................................... 31

7. Annual average discharge rates for point source discharges to streams in model area for

model scenarios for 1999, 2010, 2013, 2016, and 2017 conditions near Willow Grove and

Warminster, Pennsylvania........................................................................................................... 32

8. Model layer thicknesses and hydraulic conductivity variation with depth..................................... 38

9. Streamflow measurements used for model calibration................................................................ 48

10. Simulated groundwater budget.................................................................................................... 52
ix

Conversion Factors

U.S. customary units to International System of Units

Multiply By To obtain

Length

inch (in.) 25.4 millimeter (mm)

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Area

square foot (ft2) 0.09290 square meter (m2)

square mile (mi2) 2.590 square kilometer (km2) 

Volume

gallon (gal) 3.785 liter (L) 

gallon (gal) 0.003785 cubic meter (m3) 

million gallons (Mgal) 3,785 cubic meter (m3)

cubic foot (ft3) 0.02832 cubic meter (m3) 

Flow rate

cubic foot per second (ft3/s) 0.02832 cubic meter per second (m3/s)

gallon per minute (gal/min) 0.06309 liter per second (L/s)

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)

Hydraulic conductivity

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

International System of Units to U.S. customary units 

Multiply By To obtain

Length

millimeter (mm) 0.03937 inch (in.)

meter (m) 3.281 foot (ft) 

kilometer (km) 0.6214 mile (mi)

Area

square meter (m2) 10.76 square foot (ft2) 

square kilometer (km2) 0.3861 square mile (mi2)

Volume

liter (L) 0.2642 gallon (gal)

cubic meter (m3) 264.2 gallon (gal) 

liter (L) 61.02 cubic inch (in3) 

cubic meter (m3) 35.31 cubic foot (ft3)

Flow rate

meter per second (m/s) 3.281 foot per second (ft/s) 

cubic meter per second (m3/s) 35.31 cubic foot per second (ft3/s)

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

cubic meter per day (m3/d) 264.2 gallon per day (gal/d) 

Hydraulic conductivity

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

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as °F = (1.8 × °C) + 32.

Temperature in degrees Fahrenheit (°F) may be converted to degrees Celsius (°C) as °C = (°F – 32) / 1.8.
x

Datum

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

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

Altitude, as used in this report, refers to distance above the vertical datum.

Abbreviations

ANG Air National Guard

bls below land surface

DA drainage area

DEM digital elevation model

DOD U.S. Department of Defense

DRBC Delaware River Basin Commission

EPA U.S. Environmental Protection Agency

HAGS Horsham Air Guard Station (Willow Grove)

HWSA Horsham Water and Sewer Authority

LHA EPA Long-term Health Advisory levels for PFOS and PFOA in drinking water

NASJRB Naval Air Station/Joint Reserve Base (Willow Grove)

Navy U.S. Navy 

NAWC Naval Air Warfare Center (Warminster)

NPDES National Pollution Discharge Elimination System

PADEP Pennsylvania Department of Environmental Protection

PFAS per- and polyfluoroalkyl substances

PFOA perfluorooctanoic acid

PFOS perfluorooctane sulfonate

ppm part per million

STP sewage treatment plant

USGS U.S. Geological Survey
Groundwater Withdrawals and Regional Flow Paths at 

and near Willow Grove and Warminster, Pennsylvania 

—Data Compilation and Preliminary Simulations for 

Conditions in 1999, 2010, 2013, 2016, and 2017

By Daniel J. Goode and Lisa A. Senior

Abstract

In 2014, groundwater samples from residential and public supply wells in the vicinity of 

two former U.S. Navy bases at Willow Grove and Warminster, and an active Air National Guard 

Station at Horsham, Bucks and Montgomery Counties, Pennsylvania, were found to have 

concentrations of perfluorooctanoic acid (PFOA) and perfluorooctane sulfonate (PFOS), which 

are per- and polyfluoroalkyl substances (PFAS), above U.S. Environmental Protection Agency 

(EPA) provisional health advisory (HA) levels for drinking water. Five supply wells near the 

bases were shut down because of PFAS contamination. In 2016, after EPA established a Lifetime 

HA for PFAS in drinking water that is lower than the provisional HA in place in 2014, at least 13 

additional supply wells near the bases were shut down because of PFAS contamination. At the 

request of the U.S. Navy, and in consultation with other Federal and State agencies and local 

stakeholders, the U.S. Geological Survey used historical and recent data on well withdrawals, 

recharge rates, aquifer properties, groundwater levels, and stream base flow to evaluate regional 

groundwater-flow paths from identified areas of PFAS groundwater contamination or potential 

PFAS sources at the bases. Groundwater withdrawals near the bases from public supply and 

other large wells decreased substantially from the 1990s to 2017, increasing the proportion of 

groundwater recharge that discharged to local streams. A preliminary groundwater-flow model, 

calibrated using 1,009 groundwater levels and 17 stream base flow estimates, simulated regional 

flow paths from the bases and showed that recharge at the bases discharged to withdrawal wells 

and local streams, generally within a mile or two of the bases. Supply and remediation wells at 

the bases captured some of the recharge on base areas of possible PFAS contamination, whereas

other base recharge was simulated to flow to nearby public supply wells and streams, depending 

on water use and aquifer recharge conditions between 1999 and 2017. The locations of many 

residential wells near the bases that were identified by the Navy and Air National Guard as 

having elevated PFAS concentrations were generally consistent with the simulated flow paths 

from possible sources at the bases. However, there are some areas of observed PFAS 

contamination where no flow paths from base sources were simulated. Additionally, no data 

were available on PFAS concentrations in groundwater in some areas of simulated flow paths 

from base sources. Data and models used for this study are provided in this report and in digital 

data releases to support further investigations and model revisions. 
2

Introduction

Groundwater is a substantial source of public, domestic, and industrial water supply in 

areas underlain by the Stockton Formation, a fractured sandstone aquifer, in southern 

Montgomery and Bucks Counties, Pennsylvania, where two formerly active military bases are 

located. Both bases were active for at least 50 years from the 1940s until they were closed at the 

recommendation of the Base Realignment and Closure (BRAC) commission. The Willow Grove 

Naval Air Station Joint Reserve Base (NASJRB) in Horsham Township, Montgomery County, 

Pennsylvania, was operated by the U.S. Navy (Navy) from 1942 until September 2011, with a 

portion currently (2019) operated by the Pennsylvania Air National Guard (ANG) and designated 

Horsham Air Guard Station (HAGS) (fig. 1). At present (2019), the entire Willow Grove 

NASJRB base property remains within the control of the Navy or the ANG, with the exception 

of a 3-acre parcel transferred to the Federal Aviation Administration. The Naval Air Warfare 

Center (NAWC) Warminster (formerly the Naval Air Development Center, Johnsville) in 

Warminster and Northampton Townships, Bucks County (fig. 1) was active during 1944–96. 

Since 1996, all but about 5 acres of the former NAWC Warminster property has been transferred 

from the Navy to local municipalities, Bucks County, or private owners. 

Previous investigations identified areas of groundwater contaminated by volatile organic 

compounds at both Willow Grove NASJRB and NAWC Warminster, which subsequently were 

addressed by remediation. Ongoing remediation at present (2019) includes bioremediation at 

Willow Grove NASJRB (Naval Facilities Engineering Command, 2018) and pump and treat 

technology at NAWC Warminster (Battelle Memorial Institute, 2016).

In summer 2014, through sampling of public supply wells under the third Unregulated 

Contaminant Monitoring Rule (UCMR3) (U.S. Environmental Protection Agency, 2012) 

groundwater at and near the former Willow Grove NASJRB and former NAWC Warminster was 

found to be contaminated with perfluorooctane sulfonate (PFOS) and perfluorooctanoic acid 

(PFOA), compounds that may pose a risk to human health if present above certain concentrations 

in drinking water. In 2014, the U.S. Environmental Protection Agency (EPA) provisional health 

advisory (PHA) levels were 0.2 microgram per liter (μg/L) for PFOS and 0.4 μg/L for PFOA in 

drinking water (U.S. Environmental Protection Agency, 2014). PFOS and PFOA are part of a 

group of compounds more broadly designated as per- and polyfluoroalkyl substances (PFAS). 

Five public supply wells in Horsham and Warrington Townships were shut down in 2014 after 

UCMR3 sampling showed that produced water had PFAS concentrations that exceeded the PHA 

levels.
3 

Figure 1. Location of former Willow Grove Naval Air Station Joint Reserve Base (NASJRB), former Naval Air Warfare Center Warminster (NAWC), and 

active Horsham Air Guard Station (HAGS), land-surface altitudes, streams, and model grid boundary, Bucks and Montgomery Counties, Pennsylvania. 
4

When PFOS and PFOA were identified in samples from public supply wells in 2014, 

production wells on and near the Willow Grove NASJRB and former NAWC Warminster were 

the primary sources of public water supply for communities within and surrounding these 

facilities and private domestic wells supplied many nearby residences. Potential sources of 

PFOS and PFOA in the area include fire-suppressant compounds (fluorinated surfactants in 

aqueous film forming foams) used on and (or) near these facilities when the bases were 

operating. Some public supply wells have been shut down since summer 2014 because PFOS 

and PFOA concentrations in produced water were higher than the PHA levels of 0.2 μg/L for 

PFOS and 0.4 μg/L for PFOA in drinking water (U.S. Environmental Protection Agency, 2014). 

In May 2016, the EPA revised the PFOS and PFOA health advisories to lower concentration 

levels and established a Lifetime Health Advisory (LHA) not to exceed 70 nanograms per liter 

(ng/L) or 70 parts per trillion (ppt) (equivalent to 0.07 μg/L) for combined concentrations of 

PFOS and PFOA (U.S. Environmental Protection Agency, 2016). Groundwater from public 

supply or private wells was still a primary source of drinking water in the area in 2014, but after 

May 2016, additional (at least 13) public supply wells were shut down in Horsham, Warrington, 

and Warminster Townships because of health concerns related to the presence of PFOS and 

PFOA. Since 2016, some supply wells have remained active, or resumed active status, with 

treatment that was subsequently installed to remove PFOS and PFOA from pumped 

groundwater, as needed.

After PFAS was discovered in the area, the Navy and the ANG offered nearby residences 

with private domestic wells that yielded water with PFOS and PFOA concentrations greater than 

the LHA to be connected to public drinking water supplies. The Navy and ANG also have 

established a program to monitor PFOS and PFOA concentrations in water from nearby 

residential wells that have not been connected to public supply; these monitored wells include 

those for which combined PFOS plus PFOA concentrations were less than the LHA of 0.07 μg/L 

but higher than 0.04 μg/L.

Management and mitigation of groundwater that is contaminated with PFOS and PFOA 

on and near the Willow Grove NASJRB and the nearby former NAWC Warminster requires 

understanding of groundwater flow directions under various pumping scenarios. A regional 

understanding of groundwater flow is needed as a result of extensive pumping near and between 

the two facilities, located about 3 miles (mi) apart (fig. 1). Because of the complex aquifer 

heterogeneity, the distribution and complex pattern of pumping wells, variability of recharge,

and groundwater-surface water exchanges in the area, a tool such as numerical groundwater 

modeling is essential to interpret site monitoring data, and to understand and describe 

groundwater flow directions and associated contaminant transport under various scenarios. Such 

a tool also would be useful to inform selection of additional monitoring locations and could be 

updated as new monitoring data become available. 

Flowing groundwater can carry contaminants from source areas near the land surface to 

downgradient discharges at streams or wells by advection, by which dissolved contaminants such 

as PFAS are carried with the water. However, additional processes can affect transport of 

contaminants to varying extents, and these additional processes are not addressed by flow 

modeling alone. Diffusion is small-scale molecular spreading from areas of high concentration to 

areas of low concentration. Some contaminants will additionally adsorb to aquifer materials, 

which causes an apparent delay, or retardation, of the rate of contaminant migration relative to 

the groundwater. 
5

Contaminant transport in fractured-rock aquifers is further affected by the variability of 

the properties of fractured and unfractured strata. Transport occurs primarily in open fractures or 

fracture zones, and the variability of the apertures of fractures causes variability in local 

velocities that results in dispersion, or the separation of contaminants on flow paths having 

different velocities. Further, contaminants can diffuse from fast flowing fractures into immobile 

or nearly immobile water within the matrix of pores in unfractured rock. This matrix diffusion 

can further retard the rate of transport of contaminants relative to the groundwater in the 

fractures. 

Previous Investigations

In studies related to remedial investigations by the Navy and their contractors at the two 

former Willow Grove and Warminster bases, the U.S. Geological Survey (USGS) collected 

geophysical logs (Conger 1997, 1998, 1999; Conger and Bird, 1999; Sloto, 2008) conducted 

aquifer-interval isolation (packer) testing of wells (Sloto, 1997, 1998, 2007; Sloto and others, 

2002) and hydrologic investigations (Sloto and others, 1998; Sloto, 2002), and prepared watertable maps at or near the bases (Sloto and others, 2001; Sloto and Grazul, 1995). Sloto and others 

(1996) and Sloto (2010) investigated groundwater flow and volatile organic compound 

contaminant migration at a Superfund site near NAWC Warminster. The USGS also conducted 

township-wide hydrogeologic studies for Warminster Township (Sloto and Davis, 1982) and 

Warwick Township (Bird, 1998) and prepared a water-table map for Warwick Township 

(Rowland, 1997). Sloto and others (1996) described the hydrology of the Stockton Formation in 

Hatboro and Warminster Township. The hydrogeology of northern Bucks County east of 

Neshaminy Creek, which includes an area of diabase, was described by Sloto and Schreffler 

(1994). 

Additional investigations have been conducted by local water suppliers, by regulatory 

agencies, and by parties responsible for contamination of groundwater in the area. Many of these 

investigations are described in documents available in public record depositories. Specific 

investigations that provided data for this study are cited throughout this report. 

Purpose and Scope

This report documents hydrologic data compiled to support development and calibration 

of a numerical model to simulate regional groundwater flow under selected conditions during 

1999–2017 at and near the former Willow Grove NASJRB and NAWC Warminster in 

Montgomery and Bucks Counties, Pennsylvania. The hydrologic data include water levels in 

wells, streamflow, point discharges to streams, groundwater withdrawals (pumping wells and 

quarries), and precipitation. The report also describes model development, calibration, and 

simulations. The model simulations provide information on regional groundwater-flow paths and 

groundwater/surface-water exchanges, in relation to changes in water supply pumping and 

recharge rates. Groundwater-flow paths from possible PFAS sources at the former bases are 

simulated under steady-state conditions for the years 1999, 2010, 2013, 2016, and 2017. 

Limitations of the study findings and modeling are also discussed.
6 

Description of Study Area

The study area lies within Bucks and Montgomery Counties in southeastern 

Pennsylvania, where land use is mixed. A small part of the study area is in Philadelphia County. 

The largest land uses in the study area in 2010 were residential (about 36 percent), wooded/open/

vacant (about 26 percent), industrial/commercial/institutional/transportation (about 24 percent), 

and agricultural (about 12 percent) (Delaware Valley Regional Planning Commission, 2013; 

2019). The study area has been settled since the 18th century and has become increasingly more 

populated and urbanized since the 1960s. The former Willow Grove NASJRB is in Horsham 

Township, Montgomery County, and the former NAWC Warminster is in Warminster and 

Northampton Townships, Bucks County (fig. 1). Other municipalities adjacent to, or near (within 

2 mi of), the former bases include Warrington, Warwick, and Upper Southampton Townships, 

and Ivyland Borough in Bucks County, and Upper Moreland Township and Hatboro Borough in 

Montgomery County (fig. 1). Reported 2010 populations for these municipalities are 26,147 in 

Horsham Township, 32,682 in Warminster Township, 39,726 in Northampton Township, 23,418 

in Warrington Township, 14,437 in Warwick Township, 15,152 in Upper Southampton 

Township, 1,041 in Ivyland Borough, 24,015 in Upper Moreland Township and 7,360 in Hatboro 

Borough (U.S. Census Bureau, 2019). 

Geologic Setting

The study area lies within the Piedmont physiographic province, with most of the study 

area in the Gettysburg-Newark Lowlands Section (fig. 2) characterized by rolling topography 

and shallow valleys (Sevon, 2000). The southern portion of the study area is located in the 

Piedmont Upland and Lowland Sections. The Gettysburg-Newark Lowlands are underlain by 

sedimentary rocks of the Newark Basin, which were deposited in a rift basin during the Triassic

(about 250 million years ago) and subsequently buried, compacted, and faulted, and in places, 

intruded by diabase (Lyttle and Epstein, 1987). Sedimentary formations in the Newark Basin, in

order of decreasing age, are the Stockton Formation, Lockatong Formation, and Brunswick 

Group. The Piedmont Upland and Lowland Sections are underlain by Paleozoic metasedimentary

rocks and older (Proterozoic) metamorphic rocks.

The central part of the study area is underlain by the Triassic Stockton Formation (fig. 3), 

which consists of gray to reddish brown sandstones and conglomerates, with siltstone and shale. 

The Stockton Formation was deposited unconformably on folded and faulted metamorphic rocks 

of Paleozoic- and Proterozoic-age, which crop out along the southern border of the Stockton 

Formation and include carbonates, quartzites, schists, and gneisses (fig. 3). The Stockton 

Formation has been divided into three members in Montgomery County, which are present in a 

generally fining upward sequence: the lower member is composed of conglomerates and 

sandstones; the middle member is composed of fine- to medium-grained sandstones; and the 

upper member is composed of shales, siltstones, and fine-grained sandstones (Rima and others, 

1962). Locally, lithologies may interfinger in the Stockton Formation and beds may pinch out 

(Rima and others, 1962) or be laterally discontinuous, likely a result of the fluvial or deltaic

origin of some deposits within the Stockton Formation (Turner-Peterson and Smoot, 1985). 
7

Figure 2. Physiographic provinces of Pennsylvania and physiographic-province sections in the study 

area in southeastern Pennsylvania.

The Lockatong Formation, which overlies and is conformable with the Stockton 

Formation, consists of fine-grained lacustrine deposits (siltstones and shales) and is relatively 

resistant to erosion, forming ridges that rise above flat or rolling topography underlain by the 

Stockton Formation to the south and rocks of the overlying Brunswick Group (shales, siltstones, 

and fine-grained sandstones) to the north (fig. 3). 

Bedding within the Stockton Formation in southeastern Montgomery and Bucks Counties 

generally strikes northeast or east-northeast and dips about 5 to 18 degrees to the northwest or 

north-northwest in this region, with an average dip of about 12 degrees (Rima and others, 1962). 

Although bedding in the Stockton Formation may be laterally discontinuous in places owing to 

interbedding of lithologies, the beds have been mapped as oriented with the overall general 

northwest-dipping structure (Rima and others, 1962). In addition to the Stockton Formation, the 

other Triassic sedimentary rocks in the study area also generally strike northeast and dip 

northwest, except where interrupted by faulting (Schlische, 1992). Both major and minor faults 

are present in the area. A major regional fault, the Chalfont Fault, is about 3 mi (5 kilometers 

[km]) north of Willow Grove NASJRB (fig. 3), and a zone of intense normal faulting south of 

the Chalfont Fault has been mapped about 1 mi (1.5 km) west of the Willow Grove NASJRB 

(Schlische, 1992). A north-south trending, steeply dipping diabase dike is mapped as intruding 

the Stockton Formation on the east border of the former Willow Grove NASJRB (Rima and 

others, 1962; Lyttle and Epstein, 1987). Other smaller dikes may be present, as indicated in 

mapping by Lyttle and Epstein (1987) that shows a minor dike about 1 mi or less west of the 

north-south trending dike on the eastern border of Willow NASJRB. 

The depth to competent bedrock generally is about 40 feet (ft) or less but varies 

depending on lithology and topographic setting (Low and others, 2002). The competent bedrock 

is overlain by weathered bedrock (saprolite) and soil formed in-situ and, in a few places in the 

southern part of the study area, by thin deposits of Tertiary and Cretaceous unconsolidated 

sediments (fig. 3). 
8 

Figure 3. Geology of the study area and major faults, streams, selected U.S. Geological Survey stream gages and observation wells, and 

location of model grid, Bucks and Montgomery Counties, Pennsylvania. Geology and faults from Berg and others (1980). See figure 1 for names 

of municipalities.
9 

Figure 3. Geology of the study area and major faults, streams, selected U.S. Geological Survey stream gages and 

observation wells, and location of model grid, Bucks and Montgomery Counties, Pennsylvania. Geology and faults 

from Berg and others (1980).—Continued
10

Hydrologic Setting

The Triassic sedimentary geologic units that underlie the study area form fractured-rock 

aquifers. Recharge to the fractured-rock aquifers occurs through the overlying soil and weathered 

rock and groundwater flows through a network of fractures both parallel and orthogonal to 

bedding, which commonly results in apparent preferential flow and permeability in the strike 

direction, parallel to bedding (Longwill and Wood, 1965). Anisotropy is the characteristic of 

permeability, or hydraulic conductivity, varying with direction. Lateral and vertical changes in 

lithology in the Stockton Formation affect the water-bearing properties of the aquifer, as finergrained deposits (shales, siltstones) commonly have lower permeability than coarser-grained 

deposits (sandstones and conglomerates) in the formation (Rima and others, 1962). The 

hydrogeology of the Stockton Formation and other geologic units of the study area have been 

described in more detail by Rima and others (1962), Longwill and Wood (1965), Greenman 

(1955), and Newport (1971). Low and others (2002) provide an overview of geohydrologic 

properties of the Stockton and other formations from well records and previous investigations in 

southeastern Pennsylvania. 

Although precipitation commonly is distributed relatively evenly throughout the year, 

groundwater recharge varies seasonally, with the lowest recharge rates occurring when 

evapotranspiration rates are highest in late summer to fall, and the highest recharge rates 

occurring from winter to spring, as reflected in annual fluctuations in water levels. For example, 

in the USGS observation well BK-1020, located in the study area on NAWC Warminster (fig. 3), 

long-term (1975–2018) daily mean depth to water is greatest in the fall months and least in 

spring months (fig. 4), reflecting relatively lower and higher seasonal net recharge rates, 

respectively. Groundwater discharges locally to streams or to pumping wells. Streams draining 

the study area include the Little Neshaminy Creek, Neshaminy Creek, Pennypack Creek, 

Wissahickon Creek, and their tributaries (figs. 1 and 3). The area does not have glacial deposits, 

and no valley fill has been mapped. Streams are riffle-pool type, with some reaches underlain by 

exposed bedrock. Shallow stream fluvial deposits underlie some streams in the study area. 

Losing reaches have been identified in Little Neshaminy Creek during low-flow late summer 

months (Sloto and Davis, 1982). The former Willow Grove and Warminster bases lie on high 

ground (fig. 1) that forms topographic divides between stream basins. 

Groundwater discharge rates to streams that drain areas underlain by the dipping beds of 

the Newark Basin may be partly controlled by geologic structure and anisotropy in hydraulic 

conductivity; in similar settings, groundwater discharge to streams oriented parallel to strike may 

be relatively less than discharge to streams oriented parallel to the dip of beds (Fan and others, 

2007). 
11

Figure 4. Long-term (1976–2018) daily mean depth to water in U.S. Geological Survey observation well 

BK-1020, completed in the Stockton Formation, Warminster Township, Bucks County, Pennsylvania. See 

figure 3 for location of well BK-1020.

Hydrologic Conditions and Water Use, 1999–2017

Hydrologic data, including precipitation, streamflow, and groundwater levels, were 

compiled for periods selected for groundwater model calibration (1999, 2010, 2013, 2016) and 

for 2017. These five periods were selected based on data availability to represent changing 

pumping and climatic conditions before and after PFAS contamination in groundwater was 

found in 2014. An extensive dataset for groundwater withdrawals (U.S. Geological Survey, 

1998) and spatially distributed groundwater levels for water-table maps (Sloto and others, 2001; 

Rowland, 1997) for the mid- to late 1990s were available to characterize conditions for 1999, the 

period with the largest groundwater withdrawals in the area among time period considered 

(1999–2017). Year 2010 was selected because Willow Grove NASJRB, supplied by two 

production wells on base, was still active. Year 2013 was selected to represent conditions after 

Willow Grove NASJRB closed, but before nearby production wells were shut down owing to 

PFAS contamination in 2014. Year 2016 was selected to represent conditions when numerous 

production wells in the area were (or had been since 2014) shut down as result of PFAS 

contamination. Year 2017 was selected to represent conditions when numerous production wells 

in the area were shut down since 2014 or 2016 but other wells had been put back into production 

with treatment to remove PFAS. 

Residents, institutions, industries, and commercial activities in the study area are served 

by public supply water systems and, in some areas, rely on individual wells for water supply. 

Some public supply systems import surface water from sources outside the study area. 
12

Wastewater in much of the study area is collected by sewers and directed to a few treatment 

plants that locally discharge treated wastewater to streams. Reported water use data include 

groundwater withdrawals and point-source discharges to streams.

Precipitation 

Annual county-wide precipitation values, as determined by the National Oceanic and 

Atmospheric Administration (NOAA) (2018a,b) for Bucks and Montgomery Counties during the 

period 1990–2017 (table 1), were similar, with annual totals for Montgomery County slightly 

less (about 0.8 inch [in.]) than the totals for Bucks County. The 1990–2017 mean annual 

precipitation values of 47.56 and 46.79 in. for Bucks and Montgomery Counties, respectively, 

were about 2 in. (or about 4 percent) greater than the long-term (1901–2000) values of 45.50 and 

44.82 in., respectively. 

During 1990–2017, annual average precipitation for the two counties was more than 5 in. 

(or 10 percent) below long-term (1901–2000) mean values in 1992, 1997, 2001, and 2016 

(relatively dry years) and more than 5 in. above long-term mean values in 1993, 1996, 2003, 

2004, 2006, 2009, and 2011 (relatively wet years) (table 1; fig. 5). The two driest years were 

2001 and 2016, with average values for the two counties of 37.35 and 37.55 in., about 8 in. (or 

about 18 percent) less than the two-county average long-term mean of 45.16 in. The two wettest 

years were 1996 and 2011, with average values for the two counties of 62.53 and 64.49 in., about 

17 and 19 in. (or about 38 and 42 percent, respectively) greater than the two-county average 

long-term mean of 45.16 in. 

Figure 5. Mean annual total precipitation for Bucks and Montgomery Counties, Pennsylvania, 

1990–2017. Data from the National Oceanic and Atmospheric Administration (2018a,b). 
13

Table 1. Annual total county-wide precipitation for Bucks and Montgomery Counties, Pennsylvania, 

1990–2017. Data from the National Oceanic and Atmospheric Administration (2018a,b).

Annual precipitation, in inches

Year Bucks County Montgomery County Average of two counties

1990 49.17 46.90 48.04

1991 41.81 41.46 41.64

1992 40.97 38.83 39.90

1993 51.14 50.11 50.63

1994 49.39 48.22 48.81

1995 41.61 41.07 41.34

1996 63.00 62.05 62.53

1997 40.32 37.80 39.06

1998 41.55 40.78 41.17

1999 48.48 48.45 48.47

2000 45.45 46.14 45.80

2001 38.40 36.30 37.35

2002 44.84 44.42 44.63

2003 58.15 59.02 58.59

2004 52.33 53.87 53.10

2005 48.92 46.49 47.71

2006 53.85 53.78 53.82

2007 48.79 46.74 47.77

2008 48.11 47.13 47.62

2009 51.81 52.00 51.91

2010 44.58 45.38 44.98

2011 66.06 62.91 64.49

2012 42.08 42.11 42.10

2013 48.06 48.06 48.06

2014 47.47 47.85 47.66

2015 43.41 43.65 43.53

2016 37.81 37.28 37.55

2017 44.22 41.19 42.71

1990–2017 

mean 47.56 46.79 47.17

1901–2000 

mean 45.50 44.82 45.16
14 

Streamflow

Streamflow has been measured at USGS stream gages at selected sites within the model

area for various periods of time (table 2). During part or all of the 1999–2017 period, stream 

gages within the model area were and (or) are active at one site on Little Neshaminy Creek (gage 

01464907 at Valley Road), two sites on Neshaminy Creek (gages 01464750 near Rushland and 

01465200 near Penns Park), four sites on Pennypack Creek (gages 01467031 at Horsham, 

01467036 at Hatboro, 01467039 at Willow Grove, and 014670413 at Bethayres), one site on a 

tributary to Mill Creek (gage 01465460 Iron Works Creek at Richboro), and one site on 

Wissahickon Creek (gage 01473900 at Fort Washington) (fig. 6). Additionally, other USGS 

stream gages active during the 1999–2017 period at the boundary of, or just outside of, the model 

area are located farther downstream on Pennypack Creek (gage 01467042 at Pine Road) and both

upstream and farther downstream on Neshaminy Creek (gage 01464720 North Branch at 

Chalfont and gage 01465500 at Langhorne, respectively) (fig. 6). Two other stream gages on 

Little Neshaminy (gage 01464930 Unnamed tributary to Little Neshaminy Creek at Traymore 

and gage 01464984 Little Neshaminy Creek at Walton Road) were active for short periods until

the mid-1990s (table 2, fig. 6).

The base-flow component of streams reflects the amount of groundwater discharge to

streams and can be used to estimate recharge. Base flow at selected stream gages was estimated 

through a hydrograph separation method using daily mean streamflow and the local-minimum 

option of the computer program HYSEP (Sloto and Crouse, 1986) as coded in the USGS 

Groundwater Toolbox (Barlow and others, 2017). Daily mean streamflow data at the USGS 

stream gages were obtained from the USGS National Water Information System (NWIS) 

database (U.S. Geological Survey, 2018). Streamflow is separated into base-flow and runoff 

components, as shown, for example, for monthly values at gage 01464907 Little Neshaminy 

Creek at Valley Road (fig. 7A,B). In southeastern Pennsylvania, natural base flow typically 

varies seasonally and in relation to recharge. 
15 

Figure 6. Location of selected U.S. Geological Survey stream gages and long-term observation wells, wells with water levels used in 

model calibration, selected point-source discharges (2010 rates) to streams above and near gages where data were used in model calibration, 

and model grid. Geology and faults from Berg and others (1980) as described in figure 3.
16 

Table 2. Period of record and type of data recorded for selected U.S. Geological Survey stream gages and observation wells in the study area in and near 

Willow Grove and Warminster, Pennsylvania, 1990–2017. 

[USGS, U.S. Geological Survey; mi2, square miles; STP, sewage treatment plant; ft, feet; --, no data or not applicable]

USGS site 

identification 

number

USGS site name Draninage 

area, in mi2

Continuous

record

Data used 

for model 

calibration

Begin date End date Notes

Stream gages

Discontinued in 2015 or before

01467031

Pennypack Creek at 

Horsham 3.48 Yes Yes 3/11/2010 10/7/2015

01467036

Pennypack Creek at 

Hatboro 4.36 Yes Yes 4/2/2010 10/6/2015

01467039

Pennypack Creek near 

Willow Grove 22.2 Yes Yes 2/24/2010 10/9/2015 Below STP

014670413

Pennypack Creek at 

Bethayres 35.3 Yes Yes 12/5/2009 10/13/2015 Below STP

01464930

Unnamed Trib to Little 

Neshaminy Creek at 

Traymore 4.34 Yes No 1/1/1994 12/31/1995

01464984

Little Neshaminy Creek 

at Walton Rd nr 

Jacsksonville 40.06 Yes No 11/9/1985 9/30/1992

Active as of 2019

01473900

Wissahickon Creek at 

Fort Washington 40.8 Yes Yes 9/1/1961 Current Below STP Upper Gwynedd, Ambler

01467042

Pennypack Creek at Pine 

Road 37.9 Yes Yes1 9/8/1964 Current Below STP

01464907

Little Neshaminy Creek 

at Valley Road 26.8 Yes Yes 11/25/1998 Current Below STP

01464720

NB Neshaminy Creek at 

Chalfont 31.5 Yes No 12/19/1990 Current

01464750

Neshaminy Creek near 

Rushland 91 Yes No 12/5/1986 Current Diversions

01465500

Neshaminy Creek near 

Langhorne 210 Yes No 10/1/1934 Current

Gage moved downstream to Brownville Rd for 

record from 10/25/2014 to 8/5/2016

Peak flow and partial continuous record

01465200

Neshaminy Creek near 

Penns Park 170 Yes No 10/1/2001 Current

Stage recorder and crest-stage gage until 2018, 

then converted to continuous record

01465460

Iron Works Creek at 

Richboro 3.69 No Yes 8/23/1981 Current Wire-weight gage; stream-discharge rating
17 

Table 2. Period of record and type of data recorded for selected U.S. Geological Survey stream gages and observation wells in the study area in and near 

Willow Grove and Warminster, Pennsylvania, 1990–2017.—Continued

[USGS, U.S. Geological Survey; mi2, square miles; STP, sewage treatment plant; ft, feet; --, no data or not applicable]

USGS site 

identification 

number

USGS site name

Well 

depth, 

in ft

Type of record

Data used 

for model 

calibration

Begin date End date Notes

Observation wells

401157075032001 BK-1020 395 Continuous Yes 9/1/1975 Current Well completed in Stockton Fomation

400808075210401 MG-225 486 Continuous Yes 9/13/1956 Current Well completed in Stockton Formation

Monthly levels - Montgomery County network

401233075102701 MG-2222 185 Monthly Yes 7/1/2005 Current Montgomery County well network

400554075042201 MG-2225 -- Monthly Yes 7/1/2005 Current Montgomery County well network

401410075142601 MG-2233 -- Monthly Yes 7/1/2006 Current Montgomery County well network

400743075191001 MG-2234 -- Monthly Yes 8/1/2006 3/1/2009 Montgomery County well network

400833075033501 MG-2236 -- Monthly Yes 8/1/2006 Current Montgomery County well network

401029075161601 MG-2240 -- Monthly Yes 6/1/2009 Current Montgomery County well network

1

Stream flow data for gage 01467042 Pennypack Creek at Pine Road used to estimate missing 2010 data for Pennypack Creek gages 01467031, 01467036, and 01467039 for 

base-flow separation computations
18 

Figure 7. Data collected at gage 01464907 Little Neshaminy Creek at Valley Road and Pennypack Creek 

gages, where base flow was determined using the HYSEP (local-minimum option) hydrograph separation 

method from the U.S. Geological Survey Groundwater Toolbox (Barlow and others, 2017): A, Daily base 

flow and total streamflow, 1999–2018; B, Monthly base-flow and runoff components of streamflow, 2014–

17; C, monthly base flow and sum of reported monthly discharges from sewage treatment plants (STPs) 

above the gage, 2014–17; and D, base flow at gages 01467031 Pennypack Creek at Horsham, 01467036 

Pennypack Creek tributary at Hatboro, and 01467039 Pennypack Creek at Willow Grove with reported 

discharges from National Pollution Discharge Elimination System (NPDES) sources above Hatboro gage 

and Upper Moreland Hatboro STP above Willow Grove gage, 2010–16. Reported monthly STP discharge 

flows from the Pennsylvania Department of Environmental Protection (2019a).
19 

Figure 7. Data collected at gage 01464907 Little Neshaminy Creek at Valley Road and Pennypack Creek 

gages, where base flow was determined using the HYSEP (local-minimum option) hydrograph separation 

method from the U.S. Geological Survey Groundwater Toolbox (Barlow and others, 2017): A, Daily base 

flow and total streamflow, 1999–2018; B, Monthly base-flow and runoff components of streamflow, 2014–

17; C, monthly base flow and sum of reported monthly discharges from sewage treatment plants (STPs) 

above the gage, 2014–17; and D, base flow at gages 01467031 Pennypack Creek at Horsham, 01467036 

Pennypack Creek tributary at Hatboro, and 01467039 Pennypack Creek at Willow Grove with reported 

discharges from National Pollution Discharge Elimination System (NPDES) source above Hatboro gage 

and Upper Moreland Hatboro STP above Willow Grove gage, 2010–16. Reported monthly STP discharge 

flows from the Pennsylvania Department of Environmental Protection (2019a). —Continued 

Measured flow at the streams in the study area may include point-source discharges from 

sewage treatment plants (STPs), quarries, or other activities. These discharges are discussed in a 

separate section “Discharges to Streams.” Daily streamflow values used in hydrograph 

separation were not adjusted for discharges, and thus the estimated base flow may include 

discharges that are relatively steady. Apparent base flow may be affected by dams, the largest of 

which in the study area is the Bradford reservoir (drainage area [DA] 11 square miles [mi

2

]) 

upstream of gage 01464907 Little Neshaminy Creek at Valley Road (DA 26.8 mi2

); the presence 

of the dam reduces peak flows and extends post-storm recession, which may result in higher 

amounts of post-storm apparent base flow for a few days, although the magnitude of this possible 

effect was not determined. Measured base flow may be reduced by groundwater withdrawals 

that capture some of the recharge that would otherwise contribute to base flow. Monthly baseflow values determined from daily streamflow at gage 01464907 Little Neshaminy Creek at 

Valley Road reflect seasonal fluctuations in addition to a component of flow contributed by one 

or more STPs upstream of the stream gage (fig. 7C), which is equal to base flow during some 

low-flow conditions common in late summer/early fall. Extreme low flows were measured in 

summer 1999 at gage 01464907, but base flows subsequently rose to the near long-term average 

of about 18.8 cubic feet per second (ft3/s) (median 15.6 ft3/s) (fig. 7A) following recharge 

associated with the passage of Hurricane Floyd in mid-September 1999. Base flow at gage 

01464907 Little Neshaminy Creek at Valley Road also may be affected (decreased) by nearby 

production well pumping, as indicated by losing reaches identified in 1980 near this stream gage 

location when nearby production wells were pumping (Sloto and Davis, 1982, fig. 15, p. 27). 
20

Seasonal patterns in monthly average base flow also are indicated for 2010–15 data 

collected at Pennypack Creek stream gages and the Upper Moreland Hatboro STP (fig. 7D). The 

reported monthly average discharge from the Upper Moreland Hatboro STP is about 40 percent 

of the monthly base flow at gage 01467039 Willow Grove downstream. The larger amount of 

base flow relative to drainage area at gage 01457039 Willow Grove (DA 22.2 mi2

) compared to 

base flow at upstream gages 01467031 at Horsham (DA 3.48 mi2) and 01467036 at Hatboro (DA 

4.36 mi2

) partly reflects the contribution of the Upper Moreland Hatboro STP to flow at the 

Willow Grove stream gage. Additionally, base flow at the Pennypack Creek stream gages at 

Hatboro and Horsham may be reduced by nearby pumping, and at times is very small (less than 

0.2 ft3

/s). Base flow at gage 01467036 Hatboro during late summer 2014 was nearly equal to the 

reported monthly average discharge from a National Pollution Discharge Elimination System 

(NPDES) source (identified as “ABB”) (fig. 7D) to a tributary above that stream gage (fig. 6). 

Annual base flow was computed from daily streamflow data collected at one stream gage 

on Little Neshaminy Creek, five on Pennypack Creek, one on Iron Works Creek, and one on 

Wissahickon Creek (table 3, fig. 3; see appendix 1 for computation details). A few months of 

data for early 2010 and late 2015 were missing for four Pennypack Creek stream gages; these 

missing data were estimated using the complete record from downstream gage 01467042 

Pennypack Creek at Pine Road (table 3; appendix 1, tables 1.2 and 1.4). Additionally, base flow 

at the Iron Works Creek gage was estimated from correlation of discrete flow measurements at 

the Iron Works Creek gage with daily mean flow and monthly base flow at the Pennypack Creek 

at Horsham gage (appendix 1, tables 1.4 and 1.5). The annual base flow is reported in cubic feet 

per second and, as normalized to drainage area, in inches (table 3). For a 19-year period of record 

from 1999 through 2017 at gage 01464907 Little Neshaminy Creek at Valley Road, base flow 

ranged from 5.47 to 13.65 in., with a median of 8.61 in.; for this same period of record, base flow 

ranged from about 26 to 45 percent of total streamflow with a median of about 36 percent 

(appendix 1, table 1.1). Of the years considered for modeling, base flow was lowest in 1999 and 

greatest in 2013; base flow in 2016 was slightly (about 5 percent) less than base flow in 2010 

(table 3).

Similar to flow at gage 01464907 Little Neshaminy Creek at Valley Road, streamflow at 

the Pennypack Creek stream gages may be affected by nearby pumping and (or) STP discharges, 

including a relatively large STP discharge in Upper Moreland Township above the three most 

downstream Pennypack Creek stream gages (fig. 6). Annual base flow for 2010–15 at stream 

gages on Pennypack Creek below the Upper Moreland STP discharge (gages at 01467039 at 

Willow Grove, 014670413 at Bethayres, 01467042 at Pine Road) is from 5 to 9 in. greater than 

the two stream gages above this STP discharge (gages 01467031 at Horsham and 01467036 

Hatboro) (table 3), reflecting the contribution of STP discharge to base flow. Base flow at gage 

01473900 Wissahickon Creek at Fort Washington likely is also affected by STP discharge above 

the stream gage.
21 

Table 3. Annual average base flow at selected stream gages near Willow Grove and Warminster, Pennsylvania, 1999 and 2010–2017. Annual average base flow 

for calendar year computed using daily values and HYSEP local-minimum hydrograph separation method in the U.S. Geological Survey Groundwater Toolbox 

(Barlow and others, 2017). Estimated values (see appendix 1 for computations) in italics. 

[BF, base flow; NAD83, North American Datum 1983; lsd_ft (NAVD88), land surface datum in feet above North American Vertical Datum 1988; DA (mi2), drainage area in 

square miles; --, no data]

USGS site 

number USGS local site name Latitude (NAD83)

Longitude 

(NAD83)

lsd_ft 

(NAVD88)

DA 

(mi2) 1999 2010 2011 2012 2013 2014 2015 2016 2017

Average base flow for calendar year in cubic feet per second

01464907 Little Neshaminy Creek at Valley Road 

near Neshaminy 40.22917 -75.11972 190 26.8 10.95 15.94 25.15 13.23 17.12 18.81 16.44 15.14 16.45

01465460 Iron Works Creek at Richboro1 40.19833 -75.01083 174 3.69 -- 2.24 3.84 1.24 1.94 1.83 1.43 -- --

01467031 Pennypack Creek at Horsham2,3 40.17722 -75.12194 226 3.48 -- 2.51 4.52 1.36 2.36 2.41 1.88 -- --

01467036 Pennypack Creek tributary at Hatboro2,3 40.18083 -75.11222 219 4.36 -- 2.34 4.43 1.17 1.95 1.89 1.4 -- --

01467039 Pennypack Creek near Willow Grove2,3 40.15139 -75.07583 163 22.2 -- 24.21 35.2 18.3 21.73 23.61 18.67 -- --

014670413 Pennypack Creek at Bethayres3 40.11278 -75.07306 94 35.3 -- 37.54 43.36 29.2 33.51 37.28 26.02 -- --

01467042 Pennypack Creek at Pine Road 40.08972 -75.06917 79 37.9 -- 41.03 50.98 31.34 35.29 39.12 28.95 25.28 --

01473900 Wissahickon Creek at Fort Washington 40.12389 -75.22000 139 40.8 -- 37.82 49.62 30.87 38.83 41.78 29.16 27.23 --

Average base flow for calendar year in inches

01464907 Little Neshaminy Creek at Valley Road 

near Neshaminy 40.22917 -75.11972 190 26.8 5.55 8.07 12.74 6.70 8.67 9.53 8.33 7.67 8.33

01465460 Iron Works Creek at Richboro1 40.19833 -75.01083 174 3.69 -- 8.24 14.13 4.56 7.14 6.73 5.26 -- --

01467031 Pennypack Creek at Horsham2,3 40.17722 -75.12194 226 3.48 -- 9.79 17.63 5.30 9.21 9.40 7.33 -- --

01467036 Pennypack Creek tributary at Hatboro2,3 40.18083 -75.11222 219 4.36 -- 7.29 13.79 3.64 6.07 5.88 4.36 -- --

01467039 Pennypack Creek near Willow Grove2,3 40.15139 -75.07583 163 22.2 -- 14.80 21.52 11.19 13.29 14.44 11.42 -- --

014670413 Pennypack Creek at Bethayres3 40.11278 -75.07306 94 35.3 -- 14.44 16.67 11.23 12.89 14.34 10.01 -- --

01467042 Pennypack Creek at Pine Road 40.08972 -75.06917 79 37.9 -- 14.70 18.26 11.22 12.64 14.01 10.37 9.05 --

01473900 Wissahickon Creek at Fort Washington 40.12389 -75.22000 139 40.8 -- 12.58 16.51 10.27 12.92 13.90 9.70 9.06 --

1 Iron Works base flow estimated from discrete measurements at site correlated with daily values at Horsham; then applying that correlation to base flow at Horsham. 2 2010 partial record for Pennypack at Horsham, Hatboro, Willow Grove; missing Jan-Mar except Willow Grove missing Jan-Feb; full calendar year base flow (BF) for these sites 

estimated from correlation with 2010 Pine Rd monthly BF values (cubic feet per second). 3 2015 partial record for Pennypack at Horsham, Hatboro, Willow Grove, Bethayres; missing Oct-Dec; full calendar year base flow (BF) for these sites estimated from correlation 

with 2015 Pine Rd monthly BF values (cubic feet per second). 
22

Groundwater Levels

Wells with water-level measurements included residential wells, monitoring wells, and 

unused production wells. In the study area, residential and production wells typically are drilled 

into bedrock, with the casing extended to the depth of competent bedrock, and the remaining 

depth of the hole left open. Water-producing fractures may be present throughout the open 

interval. Monitoring wells are initially drilled as open holes but are commonly reconstructed to 

be open to only relatively short intervals of about 20 ft of aquifer thickness. Monitoring wells are 

generally shallower in depth than residential or production wells. 

Data on groundwater levels in the study area for the period of interest include one-time 

USGS measurements in wells that were used to prepare water-table maps for Warwick Township

during 1994–95 (Rowland, 1997), at and near Willow Grove during October 1999 (Sloto and 

others, 2001), and in Upper Gwynedd in 2000 (Senior and Ruddy, 2004). Other data include 

groundwater levels in wells inventoried by the USGS for various studies, monthly measurements 

in selected wells in Montgomery County (Montgomery County Department of Health and 

Human Services, Horsham Water and Sewer Authority [HWSA]), and periodic measurements in 

monitoring wells at the former NAWC Warminster (U.S. Navy and their consultants), Willow 

Grove NASJRB (U.S. Navy and consultants), HAGS (Air National Guard and consultants) and 

Raymark Superfund site (EPA and Pennsylvania Department of Environmental Protection 

[PADEP]) in Hatboro. Sources of well-construction and water-level data are identified by well in 

a data release (Senior and others, 2020). The data were compiled by the USGS during 2017–19 

from the USGS NWIS database (U.S. Geological Survey, 2018), previous reports (H&S 

Environmental, Inc., 2012; CDM Federal Programs Corporation, 2003; Gilmore and Associates, 

Inc., 2015; Koman Government Solutions, LLC, 2018a; Leidos, 2018; Resolution Consultants, 

2019), or from written communications (Toby Kessler, Gilmore Associates, Inc. for HWSA data; 

Andrew Barton, Battelle Memorial Institute, for 2015 data in Area B of NAWC Warminster; 

Colin Wade, PADEP, for water-level data at Raymark).

Additionally, the USGS collects continuous data in two long-term observation wells in 

the study area, both completed in the Stockton Formation: well BK-1020 on the east edge of the 

former NAWC Warminster base and well MG-225 in the southwest edge of the study area (figs. 

3 and 6). Monthly water levels for wells in the Montgomery County network are available in the 

USGS NWIS database and five of these wells are in the study area (fig. 6). Water-level data are 

plotted in hydrographs of groundwater levels in observation well BK-1020 (fig. 8A,B) and water 

levels in five Montgomery County network wells (fig. 8B). Among the wells with continuous or 

monthly data, annual seasonal fluctuations in groundwater levels are greatest in well BK-1020, 

ranging from 10 to 15 ft, but the pattern of seasonal fluctuations in groundwater levels is similar 

for all the wells (fig. 8B). Mean annual water-level values for well BK-1020, MG-225, and the 

Montgomery County network wells were used for model calibration datasets, as listed in Senior 

and others (2020). For other wells with more than one measurement per year, the average water 

level was used for that year model calibration.
23 

Figure 8. Groundwater levels in well BK-1020 and five Montgomery County network wells, as well as groundwater-level data used for model 

calibration in the time periods selected to represent calendar years A, 1999, and B, 2010, 2013, and 2016 conditions. Montgomery County network 

wells include MG-2222, MG-2233, MG-2234, MG-2236, and MG-2240. 
24

Groundwater-level data collected by the USGS during the 1990s have a larger spatial 

distribution in the study area than subsequent data as a result of USGS data collection for countywide assessments, water-table maps, and other studies. To represent conditions in 1999, data 

from 1990 to 2000 were compiled, with the objective of obtaining broad spatial distribution while 

minimizing periods of extreme high or low seasonal values indicated by water levels in long-term 

observation well BK-1020 (fig. 8A). Fewer and (or) less spatially distributed water level-data 

were available for years 2010, 2013, and 2016. To represent the 2016 conditions, water levels 

measured during 2015–17 were compiled, including levels measured in monitoring wells on the

former Willow Grove NASJRB, NAWC Warminster, and HAGS from 2015–17 (fig. 8B). Mean 

water levels in well BK-1020 were lower in 1999 and 2016 than in 2010 and 2013 (fig. 8A,B), 

consistent with lower base flow (and inferred recharge) for those years, respectively.

Summary statistics for a total of 670 wells with water levels used to characterize 

conditions for the 1999, 2010, 2013, and 2016 model scenarios show that the median depth to the 

top of the open interval was about 52 ft below land surface (bls) and the median depth to the 

bottom of the open interval was about 90 ft bls for wells with known or reported well construction 

information (558 wells had a known depth to the top of the interval and 602 wells had a known 

depth to the bottom of the interval) (table 4). For wells with missing construction information, the 

depth of the open interval was estimated as described in Senior and others (2020). If not known, 

the depth to the bottom of the open interval was estimated to be 150 ft bls. If not known, the depth 

to the top of the open interval was estimated to be 40 ft bls for wells greater than 80 ft in depth, 30 

ft bls for wells 65 to 79 ft in depth, 20 ft for wells 40 to 64 ft in depth, 15 ft for wells 30 to 39 ft in 

depth, 10 ft for wells 19 to 29 ft in depth, and about half the depth to the bottom of the open 

interval for wells less than 19 ft in depth. Statistics for open intervals of all wells, including 

estimated values for missing construction data, used for model input are similar to statistics for 

open intervals for wells without estimated values (table 4). Water levels ranged from about 200 ft 

below land surface to 23 ft above land surface (negative values, such as -23 ft, indicate the water 

level is above land surface). The median depth to water ranged from about 15.7 to 18.5 ft bls. 

Depth to water commonly varies by topographic setting, with the depth to water greater at hilltops 

than in valley settings (Low and others, 2002). Depth to water also may be affected by nearby 

pumping, such as in the cones of depression shown near pumping production wells in the 

Warwick Township potentiometric-surface map (Rowland, 1997).
25 

Table 4. Summary statistics for wells with water levels used in model calibration scenarios for 1999, 2010, 

2013, and 2016 conditions, including land-surface altitude, depth to top and bottom of well opening, date of

water-level measurement, and depth to water near Willow Grove and Warminster, Pennsylvania. 

[LSA_ft (NAVD88), land-surface altitude in feet above North American vertical datum of 1988, Top Open Depth 

(ft), depth to top of well opening in feet below land surface; Bottom open depth (ft), depth to top of well opening in 

feet below land surface; Date_WL, date of water-level measurement; DTW below MP (ft), depth to water, in feet 

below measuring point; MPA_ft, measuring point altitude, in feet above NAVD88; --, no data] 

Summary 

statistics

LSA_ft 

(NAVD88)

Top open 

depth (ft)

Bottom 

open depth 

(ft)

Date_WL

DTW 

below MP 

(ft)

MPA_ft 

(NAVD88)

No estimated values for well construction included in statistics

All wells with water levels used for 1999, 2010, 2013 and 2016 conditions, n= 670

Count, n= 670 558 602 -- -- --

Minimum 116 4 6 -- -- --

10th percentile 249 14 31 -- -- --

Mean 307 74 138 -- -- --

Median 313 52 90 -- -- --

90th percentile 357 169 330 -- -- --

Maximum 455 360 2148 -- -- --

Estimated values for well construction included in statistics

All wells with water levels used for 1999, 2010, 2013 and 2016 conditions, n= 670

Minimum 116 3 6 -- -- --

10th percentile 249 15 32 -- -- --

Mean 307 68 139 -- -- --

Median 313 40 100 -- -- --

90th percentile 357 157 311 -- -- --

Maximum 455 360 2148 -- -- --

Wells with water levels used for 1999 conditions, n= 368

Minimum 116 3 6 4/16/1990 -12.46 116

10th percentile 224 18 38 11/5/1993 5.00 225

Mean 303 56 158 2/23/1997 27.17 304

Median 313 40 135 5/8/1998 17.05 313

90th percentile 359 115 332 11/15/1999 68.84 359

Maximum 444 278 2148 12/11/2001 198.68 459

Wells with water levels used for 2010 conditions, n= 123

Minimum 164 6 35 1/1/2010 -18.15 164

10th percentile 285 18 45 11/2/2010 5.98 285

Mean 307 57 114 10/28/2010 18.15 307

Median 312 45 86 11/2/2010 16.59 312

90th percentile 328 109 200 11/4/2010 28.62 328

Maximum 455 189 580 12/14/2010 103.60 455
26 

Table 4. Summary statistics for wells with water levels used in model calibration scenarios for

1999, 2010, 2013, and 2016 conditions, including land-surface altitude, depth to top and bottom of well

opening, date of water-level measurement, and depth to water near Willow Grove and Warminster, 

Pennsylvania.—Continued 

[LSA_ft (NAVD88), land-surface altitude in feet above North American vertical datum of 1988, Top Open Depth 

(ft), depth to top of well opening in feet below land surface; Bottom Open Depth (ft), depth to top of well opening 

in feet below land surface; Date_WL, date of water-level measurement; DTW below MP (ft), depth to water, in feet 

below measuring point; MPA_ft, measuring point altitude, in feet above NAVD88] 

Summary 

statistics

LSA_ft 

(NAVD88)

Top open 

depth (ft)

Bottom 

open depth 

(ft)

Date_WL

DTW 

below MP 

(ft)

MPA_ft 

(NAVD88)

Estimated values for well construction included in statistics

Wells with water levels used for 2013 conditions, n= 121

Minimum 164 14 38 1/1/2013 -16.40 164

10th percentile 235 20 46 4/30/2013 6.15 235

Mean 305 61 115 10/16/2013 20.90 305

Median 311 47 87 11/14/2013 18.48 311

90th percentile 328 122 200 11/14/2013 34.47 328

Maximum 357 190 580 12/12/2013 102.00 357

Wells with water levels used for 2016 conditions, n= 397

Minimum 164 4 10 4/16/2015 -23.20 164

10th percentile 272 11 26 5/16/2016 2.33 273

Mean 312 78 109 9/1/2016 16.46 313

Median 313 52 80 5/24/2016 15.79 314

90th percentile 352 186 242 7/24/2017 30.83 354

Maximum 455 360 580 7/24/2017 101.63 455
27 

Groundwater Withdrawals 

Groundwater has been a major source of supply for public drinking water, commercial, 

and industrial uses in the study area, with withdrawal rates increasing as the area started to

urbanize in the 1950s (Sloto and Davis,1982). Additionally, pumping for remediation of 

groundwater contamination started in the 1990s at various locations, including at the Raymark 

Superfund site in Hatboro (U.S. Environmental Protection Agency, 2006; CDM Federal Programs

Corporation, 2003), at two industrial properties (CRC and the former Fisher Porter, now ABB) in

Warminster, and at the former NAWC Warminster. Other groundwater withdrawals in the study 

area are associated with quarry dewatering operations; most quarries in the study area are in the 

Lockatong Formation or Paleozoic carbonate rocks (some of which are indicated by areas of 

negative values of mapped land-surface altitude as much as -275 ft [-84 m] relative to the North 

American Vertical Datum of 1988 [NAVD88] in fig. 1). 

The study area is within the Delaware River Basin Commission (DRBC) Southeastern 

Pennsylvania Ground Water Protected Area (GWPA), where regulations that affect groundwater 

withdrawals greater than 10,000 gallons per day (gal/d) have been in effect since 1980. The 

GWPA includes all of Montgomery County, the southwestern half of Bucks County, and smaller 

parts of adjacent counties and was established by the DRBC to address concerns that groundwater

withdrawals were resulting in undesirable lower water levels and base flow (Delaware River

Basin Commission, 1980). 

Data on groundwater withdrawals in the study area were compiled from various sources 

including the DRBC, Pennsylvania Department of Environmental Protection (PADEP), and water 

suppliers. Annual average withdrawals by wells and quarries and sources of withdrawal data are 

described in more detail in the data release for this report (Senior and others, 2020). Withdrawal 

rates were compiled by the USGS for a water-use assessment (Solley and others, 1998) using 

available data in 1996 and documented by the USGS (U.S. Geological Survey, 1998) for the 

DRBC (Delaware River Basin Commission, 2019) to assist in a 1998 assessment of water

resources. Those withdrawals were used by this study for 1999 conditions. Withdrawal rates 

reported to PADEP (2019b) were used to represent the 2010, 2013, 2016, and 2017 conditions 

and are listed with well construction information by Senior and others (2020). Additional 

information about withdrawal rates for groundwater remediation extraction wells at NAWC

Warminster were compiled from various reports prepared for the Navy (Koman Government 

Solutions, LLC, 2018a,b,c) and for extraction wells at Raymark from the EPA (2006). The

compiled data represent groundwater withdrawals reported to either the DRBC or PADEP and

may not include some withdrawals not reported to these agencies, as water use less than 10,000

gal/d is not required to be reported to the DRBC. For example, residential water use (estimated to 

be about 60 gal/d per person) and some withdrawals for golf course irrigation are not reported; 

however, in the case of residential wells, withdrawals are considered relatively small and any 

consumptive effect is negligible if the residence has on-lot wastewater disposal and withdrawals

are returned to the groundwater. Quarry withdrawals for dewatering purposes commonly are 

discharged to nearby streams. Complete data on withdrawal rates for quarries are not available 

from PADEP sources (Pennsylvania Department of Environmental Protection, 2019b,c) for the

periods of interest and were estimated in some cases, where data were missing as described in

Senior and others (2020). At least two of the quarries (Glasgow Ivy Rock and Glasgow 

Montgomery) in the model area currently (2019) are used for clean fill and thus presumed to be 

inactive, with little to no dewatering. Permitted values for quarry discharge (Michael Kutney,
28

Pennsylvania Department of Environmental Protection Bureau of District Mining Operations, 

written commun., 2018) were used to confirm magnitude of likely quarry withdrawals. 

Production wells typically are constructed as open holes in the study area, and are usually 

8 to 12 in. in diameter with depths of 300 to 400 ft but can be as deep as 600 to 700 ft. Remedial 

extraction wells commonly are 6 in. in diameter and relatively shallow (less than 200 ft deep). 

For all 306 pumping wells included in the data compilation, the median depth to the top of 

opening (equivalent to the depth of casing) is 50 ft and the median depth to the bottom of the well 

is 301 ft (table 5). The distribution of well depths in the compiled data shows that most wells are 

drilled to depths of 300 to 400 ft (fig. 9A). Pumping rates range widely in wells of various depths 

but tend to be lower in shallow wells commonly drilled for remediation purposes and greater in 

deeper wells drilled for production and supply (fig. 9B). 

Figure 9. A, Distribution of well depths and B, relation between well depths and reported pumping rates 

for withdrawal wells in 1999 and 2010 in the model area near Willow Grove and Warminster, Pennsylvania. 
29 

Table 5. Summary statistics for pumping wells with withdrawals used in model scenarios for 1999, 2010, 

2013, 2016 and 2017 conditions, including land-surface altitude and depth to top and bottom of well

opening, near Willow Grove and Warminster, Pennsylvania. 

[LSA_ft (NAVD88), land-surface altitude in feet above North American vertical datum of 1988; Top Open Depth 

(ft), depth to top of well opening in feet below land surface; Bottom Open Depth (ft), depth to top of well opening 

in feet below land surface] 

Summary statistics

LSA_ft 

(NAVD88)

Top open 

depth (ft)

Bottom open 

depth (ft)

No estimated values for well construction included in statistics

All pumping wells with withdrawals during one or more years 

during 1999, 2010, 2013, 2016, and 2017

Count, n= 306 204 237

Minimum 59 10 50

10th percentile 179 19 81

Mean 260 53 347

Median 259 50 351

90th percentile 336 80 593

Maximum 499 424 763

Estimated values for well construction included in statistics

All pumping wells with withdrawals during one or more years 

during 1999, 2010, 2013, 2016, and 2017

Count, n= 306 306 306

Minimum 59 10 50

10th percentile 179 21 90

Mean 260 52 343

Median 259 50 301

90th percentile 336 70 550

Maximum 499 424 763
30

In the study area, reported annual withdrawals for supply purposes have declined since the 

1990s from rates of about 20 million gallons per day (Mgal/d) to about 10.5 Mgal/d in 2017 (table 

6; fig. 10). Total reported withdrawals for quarry dewatering in the study area has remained 

relatively constant during 1995–2017 at rates of about 9 to 10 Mgal/d, but pumping has varied by 

quarry. Although groundwater has been, and continues to be, a major source of drinking water 

supply in the area, some municipalities and water purveyors began to import surface water from 

outside the immediate area starting in the late 1990s, which accounts for some of the reported 

decline in groundwater withdrawals since 1999 (fig. 10).

Figure 10. Total annual average reported well and quarry withdrawals in the model area near 

Willow Grove and Warminster, Pennsylvania, in 1995, 1999, 2010, 2013, 2014, 2016, and 2017.

Data not compiled for years without plotted values.
31 

Table 6. Annual average withdrawals for wells and quarries in model area for selected years, 1995–2017, near Willow Grove and Warminster, 

Pennsylvania. 

Annual average

withdrawal rate, in 

million gallons per day

Annual average well withdrawal rate, in gallons per minute

Year Wells1 Quarries Number of 

wells Minimum 10th percentile Mean Median 90th percentile Maximum

1995 21.30 10.07 220 0.19 0.29 22.5 67.2 143.3 1,104.2

1999 20.27 10.07 244 0.19 0.30 13.9 57.7 130.2 1,104.2

2010 15.34 9.58 148 0.19 3.76 38.1 72.0 152.0 1,063.6

2013 14.21 9.68 149 0.05 2.59 33.2 66.2 163.6 618.6

2014 14.10 9.89 149 0.01 2.39 25.8 65.7 148.7 850.7

2016 13.05 9.64 137 1.00 3.00 26.0 66.2 146.2 983.0

2017 10.55 9.10 133 0.01 1.60 12.8 55.1 127.7 935.3

1 Data for 1995 and 1999 included more wells with lower pumping rates than data for years 2010–17; Wells pumping less than 1 cubic meter per day (220 

gallons per day or 0.15276 gallon per minute) were excluded from 1995 and 1999 totals and data used for model scenario for 1999 conditions. 
32

Discharges to Streams

Permitted discharges to streams under the NPDES in the study area include treated 

effluent from STPs, treated water from wells pumped for groundwater remediation, and quarry 

pumpage (fig. 6). Water pumped from quarries commonly is discharged to nearby streams. Data 

for annual discharge volumes were compiled from PADEP (2019a), from the Montgomery 

County Planning Commission (1994) for the year 1993, Montgomery County Planning 

Commission sewage treatment facilities status report (Drew Shaw, Montgomery County 

Planning Commission, written commun., 2019) for the year 1998, and from water-use data 

compiled by the USGS for 1995–98 (Solley and others, 1998; U.S. Geological Survey, 1998). 

Other discharges included those from groundwater remediation at NAWC Warminster and 

Raymark, which were estimated from groundwater withdrawals (Koman Government Solutions

LLC, 2018a; U.S. Environmental Protection Agency, 2006). 

Comparison of the annual average flow rates for point-source discharges to streams (table 

7) to annual base-flow values (table 3) indicates that point-source discharges likely represent 20

percent or more of annual base flow at gage 01464907 Little Neshaminy Creek, and 40 percent

or more of annual base flow at gages 01467039 Pennypack Creek at Willow Grove and

01473900 Wissahickon Creek at Fort Washington.

Table 7. Annual average discharge rates for point source discharges to streams in model area for model 

scenarios for 1999, 2010, 2013, 2016, and 2017 conditions near Willow Grove and Warminster, 

Pennsylvania. Listed rates do not include discharges from quarries. Sources of data described in Senior 

and others (2020). Location of discharges shown in figure 6.

Stream receiving point-source discharge Annual average discharge rate

1999 2010 2013 2014 2016 2017

Cubic feet per second

Little Neshaminy Creek above gage 01464907 1.86 3.11 3.01 3.16 3.20 3.20

Little Neshaminy Creek below gage 01464907 8.37 7.72 9.88 10.62 9.26 10.52

Pennypack Creek 8.63 7.90 8.59 9.74 7.96 8.27

Wissahickon Creek 16.61 15.12 15.74 16.58 15.98 15.59

Millions of gallons per day

Little Neshaminy Creek above gage 01464907 1.20 2.01 1.94 2.04 2.07 2.07

Little Neshaminy Creek below gage 01464907 5.41 4.99 6.38 6.86 5.98 6.80

Pennypack Creek 5.58 5.11 5.55 6.30 5.14 5.35

Wissahickon Creek 10.73 9.77 10.17 10.71 10.33 10.07

Total of listed discharges 22.92 21.87 24.05 25.91 23.52 24.29
33

Simulation of Regional Groundwater Flow

Regional groundwater flow was simulated using a numerical model based on field 

measurements in the study area and on previous studies in nearby areas that have similar 

hydrogeologic conditions. Calibration of the model produces estimates of the large-scale 

hydrogeologic properties of the groundwater system and recharge rates. The model was used to 

characterize the average regional flow directions as affected by changes in withdrawal rates from 

public supply wells, and by changes in recharge. This preliminary model was made with 

numerous simplifying assumptions to reduce model complexity and the number of parameters, 

while capturing the major influence of the geology on regional groundwater-flow paths. The 

simulated groundwater-flow paths from possible PFAS source locations on the bases represent 

paths along which groundwater containing dissolved PFAS may have migrated. 

Conceptual Model of the Groundwater-Flow System

The dipping Triassic sedimentary rocks form a layered, leaky aquifer system, as 

described in other reports (Rima and others, 1962; Senior and Goode, 1999; Goode and Senior, 

2000). Below a shallow weathered zone, groundwater flow is through fractures in the rock. This 

fractured rock network is commonly characterized by lower storage and higher permeability than 

the overlying weathered rock. Directions of groundwater flow are controlled by both the 

hydraulic gradient and the orientation of fractures parallel and orthogonal to bedding. The 

structural control on groundwater flow is represented by specifying three-dimensional anisotropy 

in this preliminary model that parallels the strike and dip directions. 

Recharge is related to geology, topography and the extent of urbanization (increased 

impervious area and [or] reduced permeability of soils). Base flow in selected Pennsylvania 

streams was shown to be related to surficial bedrock geology, with areas underlain by Triassic 

rocks having relatively low base flow (and inferred recharge) compared to areas underlain by 

carbonate rocks and some other geologic units (White and Sloto, 1990; Sloto and Buxton, 2007). 

Therefore, in the current groundwater model, recharge was varied spatially based on the surficial 

bedrock geology. 

The permeability of the fractured-rock aquifers generally varies with depth, with lower 

permeability in soils and weathered rock than in unweathered fractured bedrock in the Triassic 

rocks. The depth of casing represents the thickness of weathered rock and soils, with casing 

typically extending about 10 ft into competent bedrock. The median and mean casing length 

reported for pumping wells is about 50 ft (table 4). Water-producing fractures in the Stockton 

and Lockatong Formations are most commonly reported at depths of 400 ft or less (Low and 

others, 2002; Bird, 1998), but water-bearing zones are reported at depths of 500 ft or more. In the 

conceptual model for this study, water-bearing properties of the unweathered bedrock decrease 

with depth below 400 ft. 

The groundwater flow system in the study area is relatively local, with groundwater 

divides generally near surface-water divides except in areas affected by pumping or in carbonate 

rock settings. Most groundwater discharge is to nearby streams, based on water-table mapping 

and streamflow seepage data. The study area boundaries were established at remote watershed 

divides and beneath remote major streams. 
34

Model Approach and Methods

The three-dimensional numerical model, MODFLOW 6 (Langevin and others, 2017, 

2018; Provost and others, 2017) was used to simulate steady-state flow for several time periods 

with different recharge rates and well withdrawals. Model input files were produced using 

ArcGIS for Desktop 10.7 (Esri, 2019) and ModelMuse version 4 (Winston, 2019), a graphical 

pre- and post-processor for groundwater modeling. PEST++ (Welter and others, 2015, White and 

others, 2019), which is based on the methods and algorithms of PEST (Doherty, 2015), was used 

for model calibration. The Python script for PEST++ calibration of the model using FloPy 

(Bakker and others, 2016, 2018) is provided in the model archive data release (Goode and 

Senior, 2020), along with all MODFLOW 6 input and output files. MODPATH 7 (Pollock, 2016, 

2017) was used to compute flow paths (particle tracks or pathlines). 

Features, assumptions, and approximations in the regional groundwater-flow model 

included 

• flow conditions in each simulation were approximated as steady state, which represent

annual average conditions for the associated time period;

• all model layers were assumed to be confined or semi-confined in that the horizontal

transmissivity did not depend on water level;

• recharge and hydraulic conductivity were uniform within zones defined by the regional

geology; and

• flow within the fractured-rock aquifers was simulated as bulk flow distributed throughout

the model cells, based on the equivalent porous media approximation.

Model Area and Spatial Discretization

The horizontal boundaries of the model area (fig. 1) were assumed to be no-flow 

boundaries. Several of the boundary segments were located to include regional streams around 

the Willow Grove and Warminster bases, based on the assumption that groundwater flow beyond 

these streams does not substantially influence flow paths near the bases. Additional no-flow 

horizontal boundaries were located along ridges (surface-drainage watershed boundaries) and 

geologic boundaries associated with low-permeability formations. These no-flow boundaries 

were specified to include all local streams and wells that were considered to potentially affect 

flow paths near the bases. 

The discretization by vertices (DISV) grid type with quad-tree refinement of 

MODFLOW 6 was selected to allow higher grid resolution near the bases and near wells, 

streams, and diabase dikes, while reducing the numerical computational burden in other parts of 

the model domain (figs. 11 and 12). The model coordinate system is in units of meters (m), and 

input data are in units of meters and days, for example cubic meters per day for withdrawal rates. 

The largest cells are 1,640 by 1,640 ft (500 by 500 m) and the smallest are 205 by 205 ft (62.5 by 

62.5 m). The DISV quad-tree refined grid type is compatible with ModelMuse 4 and 

MODPATH 7, the latter with the additional restriction that the grid not be rotated with respect to 

the coordinate system (Pollock, 2016). The horizontal grid configuration is the same for all 

layers of the model, which is a restriction for the DISV grid type. 
35

Figure 11. ModelMuse screenshot of full model grid with model streams near Willow Grove and 

Warminster, Pennsylvania. 
36

Figure 12. ModelMuse screenshot of detail of model grid cells and withdrawal well locations close to the 

bases near Willow Grove and Warminster, Pennsylvania. 

The altitude of the top of the model was interpolated to the grid cell nodes from the 3-m 

digital elevation model (DEM; Penn State Institutes of Energy and Environment, 2015). Seven 

model layers were used, each with uniform thickness (table 8, fig. 13). Layer thickness increases 

with depth, and the overall model thickness is about 800 ft (244 m). The top two layers represent 

overburden and weathered rock and are homogeneous. Layers 3–7 represent the unweathered 

fractured rock aquifer, in which properties were specified according to the geologic unit. The 

lower surface of the bottom model layer is a no-flow boundary. 
37

Figure 13. ModelMuse screenshot of details near Willow Grove of model grid in plan (top) and cross 

section (bottom) showing seven model layers with illustrative regional groundwater-flow paths from possible 

sources to discharge cells. The color of the pathlines change with altitude from blue (highest) to brown 

(middle) to black (lowest).
38

Table 8. Model layer thicknesses and hydraulic conductivity variation with depth

[K, hydraulic conductivity; Kmax, parameter for Kx in strike direction by hydrogeologic zone; m, meters; ft, feet]

Model 

layer Name Thickness (m)

Bottom depth 

(m)

Thickness 

(ft)

Bottom depth 

(ft)

K scaling with 

depth

Homogeneous overburden and weathered rock

1 Soil 6.1 6.1 20 20 Same K as layer 2

2 Weathered 9.1 15.2 30 50 Same K as layer 1

Unweathered rock, properties by hydrogeologic zone

3 Shallow 22.9 38.1 75 125 K = Kmax / 2

4 Transition 38.1 76.2 125 250 K = Kmax

5 Intermediate 45.7 121.9 150 400 K = Kmax

6 Deep 61.0 182.9 200 600 K = Kmax / 2

7 Basement 61.0 243.9 200 800 K = Kmax / 20
39

Model Hydrogeologic Zones

The model hydrogeologic properties and recharge rates are zoned according to the 

surficial bedrock geology. The Stockton Formation, which underlies most of the model area and 

the Willow Grove NASJRB and NAWC Warminster, is separated into upper, middle, and lower 

members in Montgomery County (Rima and others, 1962; Lyttle and Epstein, 1987). The 

members have different hydrogeologic characteristics; notably, the finer grained upper member 

is generally observed to yield less water than the middle and lower members (Rima and other, 

1962). In order to approximately account for this hydrogeologic variability, a hypothetical 

hydrogeologic contact between the upper and middle members is assumed in the model area in 

Bucks County (fig. 14). The location of the actual contact is unknown. For the model, the offset 

of the upper/middle contact from the contact between the Stockton and Lockatong Formations in 

Bucks County is conceptually extrapolated from its configuration in Montgomery County, and 

with the mapped diabase dike configuration. A hypothetical contact was also added in Bucks 

County between the middle and lower members of the Stockton Formation, but the 

hydrogeologic properties of the model are the same in these members, designated as the “middle 

and lower Stockton zone.” These hypothetical contacts in Bucks County are solely for the 

purpose of assigning hydrogeologic zones for this model. 

Figure 14. ModelMuse screenshot of model hydrogeologic zones: Upper Stockton, Middle and Lower 

Stockton, Lockatong, Brunswick, Carbonate, and Other. Hypothetical contacts between the upper and 

middle members and between the middle and lower members of the Stockton Formation were added in the 

Bucks County portion of the model. Bases shown in magenta outline.
40

Aquifer Inflows and Outflows

The principal inflows to the aquifer consist of local recharge from infiltration of 

precipitation and stream loss, and the principal outflows from the aquifer consists of discharges 

to wells and quarries (as withdrawals), and to streams (as base flow). Other processes, such as 

evapotranspiration, were not explicitly considered. 

Recharge 

Groundwater recharge from the land surface and unsaturated zone was simulated into the 

top layer (1) of the model. Recharge was zoned with the surface-mapped geology, based on 

findings from previous studies in the vicinity (White and Sloto, 1990; Schreffler, 1996; Sloto and 

Buxton, 2007) and the relative rates of recharge for each geologic unit were fixed based on those 

studies (fig. 15, appendix 2, table 2.1). For this preliminary model, a single scale factor 

parameter was used to adjust the initial recharge rates proportionally for the entire model domain

for all time periods. 

Figure 15. ModelMuse screenshot of model recharge rates in 2010 and 2013 in inches per year, zoned by 

geologic unit. Recharge rates for other periods were reduced to 90 percent of the values shown for 1999 

and 70 percent for 2016 and 2017.
41

Groundwater Withdrawals

Senior and others (2020) compiled groundwater withdrawal data and developed modelinput tables for use with the MODFLOW 6 WEL package. Groundwater withdrawal wells were 

assigned to the model layer corresponding to the mid-point of the open interval of the well. 

Wells for which the mid-point is below the model bottom were assigned to the bottom (7) model 

layer. 

The overall magnitude of the groundwater withdrawals used in the model are illustrated 

for 1999 (fig. 16) and 2017 (fig. 17) conditions using the well-withdrawal footprint visualization

method of Goode (2016) as implemented in ModelMuse (Winston and Goode, 2017). The 

withdrawal rate for each well is shown as a circular area centered on the well, and the area of the 

circle is proportional to the withdrawal rate. For closely spaced wells, a merged shape (footprint) 

is drawn that has the area corresponding to the combined withdrawals for all wells within that 

merged shape, based on a scale factor, the depth-rate-index (Goode, 2016). In this example, the 

scale factor is 10 inches per year (in/year) (0.25 meters per year [m/year]), comparable to the 

average recharge rate for the calibrated model. 

Figure 16. ModelMuse screenshot of model groundwater withdrawal rates visualized as well-withdrawal 

footprints (method of Winston and Goode, 2017) for simulation of conditions in 1999. The footprint depthrate-index (scale) is 10 inches per year. Bases shown in magenta outline.
42

Figure 17. ModelMuse screenshot of model groundwater withdrawal rates visualized as well-withdrawal 

footprints (method of Winston and Goode, 2017) for simulation of conditions in 2017. The footprint depthrate-index (scale) is 10 inches per year. Bases shown in magenta outline.
43

Increases in withdrawals from wells change groundwater levels and groundwater-flow 

paths, diverting paths from streams or other wells towards or into the pumped well. Past changes 

in pumping have thus changed past flow paths. As discussed previously, pumping near the 

Willow Grove and Warminster bases has generally declined from the 1990s until the present, as 

illustrated by the difference between the well-withdrawal footprints from 1999 to 2017 (fig. 18), 

and these changes in the model will change the simulated flow paths. 

Figure 18. ModelMuse screenshot of changes in groundwater withdrawal rates in the Willow Grove and 

Warminster area visualized by the difference between well-withdrawal footprints (method of Winston and 

Goode, 2017) for 1999 (fig. 16) and 2017 (fig. 17). Light gray area shows volumetric rate of withdrawals 

common to both years. Brown area shows decrease in volumetric rate from 1999 to 2017. Blue area shows 

increase in volumetric rate from 1999 to 2017. The footprint depth-rate-index (scale) is 10 inches per year.

Bases shown in magenta outline, and wells as black squares. 
44

Flow to and from Streams

Streams draining the study area include Park Creek, Little Neshaminy Creek, Neshaminy 

Creek, Pennypack Creek, Wissahickon Creek, Mill Creek, and their tributaries (fig. 1). The 

Schuylkill River drains the west edge of the model area. A stream network was developed for 

model input from the National Hydrography Dataset (U.S. Geological Survey, 2016) augmented 

with selected tributary streams (fig. 1) from datasets provided by Robert H. Keough (Bucks 

County Planning Commission, written commun., 2016), and William McLay (Montgomery 

County Planning Commission, written commun., 2016). 

Streamflow data collected at seven USGS stream gages were used to provide information 

about groundwater discharge; these gages included one on Little Neshaminy Creek (01464907 at 

Valley Road), four on Pennypack Creek (01467031 at Horsham, 01467036 at Hatboro, 

01467039 at Willow Grove, and 014670413 at Bethayres), one site on a tributary to Mill Creek 

(01465460 Iron Works Creek at Richboro), and one site on Wissahickon Creek (01473900 at 

Fort Washington) (table 2; figs. 3 and 6). 

As discussed in the Streamflow section, the base-flow component of streamflow reflects 

the amount of groundwater discharge to streams and can be used to estimate recharge to 

groundwater. Measured flow at the streams in the study area could also include upstream point 

discharges to those streams from STPs, quarries, or other activities. Streamflow data used to 

calculate base flow from groundwater were not adjusted for these point discharges; therefore, the 

estimated base flow could be overestimated. These point discharges are specified as inflows to 

the model, however, so that simulated flow at stream gage locations includes contributions of 

both base flow from groundwater and upstream point discharges. 

The MODFLOW 6 Streamflow Routing (SFR) package was selected to represent the

stream network because it accumulates flow from upstream reaches to downstream reaches, 

permitting the use of streamflow observations in model calibration, and because it limits stream 

loss (flow from the stream into the aquifer) to the amount of flow available in the simulated 

reach. Stream stages were specified in each reach using the 3-m DEM for simplicity and 

numerical stability, rather than computing for the stream stages from a hydraulic equation. 

Specified stream stages increase in the downstream direction in some simulated reaches as a 

result of the coarse model resolution and the approximate interpolation from the DEM, but this 

inconsistency does not substantially affect simulated regional groundwater levels and flow paths. 

However, the simulated groundwater-stream interactions in and near these model cells may not 

be representative of field conditions. 

Point discharges to streams (fig. 6) are represented as SFR reaches with specified 

INFLOW to that reach, which corresponds to the annual average volumetric discharge (in cubic 

meters per day [m3

/d]). A very small streambed conductance was specified in each reach with a 

point discharge (fig. 19) to minimize flow to or from the aquifer and to route the point discharge

flow to the downstream reach. 
45

Figure 19. ModelMuse screenshot of locations of model point discharges to streams. Bases shown in 

magenta outline.
46 

Bulk Hydraulic Conductivity and Anisotropy

The regional anisotropy in hydraulic conductivity (K) reflects the preferential 

groundwater flow along the strike of dipping beds, and reduced flow in the dip direction and 

perpendicular to those beds. Hydraulic conductivity was specified through a three-dimensional

tensor. The maximum hydraulic conductivity (Kx) is oriented in strike. The medium hydraulic 

conductivity (Ky) is oriented perpendicular to Kx in the plane of the bedding, and is a proportion 

of Kx, ranging from 1/10 to 1, referred to as horizontal anisotropy. The minimum hydraulic 

conductivity (Kz) is oriented orthogonal to the plane of Kx and Ky, or perpendicular to the 

bedding, and is a proportion of Kx ranging from 1/10 to 1, referred to as vertical anisotropy. 

Anisotropy in model layers 3–7, which represent unweathered fractured rock, was 

specified such that Kx, was rotated 20 degrees north from east-west, (Winston, 2019, fig. 7). 

The dip angle for Ky is variable, with a northwest down dip of 13 degrees specified in the 

Stockton Formation and 10 degrees in the Lockatong Formation and the Brunswick Group (fig. 

20) to reflect the measured orientation of sedimentary rock beds (Greenman, 1955; Rima and 

others, 1962; Longwill and Wood, 1965). 

The Carbonate and Other zones are assumed to be isotropic in the plane of the layer (Kx

= Ky) and to have a dip angle of zero and a vertical anisotropy (Kz/Kx) of 1/3. 

Figure 20. ModelMuse screenshot showing model down dip angles to the northwest for hydraulic 

conductivity tensor in model layers 3 to 7, which represent unweathered fractured rock. Bases shown in 

magenta outline.
47 

Model Calibration

Initial model simulations were conducted using starting values for aquifer properties and 

boundary conditions based on previous model studies in the Stockton Formation and nearby 

studies (Senior and Goode, 1999, 2013, 2017). Manual adjustments were made to qualitatively 

match the simulated water levels to available water-table maps, especially at Willow Grove 

(Sloto and others, 2001) and Warminster (Tetra Tech, written commun., 2019). 

Several model properties were selected for automatic calibration through nonlinear 

regression using PEST++ (White and others, 2019), including 

• bulk aquifer hydraulic conductivity values (K), including anisotropy values;

• recharge rate scale for four time periods; and

• streambed conductance value for all model streams.

A subset of the model K values was selected for automatic calibration, based in part on 

sensitivities calculated with PEST++, and based in part on the conceptual model, hydrogeologic 

judgement based on previous studies, and field testing (appendix 2, table 2.2). Some K values, 

such as Kz, oriented orthogonal to dipping beds, were adjusted proportionally during calibration

to the Kx value for that zone, and thus the anisotropy ratio is fixed. Other values, such as all K

values for the “Other” zone (fig. 14), were fixed at initial manually adjusted values. The

orientation of the K anisotropy tensor (Provost and others, 2017; Winston, 2019) was specified 

based on field measurements of geologic strata strike and dip and was not adjusted during 

calibration (fig. 20). 

The recharge rates for all four calibration time periods, 1999, 2010, 2013, and 2016, were 

adjusted proportionally to a single model parameter, RECH. Based on annual stream base-flow

analysis, and initial manual model adjustments, recharge values were 0.9 x RECH for 1999, 1.0 x 

RECH for 2010 and 2013, and 0.7 x RECH for 2017. Recharge values varied spatially by four 

hydrogeologic zones and were scaled proportionally to the single recharge parameter (fig. 15). 

To calibrate the model, 1,009 groundwater-level measurements in 670 wells were used

(fig. 8; table 4). Groundwater levels in 1999 are from measurements taken in early October 1999 

for a synoptic water-level map (Sloto and others, 2001) and other measurements in the 1990s 

included in the USGS National Water Information System (NWIS) database. Well construction 

information is limited for many of these wells. Measurements were assigned to specific model 

layers based on the mid-point of the open interval of the well. Average values of casing depth 

and total well depth were used to determine the appropriate model layer with missing 

construction data. More recent data in 2010, 2013, and 2016 are compiled from various sources 

including the USGS NWIS database and monitoring wells on and near the former bases. 

Streamflow data for calibration were taken from estimates of annual base flow at seven USGS 

stream gages in the study area (table 9). 

All water-level data were assigned a uniform weight of 1.0. The weight for all streamflow

data was 0.003 to account for the difference in the units of the streamflow data, m3/d compared 

to the water levels (m), and to ensure that the streamflow data contributed substantially to the 

calibration objective function in the regression (Doherty and Hunt, 2010). The stream flow

model error was about 14 percent of the total model error for the final calibration.
48 

Table 9. Streamflow measurements used for model calibration. 

[--, not available]

Site number Site name (all in Pennsylvania) Short name Streamflow, in cubic meters per day

1999 2010 2013 2016

01464907 Little Neshaminy Creek at Valley Road near Neshaminy LilNeshV 26,790 38,998 41,885 37,041

01465460 Iron Works Creek at Richboro IronWorks -- 5,480 4,746 --

01467031 Pennypack Creek at Horsham PennyHor -- 6,141 5,774 --

01467036 Pennypack Creek at Hatboro PennyHat -- 5,725 4,771 --

01467039 Pennypack Creek near Willow Grove PennyWG -- 59,232 53,164 --

014670413 Pennypack Creek at Bethayres PennyBeth -- 91,844 81,985 --

01473900 Wissahickon Creek at Fort Washington WissFtWash -- 92,529 95,001 66,620
49

The calibrated model is deemed adequate for simulation of regional groundwater flow 

and flow paths. Summary graphs of the simulated and measured water levels and base flow used 

for model calibration are presented in this section. More detailed results of the final calibration 

are presented in appendix 2, including the spatial distribution of water-level residuals for the four 

calibration periods.

The overall fit between the simulated and measured water levels for the four time periods 

used for calibration is shown in figure 21. The histogram of all water-level residuals (fig. 22) 

shows a relatively narrow concentration of residuals around zero, with about 75 percent of 

residuals smaller than plus or minus 16.4 ft (5 m). The mean residual is -0.79 ft (-0.24 m), and 

the root mean square residual is 0.6 ft (0.18 m). 

Figure 21. Simulated and measured water levels in wells for optimum parameters identified by calibration

for each model period, 1999, 2010, 2013, and 2016. 
50

Figure 22. Histogram of water-level residuals for all calibration data for all model periods. 
51

Streamflow is also relatively well represented by the model (fig. 23), given the limited 

number of parameters used, and other model limitations. 

Figure 23. Simulated and measured streamflow for optimum parameters identified by calibration for all 

model periods. 
52 

Effects of Recharge and Withdrawals on Groundwater-Flow Paths

MODPATH 7 (Pollock, 2016) was used to delineate groundwater-flow paths to discharge 

areas under steady-state flow for conditions in 1999, 2010, 2013, 2016, and 2017. The 2017 

simulation uses the same recharge rate as the 2016 simulation, but well withdrawals are reduced

further (table 6; Senior and others, 2020). For simulations of recharge discharging to wells or

streams, 1 particle was used in each model cell; in these simulations, wells pumping at small 

rates may not capture any particles.

The overall water budget of the steady-state model consists of a balance between the 

inflows to and outflows from the aquifer. Inflows are primarily recharge, with some stream loss

entering the aquifer, and outflows are composed of groundwater discharge to streams (base flow)

and well withdrawals (table 10). Areas of recharge that discharge either to wells or to streams

can be visualized through particle tracking using MODPATH. 

Table 10. Simulated groundwater budget. 

Water budget component Model area budget, in inches per year

1999 2010 2013 2016 2017

Inflow Recharge 9.98 11.09 11.09 7.76 7.76

Outflow Wells 2.75 2.25 2.16 2.06 1.77

Streams1 7.22 8.83 8.93 5.70 5.99

Water budget outflow Model area outflow, as percent of inflow

Wells 27.6 20.3 19.5 26.6 22.8

Streams1 72.4 79.7 80.5 73.4 77.2

1

Net outflow to streams = total outflow to streams minus inflow from streams

PFAS contamination in groundwater has been identified at and near the Willow Grove 

and Warminster bases. To help understand groundwater migration of contaminants, regional 

flow paths are delineated from identified possible locations of PFAS contamination to either 

streams or wells for the five time periods. The locations of PFAS sources on the bases are not 

well known, but were estimated from reported 2014–17 data on elevated PFAS concentrations in 

shallow monitoring wells, sites where aircraft crashed, and other information (Tetra Tech, 2014; 

Leidos, 2018; Resolution Consultants, 2019). Results for two time periods are presented in this 

section to show examples of changes in flow paths owing to changes in recharge and well 

withdrawals. Simulated water levels, areas contributing recharge to wells, and flow paths from 

sources for all time periods are presented in appendix 3. For simulations of groundwater flow 

paths from possible PFAS sources on the bases to wells or streams, nine particles were used in 

each model source cell.
53

1999 Conditions

The 1999 conditions (based on withdrawal rates compiled for the 1990s) include the 

largest well withdrawals, so a higher proportion of recharge is captured by wells rather than 

streams compared with other time periods. 

Simulated water levels (fig. 24) reflect surface topography and groundwater discharge to 

streams under average recharge conditions but are also influenced by drawdowns near large 

production wells. Both Willow Grove NASJRB and NAWC Warminster straddle groundwater 

level highs, corresponding to topographic ridges. Groundwater flow is perpendicular to the 

water-level contours, from high to low altitudes, with some deviation caused by the horizontal

anisotropy (Ky/Kx) within bedding planes and (Kz/Kx) across the dipping bedding planes. 

Figure 24. ModelMuse screenshot of simulated potentiometric surface in model layer 1 for 1999 

conditions. The contour interval is 20 feet. 
54

In 1999 the withdrawal rates were large relative to the annual average recharge, so a large 

portion of the model area is simulated to provide recharge to wells (fig. 25). Some wells are 

deriving recharge from adjacent topographic drainage basins, which indicates that the surfacewater basins and groundwater basins do not fully coincide. 

Figure 25. ModelMuse screenshot of simulated areas that contribute recharge to discharging wells

(green) and streams (white) for 1999 conditions. Bases shown in magenta outline.
55

Regional groundwater-flow paths are delineated from possible sources at Willow Grove 

and Warminster for 1999 steady-state flow conditions (fig. 26). Simulated results suggest that 

advective transport of contaminants under 1999 conditions is limited to a mile or two from the 

bases. 

Figure 26. ModelMuse screenshot of simulated regional groundwater-flow paths (black) from locations of 

possible per- and polyfluoroalkyl substances (PFAS) contamination in groundwater for 1999 conditions. 
56

Regional groundwater-flow paths from sources at Willow Grove discharge to streams and 

wells near the base (fig. 27). Flow paths to the northwest of the base are limited to the area 

between the base and Park Creek. Several flow paths discharge to large production wells near the 

base to the southeast and northeast. The base production wells at HAGS were simulated to 

capture flow paths from possible sources at the base under 1999 conditions. 

Figure 27. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 1999 

conditions near Willow Grove. Pink squares represent locations of possible per- and polyfluoroalkyl 

substances (PFAS) contamination in groundwater. Colors along paths range from green for relatively short 

travel times to magenta for relatively long travel times. 
57

Regional groundwater-flow paths from sources at Warminster discharge to streams and 

wells near the base (fig. 28). Flow paths from sources on the east side of the base discharge 

nearly due east to Iron Works Creek, part of the Neshaminy Creek drainage. Flow paths from the 

possible source on the base’s southern border discharge to Southampton Creek in the Pennypack 

Creek drainage. Flow paths from possible sources in the northwest part of the base discharge to 

local streams and to several production wells and remediation wells under 1999 conditions. 

Figure 28. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 1999 

conditions near Warminster, with withdrawal wells in black. Pink squares represent locations of possible

per- and polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from 

green for relatively short travel times to magenta for relatively long travel times.
58

2017 Conditions

The 2017 time period has the smallest groundwater withdrawals of the simulated

scenarios (tables 6 and 10; fig. 10). This time period was not used for model calibration. 

Under 2017 conditions, regional groundwater-flow paths from sources at Willow Grove 

discharge to streams and wells near the base (fig. 29). Flow paths to the northwest of the base are 

limited to the area between the base and Park Creek. Because fewer wells were pumping in 2017 

than in 1999, most of the flow paths from the north side of the base discharge to local streams 

instead of wells. With smaller withdrawals from nearby wells, a few flow paths flow deeper in 

the aquifer and discharge to streams farther away from the base than in 1999 conditions. 

Figure 29. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2017 

conditions near Willow Grove, with withdrawal wells in black. Pink squares represent areas of possible perand polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from green

for relatively short travel times to magenta for relatively long travel times. 
59

Regional groundwater-flow paths from sources at Warminster under 2017 conditions 

discharge to streams and wells near the base (fig. 30). As in 1999, flow paths from sources on the 

east side of the base discharge nearly due east to Iron Works Creek, part of the Neshaminy Creek 

drainage. Likewise, flow paths from the possible source on the base’s southern border discharge 

to Southampton Creek in the Pennypack drainage, although at lower altitudes in the stream, 

likely a result of reduced recharge in 2017 compared to 1999 conditions. Although remediation 

well pumping is similar, other well withdrawals are reduced in 2017 compared to 1999 

conditions. Paths from the possible source in the center of the Warminster base flow to the 

northwest, north, and northeast, while under the larger withdrawal 1999 conditions paths from 

that source flowed to the northwest. The 2017 flow paths are also deeper and discharge to more 

distant locations. 

Figure 30. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2017

conditions near Warminster, with withdrawal wells in black. Pink squares represent areas of possible perand polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from green

for relatively short travel times to magenta for relatively long travel times.
60

Groundwater-Flow Paths and PFAS Contamination in Wells

The regional flow paths are different for the different simulation periods owing to 

changes in model recharge rates and well withdrawals. A composite overlay of the simulated 

regional flow paths from possible PFAS sources for all simulated conditions provides a useful 

understanding of the potential locations of groundwater contaminated by possible PFAS sources 

at the bases. Such information could also help guide further delineations of contamination in the 

study area. 

The ANG and Navy began sampling groundwater from existing residential wells and new 

monitor wells drilled on the bases, and base flow in streams at and near the bases, after PFAS 

was first detected in supply wells in 2014. The ANG and Navy also sampled residential wells for 

PFAS compounds in areas surrounding the Willow Grove and Warminster bases in 2015. In 

cooperation with the Navy, EPA sampled residential wells near Willow Grove from 2014 

through mid-2017. Maps prepared by the Navy and the ANG and their contractors showing the 

extent of private-well sampling areas and levels of PFOA and PFOS concentrations in the wellwater samples relative to drinking-water LHAs indicate most exceedances occur near the bases 

(figs. 31 and 32, U.S. Navy, written commun., 2019). 

The paths from the steady-state models do not represent the actual advective transport of 

contaminants in part because the aquifer system is not at steady state, and a moving contaminant 

will change direction as flow conditions change in response to long-term, seasonal, and shorterterm changes to recharge and pumping. For example, if a contaminant entered the aquifer in 

1999 and did not discharge from the aquifer until 2017, it would be affected by all the conditions 

represented in the simulations above, at different points along its path. 

A composite visualization of the paths for all simulated conditions was constructed to 

approximate the areas of the aquifer that may have been impacted by possible PFAS sources. 

Each set of paths for a simulated year was colored gray and made partly transparent, and all sets 

for different years were overlaid (figs. 33–35). The composite paths are darkest where paths from 

most or all the simulations overlap, light gray paths indicate areas where fewer simulations 

overlap, and the lightest gray paths are those for a single simulation, where no paths were 

simulated during other simulations. 

Many of the off-base areas where PFAS concentrations in private residential wells exceed 

drinking-water LHA levels are on or near simulated composite regional groundwater-flow paths 

from possible PFAS contaminations sources at HAGS and NASJRB Willow Grove and NAWC 

Warminster (figs. 33–35). The simulated composite regional groundwater-flow paths from 

potential PFAS contaminations sources are also consistent with the locations of numerous 

pumping wells where PFAS concentrations have exceeded LHA levels (not shown). Although 

evaluation of detailed groundwater/surface-water interaction is beyond the scope of this report, 

results from stream sampling near the Willow Grove base to date (Leidos, 2018; Resolution 

Consultants, 2019) are generally consistent with areas where simulated groundwater-flow paths 

discharge to streams.
61

Figure 31. Private-well sampling areas and levels of perfluorooctane sulfonate (PFOS) and perfluorooctanoic acid 

(PFOA) in well-water samples relative to drinking-water health advisories for Horsham Air Guard Station and the 

former Naval Air Station Joint Reserve Base Willow Grove as of September 2019 (U.S. Navy, 2019a; PPT, parts per 

trillion). 
62 

Figure 32. Private-well sampling areas and levels of perfluorooctane sulfonate (PFOS) and perfluorooctanoic acid (PFOA) in well-water samples 

relative to drinking-water health advisories for the former Naval Air Warfare Center Warminster, Bucks County, Pennsylvania as of January 2019 

(U.S. Navy, 2019b; PPT, parts per trillion).
63 

Figure 33. Composite simulated regional groundwater-flow paths from possible per- and polyfluoroalkyl substances (PFAS) contamination sources at Willow

Grove and Warminster for conditions in 1999, 2010, 2013, 2016, and 2017 and maximum levels of perfluorooctane sulfonate (PFOS) and perfluorooctanoic acid 

(PFOA) in 2014–2018 well-water samples relative to drinking-water health advisories (see figs. 31 and 32 for symbol color explanation). Water-quality data from 

the U.S. Navy (written commun., 2019). 
64 

Figure 34. Detail at Willow Grove of composite simulated regional groundwater-flow paths from possible

per- and polyfluoroalkyl substances (PFAS) contamination sources for conditions in 1999, 2010, 2013, 

2016, and 2017 and maximum levels of perfluorooctane sulfonate (PFOS) and perfluorooctanoic acid 

(PFOA) in 2014–2018 well-water samples relative to drinking-water health advisories (see figs. 31 and 32 

for symbol color explanation). Water-quality data from the U.S. Navy (written commun., 2019). 
65 

Figure 35. Detail at Warminster of composite simulated regional groundwater-flow paths from possible

per- and polyfluoroalkyl substances (PFAS) contamination sources for conditions in 1999, 2010, 2013, 

2016, and 2017 and maximum levels of perfluorooctane sulfonate (PFOS) and perfluorooctanoic acid 

(PFOA) in 2014–2018 well-water samples relative to drinking-water health advisories (see figs. 31 and 32 

for symbol color explanation). Water-quality data from the U.S. Navy (written commun., 2019). 

Some of the areas where well concentrations have exceeded drinking-water LHA levels, 

however, lie beyond the extent of the composite simulated regional groundwater-flow paths from 

potential PFAS sources at HAGS, Willow Grove NASJRB and NAWC Warminster. These areas 

include a cluster of wells on the north side of Little Neshaminy Creek and north of the Willow

Grove base and a cluster of wells just southwest of the Willow Grove base (figs. 33, 34). There

are also wells as much as about 2 mi southwest and northwest of the NAWC Warminster base 

(figs. 33, 35). Possible explanations for the lack of alignment between simulated regional flow

paths and observed, relatively elevated PFAS concentrations include model error and (or) the

presence of unidentified PFAS sources, including sources unrelated to Navy or ANG activities, 

or mechanisms of PFAS transport.
66

Stream Loss to Groundwater

When groundwater levels near streams are below the stream stage, surface water from the 

stream infiltrates to groundwater. This process is enhanced where pumping wells near streams 

are withdrawing groundwater at a high rate, creating groundwater-level drawdown near the well. 

Groundwater levels are also lowered when recharge is reduced, for example during late summer 

months when evapotranspiration is high. Contaminants in streams may be transported to the 

groundwater under these conditions. 

Although the regional model was not specifically designed to examine 

groundwater/surface-water interactions, simulated stream infiltration to groundwater is prevalent 

for the 1999 flow conditions, the period with the largest well withdrawals simulated. Stream loss 

was simulated at several locations on Park Creek and Little Neshaminy Creek downstream from 

Willow Grove NASJRB (fig. 36). These locations are mostly near pumping wells and 

downgradient of the diabase dike. The dike may contribute to stream loss by acting as a 

subsurface dam, raising groundwater levels on the upstream side of the dike, but lowering levels 

on the downstream side. Contaminants that discharge to Park Creek through base flow, or from 

surface runoff, or from point discharges, could infiltrate back into the aquifer. 

Figure 36. Detail near Willow Grove, Park Creek, and Little Neshaminy Creek of ModelMuse screenshot 

of simulated potentiometric surface in model layer 1, regional flow paths from possible per- and 

polyfluoroalkyl substances (PFAS) contamination in groundwater, and locations of stream loss for 

conditions in 1999. Contour interval is 20 ft (see fig. 26).
67

Well withdrawals in the Willow Grove area were much less in the 2017 simulation than 

in the 1999 simulation (fig. 18). As a result, less streamflow loss to groundwater was simulated 

in 2017 (fig. 37) and contaminants that enter Little Neshaminy Creek near the Willow Grove 

base would flow downstream, rather than infiltrate to groundwater from Little Neshaminy Creek. 

Figure 37. ModelMuse screenshot detail near Willow Grove, Park Creek, and Little Neshaminy Creek of 

the simulated potentiometric surface in model layer 1, regional flow paths from possible per- and 

polyfluoroalkyl substances (PFAS) contamination in groundwater, and locations of stream loss for 

conditions in 2017. The contour interval is 20 ft (see appendix 3, fig. 3.21).
68

Limitations and Suggestions for Improvements

The extent to which the simulations in this report correspond to field conditions cannot be 

fully verified or validated, owing to the lack of data and the simplifications made in the 

preliminary conceptual and numerical model developed in this report. Some of the major 

limitations, and suggestions for improvements, are described in this section. 

The model grid cells and layers are large relative to the scale of individual fractures and 

fracture zones through which the bulk of the groundwater is transported. Most of the aquifer 

space is composed of unfractured rock that does not contribute substantially to the aquifer’s 

hydraulic conductivity. The model discretizes the aquifer into large blocks (cells) that represent 

the bulk properties of the aquifers and describe the flow paths at the scale of the cells, or larger. 

Nonetheless, the model is considered to usefully approximate regional flow conditions and paths. 

The regional model can provide boundary conditions for smaller models with higher grid 

resolution. The smaller models would be useful in further exploring the relation between possible 

sources and contaminated private wells near the bases that do not appear on flow paths simulated 

with the regional model. 

The model simulates steady-state flow conditions when actual flow conditions are 

transient. The transient conditions result from both long-term changes in pumping, variations in 

recharge, and short-term changes in pumping schedules, such as withdrawal wells cycling on and 

off each day. The flow conditions, and the resulting flow paths, change with time, but the steadystate model only represents unchanging average conditions. Transient flow models could 

simulate the impact of changes in recharge, pumpage, and well cycling on flow paths. 

Available data indicate that streams lose water in selected stream reaches during lowflow conditions in late summer and where infiltration is induced by pumping wells. 

Contaminants present in surface water may enter the groundwater system through losing reaches 

but these possible contaminant pathways were not simulated by the model. Additional data on 

streamflow at more locations in the study area under various flow conditions would improve 

characterization of groundwater/surface-water interaction.

The stream stages specified in the model were interpolated from a DEM and do not 

accurately represent actual stages throughout the model domain. The representation of the 

exchange of water between groundwater and streams would be improved by more accurately 

assigning stream stages, through a better interpolation procedure or a finer model grid near 

streams. 

All water levels in the model calibration are assumed to have the same weight, and thus 

to be of equal importance to the calibration. Some of the water-level data are correlated with 

other nearby measurements, especially at contaminated sites that have been the focus of detailed 

field investigations, and it may be appropriate to account for these correlations in the 

measurement weighting scheme (Hill and Tiedeman, 2007; Doherty, 2015). 

Measurement weights could also be adjusted to account for data errors, for example those 

caused by estimating land surface at the well from a 1:24,000-scale topographic map for the 

older data records. Well-construction data are missing for some wells, and the measured water 

level may not be accurately assigned to the correct model layer without these data. Further, some 

wells may be open to a long interval, but represent the water level in a discrete fracture or 

fracture zone within that open interval. Wells for which such detailed characterization data are 

available could be carefully analyzed and assigned appropriate weights to increase their 

importance for model calibration. It may be useful to measure and compile detailed well 
69

characteristics data in important model areas. Other limitations to model results are associated 

with uncertainty and (or) incompleteness in water-use data and characterization. 

Further analysis of model residuals, especially for well and stream locations with high 

quality data, could suggest refinements to the model structure, which is presently limited to 

geologic zones based on the regional geology. The delineation of the upper Stockton from the 

middle and lower members of the Stockton Formation is hypothetical for the unmapped Bucks 

County portion of the model. Further geologic investigations and mapping could guide additional 

refinements of the conceptual and numerical models, especially near areas of interest. 

The transport of contaminants in groundwater in fractured-rock aquifers is affected by a 

number or processes that are not included in the evaluation of flow paths in this report. These 

processes include dispersion, diffusion, sorption, and reactions. Also, additional information is 

needed about the distribution and magnitude of PFAS sources areas and transport mechanisms. 

Summary

In 2014, groundwater samples from residential and public supply wells in the vicinity of 

two former U.S. Navy bases at Willow Grove and Warminster and an active Air National Guard 

Station at Horsham, in Montgomery and Bucks Counties, southeastern Pennsylvania, were found 

to have concentrations of perfluorooctanoic acid (PFOA) and perfluorooctane sulfonate (PFOS), 

which are per- and polyfluoroalkyl substances (PFAS), above U.S. Environmental Protection 

Agency (EPA) provisional health advisory (HA) levels. In 2016, the EPA established a new 

lower, lifetime HA of 70 parts per trillion for combined concentrations of PFOS and PFOA, 

which resulted in the shutdown of additional supply wells. Groundwater has been, and remains, a 

major source of drinking water supply in the area. At the request of the U.S. Navy, and in 

consultation with other Federal and State agencies and local stakeholders, the U.S. Geological 

Survey developed a numerical model to simulate regional groundwater flow and flow paths from 

potential PFAS sources at the bases under selected conditions during 1999–2017 at and near the 

former Willow Grove NASJRB and NAWC Warminster to provide information to help address 

the PFAS contamination in groundwater. 

The bases and surrounding areas are underlain by Triassic geologic strata that form 

fractured-rock aquifers, which are recharged locally by precipitation. The sedimentary Triassic 

rocks (sandstones and shales), principally the Stockton Formation and overlying Lockatong 

Formation in the study area, generally strike east-northeast and dip to the northwest. 

Groundwater discharges locally to streams and wells. To support the model development and 

calibration, the U.S. Geological Survey compiled hydrologic data including water levels in wells, 

streamflow, point discharges to streams, groundwater withdrawals (pumping wells and quarries), 

and precipitation for selected periods between 1999–2017. Stream base flow, which can be used 

to estimate recharge and was used for model calibration, was estimated by hydrograph separation 

methods at seven stream gages in the study area. Data from years 1999, 2010, 2013, and 2016 

were used for model calibration and data for year 2017 were used for simulation only. 

Groundwater withdrawals near the bases from public supply and other large wells decreased 

substantially from the 1990s to 2017, increasing the proportion of groundwater recharge that 

discharged to local streams.
70

A three-dimensional numerical model, MODFLOW 6, was used to simulate steady-state 

flow for the selected time periods with different recharge rates and well withdrawals. The model 

structure included seven layers of varying hydraulic properties and three-dimensional hydraulic 

conductivity tensor to estimate the anisotropic character of the dipping-bed fractured-rock 

aquifer in the area. The model was calibrated using 1,009 groundwater levels and 17 stream 

base-flow estimates using the parameter estimation software PEST++. The model simulations 

provide information on regional groundwater-flow paths and groundwater/surface-water 

exchanges, in relation to changes in water supply pumping and recharge rates. 

Regional groundwater-flow paths from possible PFAS sources at the former bases were 

simulated for 1999, 2010, 2013, 2016, and 2017 conditions. Simulations showed that recharge at 

the bases discharged to withdrawal wells and local streams, generally within a mile or two of the 

bases. Supply and remediation wells at the bases captured some of the recharge on base areas of 

possible PFAS contamination, whereas other base recharge was simulated to flow to nearby 

public supply wells and streams, depending on water use and aquifer recharge conditions 

between 1999 and 2017. Locations of many residential wells near the bases identified by the U.S. 

Navy and Air National Guard as having elevated PFAS concentrations are generally consistent 

with the simulated flow paths from possible sources at the bases. However, there are some areas 

of observed PFAS contamination where no flow paths from base sources were simulated, 

indicating the presence of unknown PFAS sources, unidentified transport processes, and (or) 

model limitations. Additionally, no data are available on PFAS concentrations in groundwater in 

some areas of simulated flow paths from base sources. 

The preliminary model used here has limitations that should be considered when using 

these results. Further investigations to improve this evaluation could include additional data 

collection, such as synoptic water levels with well construction information, and refined 

simulation with transient flow conditions. 
71

References Cited

Bakker, M., Post, V., Langevin, C.D., Hughes, J.D., White, J.T., Starn, J.J. and Fienen, M.N., 

2016, Scripting MODFLOW model development using Python and FloPy: Groundwater, v. 54, 

p. 733–739, accessed October 19, 2018, at https://doi.org/10.1111/gwat.12413.

Bakker, M., Post, V., Langevin, C.D., Hughes, J.D., White, J.T., Starn, J.J., and Fienen, M.N., 

2018, FloPy—Python package for creating, running, and post-processing MODFLOW-based 

models (ver. 3.2.10): U.S. Geological Survey software release, accessed October 19, 2018, at

https://doi.org/10.5066/F7BK19FH. 

Barlow, P.M., Cunningham, W.L., Zhai, T., and Gray, M., 2017, The groundwater toolbox—A 

graphical and mapping interface for analysis of hydrologic data (ver. 1.3.1): U.S. Geological 

Survey software release, accessed March 13, 2018, at https://doi.org/10.5066/F7R78C9G. 

Battelle Memorial Institute, 2016, Final fourth five-year review report for former NAWC 

Warminster, PA: Prepared under Contract No. N62583-11-D-0515 Task Order No. 0063 for 

Naval Facilities Engineering Command BRAC Program Management Office East, November 

2016, accessed April 19, 2019, at https://www.navfac.navy.mil/content/navfac/en/

products_and_services/ev/products_and_services/env_restoration/administrative_records.html

?fromDate=DD-MON-YYYY&toDate=DD-MONYYYY&p_instln_id=WARMINSTER_NAWC&basic=&title=&sites=&author=&keywords= 

Berg, T.M., Edmunds, W.E., Geyer, A.R., and others, comps., 1980, Geologic map of 

Pennsylvania: Pennsylvania Geological Survey, 4th Series, map 1, 2d ed., scale 1:250,000 

Bird, P.H., 1998, Geohydrology and ground-water quality of Warwick Township, Bucks County, 

Pennsylvania: U.S. Geological Survey Water-Resources Investigations Report 97–4267, 37 p., 

accessed April 19, 2019, at https://pubs.er.usgs.gov/publication/wri974267. 

CDM Federal Programs Corporation, 2003, Annual groundwater monitoring report for May 

2003 sampling event, Raymark Superfund Site, Hatboro, Pennsylvania: Report prepared for 

EPA under contract no. 68-S7-3003, work assignment no. 020-RALR-0378, document control 

no. 3232-020-RT-OTHR-01880, July 9, 2003, accessed November 21, 2019, at 

https://semspub.epa.gov/work/03/2073046.pdf. 

Conger, R.W., 1997, Evaluation of geophysical logs, phase I, at Willow Grove Naval Air 

Station, Montgomery County, Pennsylvania: U.S. Geological Survey Open-File Report 97–

631, 49 p., accessed July 9, 2017, at https://doi.org/10.3133/ofr97631. 

Conger, R.W., 1998, Identification of water-bearing zones by the use of geophysical logs and 

borehole television surveys, collected February to September 1997, at the Former Naval Air 

Warfare Center, Warminster, Bucks County, Pennsylvania: U.S. Geological Survey Open-File 

Report 98–86, 26 p., accessed July 9, 2017, at https://pubs.er.usgs.gov/publication/ofr9886. 

Conger, R.W., 1999, Evaluation of geophysical logs, phase II, at Willow Grove Naval Air 

Station Joint Reserve Base, Montgomery County, Pennsylvania: U.S. Geological Survey 

Open-File Report 99–73, 14 p., accessed July 9, 2017, at https://doi.org/10.3133/ofr9973. 

Conger, R.W., and Bird, P.H., 1999, Identification of water-bearing fractures by the use of 

geophysical logs, May to July 1998, former Naval Air Warfare Center, Bucks County, 

Pennsylvania: U.S. Geological Survey Open-File Report 99–215, 27 p., accessed July 9, 2017, 

at https://pubs.er.usgs.gov/publication/ofr99215. 

Delaware River Basin Commission, 1980, Resolution no. 80–18: Delaware River Basin 

Commission web page, accessed November 10, 2019, at 

https://www.state.nj.us/drbc/library/documents/Res1980-18gwpa.pdf. 
72

Delaware River Basin Commission, 2019, Southeastern Pennsylvania ground water protected 

area (GWPA) GIS data sets: Delaware River Basin Commission web page, accessed 

November 1, 2019, at https://www.nj.gov/drbc/programs/project/pr/gwpa-data.html. 

Delaware Valley Regional Planning Commission, 2013, Land use in the Delaware Valley, 2010, 

analytical data report 020: Delaware Valley Regional Planning Commission Publication 

ADR020, 88 p., accessed November 27, 2019, at https://www.dvrpc.org/Products/ADR020/. 

Delaware Valley Regional Planning Commission, 2019, Land use in the Delaware Valley, 2010, 

digital data: Delaware Valley Regional Planning Commission excel file, accessed November 

27, 2019, at https://www.dvrpc.org/Products/ADR020/. 

Doherty, J., 2015, Calibration and uncertainty analysis for complex environmental models: 

Brisbane, Australia, Watermark Numerical Computing, 227 p. 

Doherty, J.E., and Hunt, R.J., 2010, Approaches to highly parameterized inversion—A guide to 

using PEST for groundwater-model calibration: U.S. Geological Survey Scientific 

Investigations Report 2010–5169, 59 p. 

Esri, 2019, ArcGIS desktop release 10.7: Esri, Redlands, Calif., accessed January 3, 2019, at 

http://desktop.arcgis.com/en/. 

Fan, Y., Toran, L., and Schlische, R.W., 2007, Groundwater flow and groundwater‐stream

interaction in fractured and dipping sedimentary rocks—Insights from numerical models: 

Water Resources Research, v. 43, no. 1, 13 p., accessed July 9, 2017, at 

https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2006WR004864. 

Gilmore and Associates, Inc., 2015, Annual hydrogeologic report no. 11, DRBC docket no. D97-16 CP-3, Horsham Water and Sewer Authority, Horsham, Pennsylvania: Report prepared

for Horsham Water and Sewer Authority, File no. 15–01135, May 17, 2016, 117 p.

Goode, D.J., 2016, Map visualization of groundwater withdrawals at the sub-basin scale: 

Hydrogeology Journal, v. 24, p. 1057–1065, accessed July 9, 2017, at

https://doi.org/10.1007/s10040-016-1379-x. 

Goode, D.J., and Senior, L.A., 2000, Simulation of aquifer tests and ground-water flowpaths at 

the local scale in fractured shales and sandstones of the Brunswick Group and Lockatong 

Formation, Lansdale, Montgomery County, Pennsylvania: U.S. Geological Survey Open-File 

Report 2000–97, 46 p., accessed July 9, 2017, at https://doi.org/10.3133/ofr0097. 

Goode, D.J., and Senior, L.A., 2020, MODFLOW 6 and MODPATH 7 model data sets used to 

evaluate groundwater flow in the vicinity of Horsham and Warminster, Bucks and 

Montgomery Counties, Pennsylvania—Preliminary simulations for conditions in 1999, 2010, 

2013, 2016, and 2017: U.S. Geological Survey data release, 

https://doi.org/10.5066/P9K36P5S. 

Greenman, D.W., 1955, Groundwater resources of Bucks County, Pennsylvania: Pennsylvania 

Geological Survey, 4th Series, Bulletin W11, 66 p., 1 plate. 

H&S Environmental, Inc., 2012, Semi-annual monitoring report for the fall (November 2011) 

performance monitoring event for long-term environmental monitoring at Operable Units 1A, 

3, and 4 at the former NAWC, Warminster, Pennsylvania: Report prepared for Naval Facilities 

Engineering Command Mid-Atlantic BRAC Program Management Office Northeast under 

contract no. N40085-10-D-9409, task order no. 004, June 2012, 234 p. 

Hill, M.C., and Tiedeman, C.R., 2007, Effective groundwater model calibration, with analysis of 

sensitivities, predictions, and uncertainty: New York, Wiley and Sons, 455 p. 

Koman Government Solutions LLC, 2018a, Annual monitoring report for the spring (May 2017) 

performance monitoring event for long-term environmental monitoring at Operable Units 1A, 
73 

3, and 4 at the former NAWC, Warminster, Pennsylvania: Report prepared for Naval Facilities 

Engineering Command Mid-Atlantic BRAC Program Management Office Northeast under 

Contract No. N40085-10-D-9409, contract task order no. 004, April 2019, accessed 

November 9, 2019, at 

https://www.navfac.navy.mil/niris/MID_ATLANTIC/WARMINSTER_NAWC/N62269_0012

27.pdf.

Koman Government Solutions LLC, 2018b, Former NAWC, Warminster, PA December 2017 

groundwater treatment plant report: Report prepared for Naval Facilities Engineering 

Command Mid-Atlantic BRAC Program Management Office Northeast under Contract No. 

N40085-10-D-9409, contract task order no. 004, January 2018, 136 p. 

Koman Government Solutions LLC, 2018c, Semi-annual monitoring report for the Fall 

(November 2016) performance monitoring event for long-term environmental monitoring at 

Operable Units 1A, 3, and 4 at the former NAWC, Warminster, Pennsylvania: Report prepared 

for Naval Facilities Engineering Command Mid-Atlantic BRAC Program Management Office 

Northeast under Contract No. N40085-10-D-9409, contract task order no. 004, April 2018, 

278 p. 

Langevin, C.D., Hughes, J.D., Banta, E.R., Niswonger, R.G., Panday, S., and Provost, A.M., 

2017, Documentation for the MODFLOW 6 Groundwater Flow Model: U.S. Geological 

Survey Techniques and Methods, book 6, chap. A55, 197 p., accessed July 9, 2018, at 

https://doi.org/10.3133/tm6A55. 

Langevin, C.D., Hughes, J.D., Banta, E.R., Provost, A.M., Niswonger, R.G., and Panday, S., 

2018, MODFLOW 6 Modular Hydrologic Model (ver. 6.0.3): U.S. Geological Survey 

software release, accessed August 9, 2018, at https://doi.org/10.5066/F76Q1VQV. 

Leidos, 2018, Groundwater monitoring report for the perfluorinated compound facility 

investigation at the Horsham Air Guard Station (111th Attack Wing), Horsham, Pennsylvania: 

Report prepared for NGB/A4OR Joint Base Andrews, contract W9133L-14-D-0007, delivery 

order number 0005, August 2018, accessed November 10, 2019, at http://afcec.publicadminrecord.us.af.mil/Search.aspx. 

Longwill, S.M., and Wood, C.R., 1965, Ground-water resources of the Brunswick Formation in 

Montgomery and Berks Counties, Pennsylvania: Pennsylvania Geological Survey, 4th Series, 

Bulletin W22, 59 p., 1 plate.

Low, D.J, Hippe, D.J., and Yannacci, D., 2002, Geohydrology of southeastern Pennsylvania: 

U.S. Geological Survey Water-Resources Investigation Report 2002–4166, 347 p. [Also 

available at https://doi.org/10.3133/wri004166.]

Lyttle, P.T., and Epstein, J.B., 1987, Geologic map of the Newark 1 degree x 2 degrees 

quadrangle, New Jersey, Pennsylvania and New York: U.S. Geological Survey Miscellaneous 

Investigations Map I-1715, scale 1:250,000. [Also available at https://doi.org/10.3133/i1715; 

Georeferenced version of Map I-1715 available at 

https://ngmdb.usgs.gov/Prodesc/proddesc_9892.htm.] 

Montgomery County Planning Commission, 1994, 1993 Status report sewage treatment facilities 

in Montgomery County: Norristown, Pennsylvania.

National Oceanic and Atmospheric Administration, 2018a, Climate at a glance, county time 

series, Bucks County, Pennsylvania: National Centers for Environmental Information web 

page, accessed December 22, 2018, at https://www.ncdc.noaa.gov/cag/.
74 

National Oceanic and Atmospheric Administration, 2018b, Climate at a glance, county time 

series, Montgomery County, Pennsylvania: National Centers for Environmental Information 

web page, accessed December 22, 2018, at https://www.ncdc.noaa.gov/cag/.

Naval Facilities Engineering Command, 2018, Final second five-year review for former Naval 

Air Station Joint Reserve Base Willow Grove Horsham, Pennsylvania: Naval Facilities 

Engineering Command BRAC Program Management Office East, September 2018, accessed 

April 19, 2019, at https://www.navfac.navy.mil/niris/MID_

ATLANTIC/WILLOW_GROVE_NAS/N00158_000760.pdf.

Newport, T.G., 1971, Ground-water resources of Montgomery County, Pennsylvania: 

Pennsylvania Geological Survey, 4th Series, Bulletin W29, 83 p., 2 plates. 

Penn State Institutes of Energy and Environment, 2015, PAMAP Program 3m county mosaic 

digital elevation model: Pennsylvania Spatial Data Access web page, accessed May 2016 from 

ftp://ftp.pasda.psu.edu/pub/pasda/pamap/LidarMosaics/CountyMosaics/county_DEM_3M/. 

Pennsylvania Department of Environmental Protection, 2019a, Wastewater reports, eDMR data: 

Pennsylvania Department of Environmental Protection Bureau of Clean Water web page, 

accessed March 2019 at 

https://www.dep.pa.gov/Business/Water/CleanWater/Pages/Wastewater-Reports.aspx.

Pennsylvania Department of Environmental Protection, 2019b, Water reports, water use reports, 

water quantity report by source: Pennsylvania Department of Environmental Protection web 

page, accessed November 2019 at 

https://www.dep.pa.gov/DataandTools/Reports/Pages/Water.aspx. 

Pennsylvania Department of Environmental Protection, 2019c, Water reports, water use reports, 

water use report by facility: Pennsylvania Department of Environmental Protection web page, 

accessed November 2019 at https://www.dep.pa.gov/DataandTools/Reports/Pages/Water.aspx. 

Pennsylvania Department of Conservation & Natural Resources Bureau of Topographic and 

Geological Survey, 2018, Pennsylvania Groundwater Information System: Pennsylvania 

Department of Conservation & Natural Resources web page, accessed November 2018 at 

https://www.dcnr.pa.gov/Conservation/Water/Groundwater/PAGroundwaterInformationSyste

m/Pages/default.aspx. 

Pollock, D.W., 2016, User guide for MODPATH version 7—A particle-tracking model for 

MODFLOW: U.S. Geological Survey Open-File Report 2016–1086, 35 p., accessed July 9, 

2018, at https://doi.org/10.3133/ofr20161086.

Pollock, D.W., 2017, MODPATH v7.2.01—A particle-tracking model for MODFLOW: U.S. 

Geological Survey software release, accessed November 2017 at

https://doi.org/10.5066/F70P0X5X. 

Provost, A.M., Langevin, C.D., and Hughes, J.D., 2017, Documentation for the "XT3D" option 

in the Node Property Flow (NPF) Package of MODFLOW 6: U.S. Geological Survey 

Techniques and Methods, book 6, chap. A56, 40 p., accessed July 9, 2018, at

https://doi.org/10.3133/tm6A56. 

Resolution Consultants, 2019, Remedial investigation report, per- and polyfluoralkyl substances 

investigation activities NASJRB Willow Grove, PA: Report prepared for Naval Facilities 

Engineering Command (NAVFAC) Mid-Atlantic, contract N62470-11-D-8013, Final version 

1, September 20, 2019, accessed November 10, 2019, at 

https://www.navfac.navy.mil/products_and_services/ev/products_and_services/env_restoration

/administrative_records.html?fromDate=DD-MON-YYYY&toDate=DD-MONYYYY&p_instln_id=WILLOW_GROVE_NAS&basic=&title=&sites=&author=&keywords=. 
75

Rima, D.R., Meisler, H., and Longwill, S., 1962, Geology and hydrology of the Stockton 

Formation in southeastern Pennsylvania: Pennsylvania Geological Survey, 4th Series, Bulletin 

W14, 110 p., 4 plates. 

Rowland, C.J., 1997, Altitude and configuration of the potentiometric surface in Warwick 

Township, Bucks County, Pennsylvania, September 1994 through May 1995: U.S. Geological 

Survey Open-File Report 97–554, 1 plate, accessed July 9, 2017, at 

https://doi.org/10.3133/ofr97554. 

Schlische, R., 1992, Structural and stratigraphic development of the Newark extensional basin, 

eastern North America—Evidence for the growth of the basin and its bounding structures: 

Geological Society of America Bulletin, v. 104, p. 1246–1263. 

Schreffler, C.L., 1996, Water-use analysis program for the Neshaminy Creek basin, Bucks and 

Montgomery counties, Pennsylvania: U.S. Geological Survey Water-Resources Investigations 

Report 96–4127, 85 p. 

Senior, L.A., and Goode, D.J., 1999, Ground-water system, estimation of aquifer hydraulic 

properties, and effects of pumping on ground-water flow in Triassic sedimentary rocks in and 

near Lansdale, Pennsylvania: U.S. Geological Survey Water-Resources Investigations Report 

99–4228, 112 p. 

Senior, L.A., and Goode, D.J., 2013, Investigations of groundwater system and simulation of 

regional groundwater flow for North Penn Area 7 Superfund site, Montgomery County, 

Pennsylvania (ver. 1.1, April 2015): U.S. Geological Survey Scientific Investigations Report 

2013–5045, 95 p., accessed July 9, 2017, at https://pubs.usgs.gov/sir/2013/5045/. 

Senior, L.A., and Goode, D.J., 2017, Effects of changes in pumping on regional groundwaterflow paths, 2005 and 2010, and areas contributing recharge to discharging wells, 1990–2010, 

in the vicinity of North Penn Area 7 Superfund site, Montgomery County, Pennsylvania: U.S. 

Geological Survey Scientific Investigations Report 2017–5014, 36 p., accessed July 9, 2017, at 

https://doi.org/10.3133/sir20175014. 

Senior, L.A., Ludlow, R.A., and Zarr, L.F, 2020, Groundwater levels, groundwater withdrawals, 

and point-source discharges to streams in the vicinity of Willow Grove and Warminster, Bucks 

and Montgomery Counties, Pennsylvania, for selected years during 1999–2017: U.S. 

Geological Survey data release, https://doi.org/10.5066/P9ZGEI67. 

Senior, L.A., and Ruddy, A.J., 2004, Altitude and configuration of the water-level surface in 

Mesozoic sedimentary rocks at and near the North Penn Area 7 Superfund Site, Upper 

Gwynedd Township, Montgomery County, Pennsylvania, December 4–6, 2000: U.S. 

Geological Survey Open-File Report 2004–1006, 1 sheet, 1:24,000, accessed July 9, 2017, at 

https://pubs.er.usgs.gov/publication/ofr041006.

Sevon, W.D., 2000, comp., Physiographic provinces of Pennsylvania: Pennsylvania Geological 

Survey Map 13, scale 1:2,000,000, 4th ed., accessed January 31, 2019, at 

http://www.docs.dcnr.pa.gov/cs/groups/public/documents/document/dcnr_016202.pdf. 

Sloto, R.A., 1997, Results of borehole geophysical logging and aquifer-isolation tests conducted 

in the John Wagner and Sons, Inc former production well, Ivyland, Pennsylvania: U.S. 

Geological Survey Water-Resources Investigations Report 1997–4095, 18 p., accessed July 9, 

2017, at https://pubs.er.usgs.gov/publication/wri974095. 

Sloto, R.A., 1998, Results of borehole geophysical logging and hydraulic tests conducted in Area 

D supply wells, former U.S. Naval Air Warfare Center, Warminster, Pennsylvania: U.S. 

Geological Survey Water-Resources Investigations Report 98–4129, 47 p. 
76

Sloto, R.A., 2002, Hydrogeological investigation at Site 5, Willow Grove Naval Air Station/Joint 

Reserve Base, Horsham Township, Montgomery County, Pennsylvania: U.S. Geological 

Survey Water-Resources Investigations Report 2001–4263, 37 p., accessed July 9, 2017, at 

https://pubs.er.usgs.gov/publication/wri014263. 

Sloto, R.A., 2007, Interpretation of borehole geophysical logs, aquifer-isolation tests, and waterquality data for sites 1, 3, and 5 at the Willow Grove Naval Air Station/Joint Reserve Base, 

Horsham Township, Montgomery County, Pennsylvania—2005: U.S. Geological Survey 

Scientific Investigations Report 2006–5243, 74 p., accessed July 9, 2017, at 

https://doi.org/10.3133/sir20065243. 

Sloto, R.A., 2008, Interpretation of borehole geophysical logs at Area C, former Naval Air 

Warfare Center, Warminster Township, Bucks County, Pennsylvania, 2007: U.S. Geological 

Survey Open-File Report 2008–1207, 20 p., accessed July 9, 2017, at 

https://doi.org/10.3133/ofr20081207. 

Sloto, R.A., 2010, Changes in groundwater flow and volatile organic compound concentrations 

at the Fischer and Porter Superfund Site, Warminster Township, Bucks County, Pennsylvania, 

1993–2009: U.S. Geological Survey Scientific Investigations Report 2010–5054, 115 p. 

Sloto, R.A., and Buxton, D.E., 2007, Estimated ground-water availability in the Delaware River 

Basin, 1997–2000 (ver. 1.1, May 2007): U.S. Geological Survey Scientific Investigations 

Report 2006–5125, 67 p. 

Sloto, R.A., Conger, R.W., and Grazul, K.E., 1998, Geohydrology and distribution of volatile 

organic compounds in ground water in the Casey Village area, Bucks County, Pennsylvania: 

U.S. Geological Survey Water-Resources Investigations Report 1998–4010, 81 p., accessed 

July 9, 2017, at https://pubs.er.usgs.gov/publication/wri984010. 

Sloto, R.A., and Crouse, M.Y., 1986, HYSEP—A computer program for streamflow hydrograph 

separation and analysis: U.S. Geological Survey Water-Resources Investigations Report 86–

4040, 46 p., accessed March 28, 2019, at https://pubs.er.usgs.gov/publication/wri964047. 

Sloto, R.A., and Davis, D.K., 1982, Effect of urbanization on the water resources of Warminster 

Township, Bucks County, Pennsylvania: U.S. Geological Survey Water-Resources 

Investigations Report 82–4010, 78 p., accessed March 28, 2019, at 

https://pubs.er.usgs.gov/publication/wri824020.

Sloto, R.A., Goode, D.J., and Frasch, S.M., 2002, Interpretation of borehole geophysical logs, 

aquifer-isolation tests, and water quality, supply wells 1 and 2, Willow Grove Naval Air 

Station/Joint Reserve Base, Horsham Township, Montgomery County, Pennsylvania: U.S. 

Geological Survey Water-Resources Investigations Report 2001–4264, 64 p., accessed July 9, 

2017, at https://pubs.er.usgs.gov/publication/wri014264. 

Sloto, R.A., Goode, D.J., and Way, J.H., 2001, Altitude and configuration of the regional 

potentiometric surface, Willow Grove Naval Air Station/Joint Reserve Base and vicinity, 

Horsham Township, Montgomery County, Pennsylvania, October 7–8, 1999: U.S. Geological 

Survey Open-File Report 2001–0149, 13 p., accessed July 9, 2017, at 

https://pubs.er.usgs.gov/publication/ofr01149. 

Sloto, R.A., and Grazul, K.E., 1995, Altitude and configuration of the potentiometric surface, 

Casey Village, Warminster and Upper Southampton Townships, Bucks County, Pennsylvania, 

August 3, 1995: U.S. Geological Survey Open-File Report 95–717, 2 maps on one sheet, 

accessed July 9, 2017, at https://doi.org/10.3133/ofr95717. 

Sloto, R.A., Macchioroli, P., and Towle, M.T., 1996, Geohydrology of the Stockton Formation 

and cross-contamination through open boreholes, Hatboro Borough and Warminster 
77 

Township, Pennsylvania: U.S. Geological Survey Water-Resources Investigations Report 86–

4047, 49 p., 2 plates, accessed March 28, 2019, at 

https://pubs.er.usgs.gov/publication/wri964047.

Sloto, R.A., and Schreffler, C.L., 1994, Hydrogeology and ground-water quality of northern 

Bucks County, Pennsylvania: U.S. Geological Survey Water-Resources Investigations Report 

94–4109, 85 p., 1 plate. [Also available at https://doi.org/10.3133/wri944109.]

Solley, W.B., Pierce, R.R., and Perlman, H.A., 1998, Estimated use of water in the United States 

in 1995: U.S. Geological Survey Circular 1200, 71 p., accessed November 1, 2019, at 

https://doi.org/10.3133/cir1200. 

Tetra Tech, 2014, Evaluation of potential source of perfluorinated compounds, former Naval Air 

Warfare Center Warminster, Pennsylvania: Report prepared for Naval Facilities Engineering 

Command Mid-Atlantic, Norfolk, Virginia, contract number N62470-08-D-1001, contract task 

orderWE23, report N62269.AR.001000 NAWC WARMINSTER 5090.3a, November 2014, 

accessed November 11, 2019, at https://www.navfac.navy.mil/niris/MID_

ATLANTIC/WARMINSTER_NAWC/N62269_001000.pdf. 

Turner-Peterson, C.E., and Smoot, J.P., 1985, New thoughts on facies relationships in the 

Triassic Stockton and Lockatong Formations, Pennsylvania and New Jersey, in Gilpin, G.R., 

Jr., and Froehlich, A.J., eds., Proceedings of the second U.S. Geological Survey workshop on 

the early Mesozoic basins of the Eastern United States: U.S. Geological Survey Circular 946, 

p. 10–17, accessed July 9, 2017, at https://doi.org/10.3133/cir946.

U.S. Census Bureau, 2019, Quick facts: U.S. Census Bureau web page, accessed July 9, 2019, at

https://www.census.gov/quickfacts/fact/table/northamptontownshipbuckscountypennsylvania,h

atboroboroughpennsylvania,warringtontownshipbuckscountypennsylvania,warwicktownshipb

uckscountypennsylvania,warminstertownshipbuckscountypennsylvania,horshamtownshipmont

gomerycountypennsylvania/PST045218. 

U.S. Environmental Protection Agency, 2006, Pilot region-based optimization program for fundlead sites in EPA Region 3, Site optimization tracker: Raymark, Hatboro, Pennsylvania: EPA 

report 542-R-06-006f, December 2006, 53 p. 

U.S. Environmental Protection Agency, 2012, Third unregulated contaminant monitoring rule: 

U.S. Environmental Protection Agency web page, accessed January 31, 2019, at 

https://www.epa.gov/dwucmr/third-unregulated-contaminant-monitoring-rule. 

U.S. Environmental Protection Agency, 2014, Peer review of health effects documents for PFOA 

and PFOS: U.S. Environmental Protection Agency, Office of Water, EPA-820-F-14-004, 2 p., 

accessed April 10, 2015, at http://water.epa.gov/drink/standards/upload/Peer-Review-ofHealth-Effects-Documents-for-PFOA-and-PFOS-Factsheet-February-2014.pdf. 

U.S. Environmental Protection Agency, 2016, Fact sheet PFOA & PFOS drinking water health 

advisories: U.S. Environmental Protection Agency web page, EPA 800-F-16-003, 5 p., 

accessed March 28, 2018, at https://www.epa.gov/sites/production/files/2016-

06/documents/drinkingwaterhealthadvisories_pfoa_pfos_updated_5.31.16.pdf. 

U.S. Geological Survey, 1998, Southeastern Pennsylvania ground water protected area: U.S. 

Geological Survey Open-File Report 98–571, 1 CD, accessed August 12, 2019, at 

https://pubs.er.usgs.gov/publication/ofr98571. 

U.S. Geological Survey, 2005, USGS National Hydrography Dataset (NHD)–High-resolution 

flowline features: U.S. Geological Survey National Hydrography Dataset, accessed June 3, 

2009, at https://www.pasda.psu.edu/uci/DataSummary.aspx?dataset=364. 
78

U.S. Geological Survey, 2016, USGS National Hydrography Dataset (NHD) downloadable data 

collection: U.S. Geological Survey National Hydrography Dataset, accessed November 7, 

2019, at https://www.usgs.gov/core-science-systems/ngp/tnm-delivery/. 

U.S. Geological Survey, 2018, USGS water data for the Nation: U.S. Geological Survey 

National Water Information System database, accessed September 16, 2018, 

at https://doi.org/10.5066/F7P55KJN. 

U.S. Navy, 2019a, Naval Facilities Engineering Command Former Naval Air Station Joint 

Reserve Base Willow Grove documents: Naval Facilities Engineering Command web page, 

accessed November 5, 2019, at 

https://www.bracpmo.navy.mil/brac_bases/northeast/reserve_base_willow_grove/documents.h

tml.

U.S. Navy, 2019b, Naval Facilities Engineering Command Former Naval Air Warfare Center 

Warminster documents: Naval Facilities Engineering Command web page, accessed 

November 5, 2019, at 

https://www.bracpmo.navy.mil/brac_bases/northeast/former_warfare_center_warminster/docu

ments.html.

Welter, D.E., White, J.T., Hunt, R.J., and Doherty, J.E., 2015, Approaches in highly 

parameterized inversion—PEST++ Version 3, a Parameter ESTimation and uncertainty 

analysis software suite optimized for large environmental models: U.S. Geological Survey 

Techniques and Methods, book 7, section C12, 54 p. 

White, J., Welter, D., and Doherty, J., 2019, PEST++ Version 4: 170 p., accessed November 26, 

2019, at http://www.pesthomepage.org/PEST_and_PESTPP.php. 

White, K.E., and Sloto, R.A., 1990, Base-flow frequency characteristics of selected Pennsylvania 

streams: U.S. Geological Survey Water-Resources Investigations Report 90–4160, 67 p. 

Winston, R.B., 2019, ModelMuse version 4—A graphical user interface for MODFLOW 6: U.S. 

Geological Survey Scientific Investigations Report 2019–5036, 10 p. 

Winston, R.B., and Goode, D.J., 2017, Visualization of groundwater withdrawals: U.S. 

Geological Survey Open-File Report 2017–1137, 8 p. 
79 

Appendix 1. Annual Base Flow as Determined from 

Measured Streamflow at SelectedGages and 

Estimated for Missing Streamflow Records During 

2010–2015 

Groundwater base-flow data were compiled from records at selected stream gages near 

Willow Grove and Warminster, Pennsylvania (U.S. Geological Survey, 2018). Base flow was 

estimated using daily mean streamflow values and HYSEP local-minimum hydrograph 

separation method in U.S. Geological Survey (USGS) Groundwater Toolbox (Barlow and 

others, 2017). The annual average base flow was computed from the average of monthly 

average base-flow values (tables 1.1–1.3). 

Missing values for stream flow at several stream gages were estimated by regression 

(table 1.4, fig. 1.1). USGS continuous gages with missing data were: 01467031 Pennypack 

Creek at Horsham; 01467036 Pennypack tributary at Hatboro; 01467039 Pennypack Creek at 

Willow Grove; and 014670413 Pennypack Creek at Bethayres. Regressions with 2010 monthly 

base flow (in cubic feet per second) at Pine Road gage was selected to fill missing monthly data 

for January-March 2010 at three Pennypack Creek gages: Horsham, Hatboro, Willow Grove. 

Regression with 2015 monthly base flow (in cubic feet per second) at Pine Road gage was 

selected to fill missing monthly data for October-December 2015 at four Pennypack Creek 

gages: Horsham, Hatboro, Willow Grove, Bethayres. Regressions with instantaneous values at 

Iron Works gage and daily values at Horsham gage (table 1.5, figs. 1.2 and 1.3) were used to 

estimate monthly average base flow at Iron Works gage from monthly average base flow at 

Horsham gage (table 1.6). 

Regression was used to estimate missing data for four gages and, for Iron Works Creek, 

to estimate average baseflow values. Regressions with 2010 monthly base flow (in cubic feet per 

second) at Pine Road gage was selected to fill missing monthly data for January–March 2010 at 

three Pennypack Creek gages: Horsham, Hatboro, Willow Grove. Regressions with 2015 

monthly baseflow (in cubic feet per second) at Pine Road gage was selected to fill missing 

monthly data for October–December 2015 at four Pennypack Creek gages: Horsham, Hatboro, 

Willow Grove, Bethayres. See tables 1.2 and 1.4 for more information about gages. 

References 

Barlow, P.M., Cunningham, W.L., Zhai, T., and Gray, M., 2017, The groundwater toolbox—A 

graphical and mapping interface for analysis of hydrologic data (ver. 1.3.1): U.S. Geological 

Survey software release, accessed March 13, 2018, at https://doi.org/10.5066/F7R78C9G. 

U.S. Geological Survey, 2018, USGS water data for the Nation: U.S. Geological Survey 

National Water Information System database, accessed September 16, 2018, at 

https://doi.org/10.5066/F7P55KJN. 
Year count

Calendar

year

Streamflow

CFS in CFS

Base flow

in CFS

Runoff

in BFP

1 1999 42.8 21.68 10.95 5.55 31.85 16.13 25.6

2 2000 39.2 19.92 15.75 8.00 23.48 11.93 40.1

3 2001 46.0 23.28 14.49 7.34 31.46 15.94 31.5

4 2002 28.1 14.22 10.81 5.47 17.26 8.74 38.5

5 2003 75.8 38.41 26.96 13.65 48.88 24.76 35.5

6 2004 60.1 30.55 22.99 11.68 37.16 18.87 38.2

7 2005 54.3 27.49 18.65 9.45 35.63 18.04 34.4

8 2006 66.4 33.62 23.04 11.67 43.34 21.95 34.7

9 2007 51.6 26.14 16.99 8.61 34.61 17.53 32.9

10 2008 46.6 23.68 18.06 9.17 28.56 14.51 38.7

11 2009 60.0 30.40 23.29 11.79 36.74 18.61 38.8

12 2010 46.8 23.71 15.94 8.07 30.87 15.64 34.1

13 2011 91.4 46.29 25.15 12.74 66.24 33.55 27.5

14 2012 30.7 15.59 13.23 6.72 17.47 8.87 43.1

15 2013 49.7 25.20 17.12 8.67 32.62 16.52 34.4

16 2014 52.6 26.63 18.81 9.53 33.77 17.11 35.8

17 2015 40.1 20.34 16.44 8.33 23.71 12.01 41.0

18 2016 33.7 17.14 15.14 7.69 18.61 9.45 44.9

19 2017 37.8 19.15 16.45 8.33 21.35 10.81 43.5

Minimum -- 28.1 14.22 10.81 5.47 17.26 8.74 25.6

Maximum -- 91.4 46.29 26.96 13.65 66.24 33.55 44.9

Median -- 46.8 23.71 16.99 8.61 31.85 16.13 35.8

Average -- 50.2 25.44 17.91 9.08 32.30 16.37 36.5

Table 1.1. Annual base flow (calendar year January−December) as determined by hydrograph separation 

methods for daily streamflow values at U.S. Geological Survey stream gage 01464907 Little Neshaminy Creek at 

Valley Road, 1999−2017, drainage area 26.8 square miles. HYSEP local-minimum hydrograph separation 

method in USGS Groundwater Toolbox (Barlow and others, 2017).

[CFS, annual average flow, in cubic feet per second; in, annual average flow per drainage area, in inches; BFP, 

Base-flow percentage (ratio of base flow to streamflow multiplied by 100); --, not applicable]

80
Streamflow Base flow Runoff

Year count Date CFS in CFS in CFS in BFP

Station 01467031 Pennypack Creek at Horsham; drainage area 3.5 square miles

1

2

3

4

2011

2012

2013

2014

10.4

2.7

4.8

5.6

40.48

10.50

18.88

21.82

4.52

1.36

2.36

2.41

17.63

5.30

9.20

9.39

5.86

1.33

2.48

3.19

22.85

5.20

9.67

12.43

43.6

50.5

48.8

43.0

Station 01467036 Pennypack tributary at Hatboro; drainage area 4.4 square miles

1

2

3

4

2011

2012

2013

2014

11.4

3.5

5.6

5.1

35.41

10.78

17.55

16.00

4.43

1.17

1.95

1.89

13.80

3.65

6.08

5.88

6.94

2.28

3.68

3.25

21.61

7.12

11.47

10.12

39.0

33.9

34.7

36.8

Station 01467039 Pennypack Creek at Willow Grove; drainage area 22.2 square miles

1

2

3

4

2011

2012

2013

2014

72.5

32.6

41.8

45.3

44.33

19.96

25.56

27.68

35.20

18.30

21.73

23.61

21.52

11.22

13.29

14.43

37.30

14.25

20.07

21.67

22.81

8.74

12.27

13.25

48.6

56.2

52.0

52.1

Station 014670413 Pennypack Creek at Beth Ayres; drainage area 35.3 square miles

1

2

3

4

5

2010

2011

2012

2013

2014

71.1

111.2

49.9

61.8

68.5

27.33

42.77

19.24

23.78

26.32

37.54

43.36

29.20

33.51

37.28

14.44

16.68

11.26

12.89

14.33

33.53

67.86

20.71

28.33

31.18

12.89

26.10

7.98

10.89

11.99

52.8

39.0

58.5

54.2

54.5

Station 01467042 PennypackCreek at Pine Road; drainage area 37.9 square miles

1

2

3

4

5

6

7

8

9

2008

2009

2010

2011

2012

2013

2014

2015

2016

61.2

80.7

79.8

126.3

53.0

63.7

72.2

53.5

42.2

22.00

28.91

28.57

45.22

19.04

22.81

25.85

19.17

15.14

32.16

38.34

41.03

50.98

31.34

35.29

39.12

28.95

25.28

11.55

13.73

14.70

18.26

11.26

12.64

14.01

10.37

9.06

29.08

42.37

38.75

75.28

21.68

28.39

33.06

24.56

16.94

10.44

15.18

13.88

26.96

7.78

10.17

11.84

8.80

6.08

52.5

47.5

51.4

40.4

59.1

55.4

54.2

54.1

59.9

Table 1.2. Annual base flow (calendar year January December) as determined by hydrograph separation 

methods for available complete years of daily streamflow values at U.S. Geological Survey stream gages 

on Pennypack Creek, 2008-16. HYSEP local-minimum hydrograph separation method in USGS Groundwater 

Toolbox (Barlow and others, 2017).

[CFS, annual average flow, in cubic feet per second; in, annual average flow per unit area, in inches; BFP, Baseflow percentage (ratio of base flow to streamflow multiplied by 100)]

81
Year count

Calendar

year

Streamflow

CFS in

Base flow

CFS in

Runoff

CFS in BFP

1 2001 69.1 22.93 34.25 11.37 34.86 11.57 49.6

2 2002 50.3 16.68 24.24 8.04 26.03 8.64 48.2

3 2003 115.8 38.43 48.85 16.21 66.94 22.22 42.2

4 2004 97.1 32.30 43.50 14.48 53.56 17.83 44.8

5 2005 76.8 25.49 38.36 12.73 38.44 12.76 49.9

6 2006 93.5 31.04 40.64 13.49 52.89 17.55 43.5

7 2007 81.9 27.18 34.88 11.58 47.02 15.61 42.6

8 2008 66.5 22.14 28.59 9.51 37.95 12.63 43.0

9 2009 97.9 32.50 42.26 14.03 55.65 18.47 43.2

10 2010 83.5 27.73 37.82 12.55 45.72 15.17 45.3

11 2011 130.9 43.43 49.62 16.47 81.24 26.96 37.9

12 2012 58.4 19.43 30.87 10.28 27.49 9.15 52.9

13 2013 78.2 25.95 38.83 12.89 39.36 13.06 49.7

14 2014 87.8 29.13 41.78 13.87 45.99 15.26 47.6

15 2015 57.8 19.17 29.16 9.68 28.59 9.49 50.5

16 2016 48.0 15.96 27.23 9.06 20.74 6.90 56.8

Table 1.3. Annual base flow (calendar year January−December) as determined by hydrograph separation 

methods for daily streamflow values at U.S. Geological Survey stream gage 01473900, Wissahickon Creek at 

Fort Washington, 2001−16, drainage area 40.9 square miles. HYSEP local-minimum hydrograph separation 

method in USGS Groundwater Toolbox (Barlow and others, 2017).

[CFS, annual average flow, in cubic feet per second; in, annual average flow per unit area, in inches; BFP, 

Base-flow percentage (ratio of base flow to streamflow multiplied by 100)]

82
Table 1.4. Regressions with monthly base flow at U.S. Geological Survey gage 01467042 Pennypack Creek at Pine Road, at Philadelphia, Pennysylvania to fill missing monthly data for January−March 

2010 and October−December 2015 and estimate annual average 2010 and 2015 base flow at four Pennypack Creek USGS continuous gages, 2010−15. Base flow determined using daily streamflow 

values and HYSEP local-minimum hydrograph separation method in USGS Groundwater Toolbox (Barlow and others, 2017). Estimated values in bold italics. Annual average base flow computed from 

average of monthly average base-flow values. Station headings shaded in colors shown for sites in regression plots (fig. 1.1).

2 [M-YY, month-year; DA, drainage area; mi , square miles; CFS, in cubic feet per second; BF, base flow; BFP, Base-flow percentage (ratio of base flow to streamflow multiplied by 100)]

01467031 Pennypack 01467036 Pennypack tributary 01467039 Pennypack Creek 014670413 Pennypack Creek 

01467042 Pennypack Creek at Pine Road Creek at Horsham at Hatboro at Willow Grove at Bethayres

(DA 37.9 mi2) (DA 3.5 mi2

) (DA 4.4 mi2

) (DA 22.2 mi2) (DA 35.3 mi2

)

Equation1 Equation1 Equation1 Equation 

Monthly Monthly Monthly Annual Monthly used to Annual Monthly used to Annual Monthly used to Annual Monthly used to 

average average average average estimate average average estimate average average estimate average average estimate average Date stream flow base flow runoff Monthly base flow base flow monthly BF at base flow base flow monthly BF at base flow base flow monthly BF at base flow base flow monthly BF at 

(M-YY) (CFS) (CFS) (CFS) average BFP Year (CFS) (CFS) Hatboro (CFS) (CFS) Horsham (CFS) (CFS) Willow Grove (CFS) (CFS) Beth Ayres

J-10 102.0 53.25 48.78 52.2 3.03 [1] 4.02 [2] a 31.16 [3] 50.79

F-10 97.2 65.29 31.88 67.2 3.71 [1] 5.37 [2] a 38.08 [3] 60.36

M-10 250.9 94.41 156.50 37.6 5.35 [1] 8.63 [2] a 55.63 94.11

A-10 89.1 65.75 23.33 73.8 3.68 5.31 37.44 70.61

M-10 52.2 40.92 11.31 78.3 2.01 2.16 22.50 35.12

J-10 37.9 25.91 12.01 68.3 1.09 0.84 15.45 22.82

J-10 81.3 25.04 56.27 30.8 1.21 0.55 15.30 23.65

A-10 32.3 19.82 12.49 61.3 0.89 0.33 11.94 17.21

S-10 29.3 18.56 10.74 63.3 0.75 0.10 10.84 13.55

O-10 86.6 29.85 56.80 34.5 2.24 0.91 17.75 23.02

N-10 36.5 25.44 11.07 69.7 2.19 0.71 16.34 18.06

D-10 59.5 29.57 29.95 49.7 2010 2.34 1.97 2.51 1.14 24.21 18.11 37.65 22.55

J-11 43.3 35.15 8.15 81.2 1.48 1.20 18.46 26.56

F-11 82.2 55.43 26.81 67.4 5.13 4.17 38.34 50.52

M-11 137.9 70.55 67.37 51.2 5.78 6.82 46.01 63.17

A-11 133.0 64.01 69.00 48.1 3.66 4.83 37.34 55.13

M-11 78.7 58.47 20.26 74.3 3.34 3.97 34.00 37.79

J-11 46.9 34.20 12.67 73.0 2.71 1.55 19.54 31.48

J-11 42.4 19.05 23.30 45.0 0.66 0.39 12.76 18.35

A-11 321.4 34.08 287.28 10.6 3.77 5.03 32.08 37.27

S-11 306.4 70.80 235.64 23.1 14.98 11.05 87.37 63.15

O-11 82.0 54.41 27.61 66.3 4.61 3.95 31.89 44.86

N-11 101.6 51.09 50.48 50.3 3.52 5.17 30.38 42.43

D-11 137.7 65.46 72.27 47.5 2011 4.46 3.88 4.53 6.25 35.32 35.62 43.47 50.95

1Equations relating base flow at gage with missing values (y) to base flow at Pennypack Creek at Pine Road (x): [4] y = 0.0753x - 0.7773; [5] y = 0.0925x - 0.8041; [6] y = 0.6264x +0.5025; [7] y = 0.9206x - 0.6814.

a

Regression equation for period April 2010 - August 2011 at Horsham; regression equation for period April 2010 - December 2010 of y = 0.1097x-2.0833 results in estimated monthly average base flow of

3.76 cfs for January 2010, 5.08 cfs for February 2010, 8.27 cfs for March 2010, with estimated annual average of 2.43 cfs at Horsham (3 percent smaller than estimated value of 2.51 cfs used for model input).

83
Table 1.4. Regressions with monthly base flow at U.S. Geological Survey gage 01467042 Pennypack Creek at Pine Road, at Philadelphia, Pennysylvania to fill missing monthly data for January−March 

2010 and October−December 2015 and estimate annual average 2010 and 2015 base flow at four Pennypack Creek USGS continuous gages, 2010−15. Base flow determined using daily streamflow 

values and HYSEP local-minimum hydrograph separation method in USGS Groundwater Toolbox (Barlow and others, 2017). Estimated values in bold italics. Annual average base flow computed from 

average of monthly average base-flow values.Station headings shaded in colors shown for sites in regression plots (fig. 1.1).--Continued

2 [M-YY, month-year; DA, drainage area; mi , square miles; CFS, cubic feet per second; BF, base flow; BFP, Base-flow percentage (ratio of base flow to streamflow multiplied by 100)]

01467042 Pennypack Creek at Pine Road

(DA 37.9 mi2)

 01467031 Pennypack 

Creek at Horsham 

(DA 3.5 mi2)

01467036 Pennypack tributary 

at Hatboro

(DA 4.4 mi2)

01467039 Pennypack Creek 

at Willow Grove

(DA 22.2 mi2)

014670413 Pennypack Creek 

at Bethayres

(DA 35.3 mi2)

 Date

(M-YY)

Monthly Monthly Monthly 

average average average 

stream flow base flow Monthly runoff 

(CFS) (CFS) (CFS) average BFP Year

Equation 

Annual Monthly used to 

average average estimate 

base flow base flow monthly BF at 

(CFS) (CFS) Hatboro

Equation 

Annual Monthly used to 

average average estimate 

base flow base flow monthly BF at 

(CFS) (CFS) Horsham

Equation 

Annual Monthly used to 

average average estimate 

base flow base flow monthly BF at 

(CFS) (CFS) Willow Grove

Equation 

Annual Monthly used to 

average average estimate 

base flow base flow monthly BF at 

(CFS) (CFS) Beth Ayres

J-12 92.8 59.53 33.28 64.1 2.82 4.66 33.67 48.37

F-12 54.8 49.48 5.35 90.2 1.97 2.45 27.16 37.54

M-12 45.5 36.96 8.58 81.2 1.52 1.76 22.22 33.74

A-12 49.6 28.76 20.82 58.0 1.09 1.21 16.19 30.61

M-12 62.7 35.22 27.44 56.2 1.48 1.72 19.73 39.52

J-12 37.8 30.10 7.65 79.7 0.90 1.17 16.20 29.64

J-12 39.3 21.37 17.91 54.4 0.53 0.33 11.21 20.04

A-12 48.9 22.25 26.67 45.5 1.01 0.37 14.28 21.05

S-12 57.6 23.42 34.22 40.6 0.58 0.36 14.60 22.15

O-12 48.6 21.27 27.37 43.7 0.44 0.33 12.11 20.45

N-12 25.1 21.49 3.64 85.5 0.32 0.64 13.95 20.36

D-12 72.2 26.74 45.47 37.0 2012 1.17 1.38 1.36 1.28 18.32 18.50 29.21 27.05

J-13 65.9 37.26 28.62 56.6 1.68 1.94 21.05 37.87

F-13 51.9 41.07 10.86 79.1 2.45 2.78 25.72 37.13

M-13 63.0 43.13 19.85 68.5 2.55 3.47 27.95 40.69

A-13 52.8 38.48 14.35 72.8 2.19 3.07 23.42 35.07

M-13 52.3 35.92 16.43 68.6 1.38 2.64 22.03 33.20

J-13 136.5 50.94 85.59 37.3 3.44 4.18 31.64 49.16

J-13 62.0 34.23 27.75 55.2 2.46 2.24 21.78 35.74

A-13 76.0 31.32 44.73 41.2 1.91 1.97 18.30 30.60

S-13 35.5 25.00 10.46 70.5 1.03 1.43 14.21 21.69

O-13 41.1 23.77 17.29 57.9 0.82 1.00 14.25 20.05

N-13 52.4 21.47 30.91 41.0 0.77 0.84 12.99 19.41

D-13 74.3 41.25 33.08 55.5 2013 1.96 2.79 2.36 2.79 21.76 27.73 33.52 41.61

84
Table 1.4. Regressions with monthly base flow at U.S. Geological Survey gage 01467042 Pennypack Creek at Pine Road, at Philadelphia, Pennysylvania to fill missing monthly data for January−March 

2010 and October−December 2015 and estimate annual average 2010 and 2015 base flow at four Pennypack Creek USGS continuous gages, 2010−15. Base flow determined using daily streamflow 

values and HYSEP local-minimum hydrograph separation method in USGS Groundwater Toolbox (Barlow and others, 2017). Estimated values in bold italics. Annual average base flow computed from 

average of monthly average base-flow values. Station headings shaded in colors shown for sites in regression plots (fig. 1.1).--Continued

2 [M-YY, month-year; DA, drainage area; mi , square miles; CFS, in cubic feet per second; BF, base flow; BFP, Base-Flow Percentage (ratio of base flow to streamflow multiplied by 100)]

01467042 Pennypack Creek at Pine Road

(DA 37.9 mi2)

 01467031 Pennypack 

Creek at Horsham 

(DA 3.5 mi2)

01467036 Pennypack tributary 

at Hatboro

(DA 4.4 mi2)

01467039 Pennypack Creek 

at Willow Grove

(DA 22.2 mi2)

014670413 Pennypack Creek 

at Beth ayres

(DA 35.3 mi2)

 Date 

(M-YY)

Monthly Monthly Monthly 

average average average 

stream flow base flow runoff Monthly 

(CFS) (CFS) (CFS) average BFP Year

Equation1

Annual Monthly used to 

average average estimate 

base flow base flow monthly BF at 

(CFS) (CFS) Hatboro

Equation1

Annual Monthly used to 

average average estimate 

base flow base flow monthly BF at 

(CFS) (CFS) Horsham

Equation1

Annual Monthly used to 

average average estimate 

base flow base flow monthly BF at 

(CFS) (CFS) Willow Grove

Equation1

Annual Monthly used to 

average average estimate 

base flow base flow monthly BF at 

(CFS) (CFS) Beth Ayres

J-14 79.4 53.52 25.90 67.4 3.54 3.72 33.57 52.75

F-14 112.3 52.10 60.22 46.4 3.26 4.15 38.36 47.98

M-14 109.1 56.41 52.72 51.7 3.68 5.71 38.58 50.07

A-14 148.4 64.51 83.93 43.5 3.86 4.29 37.59 64.21

M-14 138.0 66.91 71.04 48.5 3.18 3.54 33.88 59.97

J-14 69.2 43.86 25.35 63.4 1.07 2.55 22.38 45.26

J-14 47.6 28.49 19.13 59.8 0.45 1.31 14.49 27.82

A-14 29.3 21.34 7.95 72.9 0.23 0.65 11.70 19.82

S-14 19.6 17.12 2.48 87.3 0.21 0.39 9.63 16.13

O-14 27.9 18.28 9.61 65.5 0.40 0.51 9.89 17.53

N-14 41.1 21.51 19.61 52.3 1.07 0.81 12.77 20.44

D-14 47.7 26.40 21.34 55.3 2014 1.90 1.80 2.42 1.36 23.69 21.46 37.35 26.27

J-15 71.5 32.33 39.17 45.2 1.87 1.78 22.18 29.19

F-15 55.6 40.54 15.07 72.9 2.24 2.35 20.59 32.62

M-15 121.7 60.96 60.76 50.1 4.41 5.42 42.83 58.75

A-15 70.2 43.57 26.60 62.1 2.09 3.20 27.02 39.76

M-15 34.6 32.32 2.30 93.4 1.10 2.06 17.78 26.14

J-15 50.9 24.92 25.94 49.0 0.48 1.35 15.41 20.70

J-15 35.3 19.66 15.69 55.6 0.38 1.14 13.71 19.19

A-15 37.9 13.89 24.06 36.6 0.57 0.78 11.81 14.30

S-15 35.3 14.35 20.95 40.7 1.15 0.83 10.17 13.32

O-15

N-15

D-15

35.9

32.2

60.5

17.91

16.57

30.99

17.95

15.60

29.47

49.9

51.5

51.3 2015 1.41

0.57

0.47

1.56

[4]

[4]

[4] 1.88

0.85

0.73

2.06

[5]

[5]

[5] 18.67

11.72

10.88

19.91

[6]

[6]

[6] 26.02

15.81

14.57

27.85

[7]

[7]

[7]

1Equations relating base flow at gage with missing values (y) to base flow at Pennypack Creek at Pine Road (x): [4] y = 0.0753x - 0.7773; [5] y = 0.0925x - 0.8041; [6] y = 0.6264x +0.5025; [7] y = 0.9206x - 0.6814.

85
Measurement date 

and time

Time zone 

code

Gage height, in 

feet

Iron Works 

instantaneous 

discharge, in CFS

Horsham mean 

daily flow, in CFS

5/18/2010 13:25 EDT 1.83 9.3 8.32

8/17/2010 13:58 EDT 1.54 0.37 1.52

10/5/2010 13:05 EDT 1.66 2.38 5.06

2/7/2011 14:30 EST 1.75 4.54 6.69

5/11/2011 12:18 EDT 1.70 2.77 2.2

10/6/2011 14:10 EDT 1.73 4.15 4.26

2/6/2012 10:43 EST 1.72 3.78 3.39

5/16/2012 10:04 EDT 1.69 2.88 8.08

8/27/2012 12:51 EDT 1.84 7.76 1.61

10/10/2012 12:43 EDT 1.62 1.12 0.34

2/4/2013 16:07 EST 1.72 4.73 2.98

6/26/2013 12:35 EDT 1.70 3.48 4.56

10/25/2013 11:07 EDT 1.61 1.12 1.01

11/19/2013 14:49 EST 1.59 0.71 0.66

4/7/2014 13:21 EDT 1.78 4.43 6.64

8/29/2014 13:45 EDT 1.60 0.22 0.45

10/28/2014 14:19 EDT 1.60 0.62 0.73

1/23/2015 12:51 EST 1.71 3.93 2.93

1/23/2015 13:29 EST 1.71 5.08 2.93

1/23/2015 13:47 EST 1.71 5.52 2.93

5/7/2015 13:04 EDT 1.68 3.04 3.06

9/23/2015 11:58 EDT 1.55 0.47 0.86

Table 1.5. Instantaneous discharge (flow) measurements at U.S. Geological Survey gage 

01465460 Iron Works Creek at Richboro and daily mean flow at gage 01467031 Pennypack 

Creek at Horsham, Pennsylvania, 2010−15 (U.S. Geological Survey, 2018).

[CFS, cubic feet per second; EST, eastern standard time; EDT, eastern daylight time]

86
Table 1.6. Estimated base flow at U.S. Geological Survey (USGS) gage 01465460 Iron Works Creek at Richboro, Pennsylvania from 

relation to base flow at USGS gage 01467031 Pennypack Creek at Horsham, Pennsylvania, 2010−15. Base flow determined using 

daily streamflow values and HYSEP local-minimum hydrograph separation method in USGS Groundwater Toolbox (Barlow and others, 

2017). Estimated values in italics. Annual average base flow computed from average of monthly average base-flow values.

[DA, drainage area; mi2, square miles; CFS, cubic feet per second; BF, base flow; BFP, Base-flow percentage (ratio of base flow 

to streamflow multiplied by 100)]

 Date 

(M-YY)

01467042 Pennypack Creek at Pine Road 

(DA 37.9 mi2)

Monthly Monthly Monthly 

average average average 

stream flow base flow runoff Monthly 

(CFS) (CFS) (CFS) average BFP Year

 01467031 Pennypack Creek at 

Horsham 

(DA 3.5 mi2)

Annual Monthly Equation1 used 

average average to estimate 

base flow base flow monthly BF at 

(CFS) (CFS) Hatboro

01465460 Iron Works Creek at 

Richboro

(DA 3.69 mi2)

Estimated 

Equation1 used Annual monthly to estimate 

average average monthly BF at base flow base flow Iron Works

(CFS) (CFS)

J-10 102.0 53.25 48.78 52.2 3.03 [1] 2.85 [8]

F-10 97.2 65.29 31.88 67.2 3.71 [1] 3.38 [8]

M-10 250.9 94.41 156.50 37.6 5.35 [1] 4.62 [8]

A-10 89.1 65.75 23.33 73.8 3.68 3.36 [8]

M-10 52.2 40.92 11.31 78.3 2.01 2.01 [8]

J-10 37.9 25.91 12.01 68.3 1.09 1.19 [8]

J-10 81.3 25.04 56.27 30.8 1.21 1.30 [8]

A-10 32.3 19.82 12.49 61.3 0.89 1.00 [8]

S-10 29.3 18.56 10.74 63.3 0.75 0.87 [8]

O-10 86.6 29.85 56.80 34.5 2.24 2.20 [8]

N-10

D-10

36.5

59.5

25.44

29.57

11.07

29.95

69.7

49.7 2010 2.34

2.19

1.97 2.24

2.16

1.97

[8]

[8]

J-11 43.3 35.15 8.15 81.2 1.48 1.55 [8]

F-11 82.2 55.43 26.81 67.4 5.13 4.46 [8]

M-11 137.9 70.55 67.37 51.2 5.78 4.93 [8]

A-11 133.0 64.01 69.00 48.1 3.66 3.34 [8]

M-11 78.7 58.47 20.26 74.3 3.34 3.09 [8]

J-11 46.9 34.20 12.67 73.0 2.71 2.59 [8]

J-11 42.4 19.05 23.30 45.0 0.66 0.78 [8]

A-11 321.4 34.08 287.28 10.6 3.77 3.43 [8]

S-11 306.4 70.80 235.64 23.1 14.98 11.09 [8]

O-11 82.0 54.41 27.61 66.3 4.61 4.07 [8]

N-11 101.6 51.09 50.48 50.3 3.52 3.23 [8]

D-11 137.7 65.46 72.27 47.5 2011 4.46 3.88 3.84 3.51 [8]

J-12 92.8 59.53 33.28 64.1 2.82 2.68 [8]

F-12 54.8 49.48 5.35 90.2 1.97 1.97 [8]

M-12 45.5 36.96 8.58 81.2 1.52 1.58 [8]

A-12 49.6 28.76 20.82 58.0 1.09 1.19 [8]

M-12 62.7 35.22 27.44 56.2 1.48 1.55 [8]

J-12 37.8 30.10 7.65 79.7 0.90 1.01 [8]

J-12 39.3 21.37 17.91 54.4 0.53 0.65 [8]

A-12 48.9 22.25 26.67 45.5 1.01 1.12 [8]

S-12 57.6 23.42 34.22 40.6 0.58 0.70 [8]

O-12 48.6 21.27 27.37 43.7 0.44 0.55 [8]

N-12 25.1 21.49 3.64 85.5 0.32 0.42 [8]

D-12 72.2 26.74 45.47 37.0 2012 1.17 1.38 1.24 1.46 [8]

1 See table 3.4 for Equations [1] and [4] relating base flow at gage with missing values (x) to base flow at Pennypack Creek at Pine Road. 

Equation [8] estimating monthly base flow at Iron Works gage (y) using monthly base flow at Pennypack Creek at Horsham (x): 

 y = 1.1084x0.8509 , where equation based on regression of instantaneous flow at Iron Works (IW) with Horsham daily values (dv) in CFS;

for v = dv Horsham, w = IW instantaneous flow, w = 1.1084v0.8509 , R2 = 0.5863; see table 1.5 and figure 1.2.

87
Table 1.6. Estimated base flow at U.S. Geological Survey (USGS) gage 01465460 Iron Works Creek at Richboro, Pennsylvania from 

relation to base flow at USGS gage 01467031 Pennypack Creek at Horsham, Pennsylvania, 2010−15. Base flow determined using 

daily streamflow values and HYSEP local-minimum hydrograph separation method in USGS Groundwater Toolbox (Barlow and others, 

2017). Estimated values in italics. Annual average base flow computed from average of monthly average base-flow values. --Continued

2 [DA, drainage area; mi , square miles; CFS, cubic feet per second; BF, base flow; BFP, Base-flow percentage (ratio of base flow 

to streamflow multiplied by 100)]

01467042 Pennypack Creek at Pine Road 

(DA 37.9 mi2

)

 01467031 Pennypack Creek at 

Horsham 

(DA 3.5 mi2)

01465460 Iron Works Creek at 

Richboro

(DA 3.69 mi2)

 Date 

(M-YY)

Monthly Monthly Monthly 

average average average 

stream flow base flow runoff Monthly 

(CFS) (CFS) (CFS) average BFP Year

Annual Monthly Equation1 used 

average average to estimate 

base flow base flow monthly BF at 

(CFS) (CFS) Hatboro

Estimated Equation1 used 

Annual monthly to estimate 

average average monthly BF at 

base flow base flow Iron Works

(CFS) (CFS)

J-13 65.9 37.26 28.62 56.6 1.68 1.72 [8]

F-13 51.9 41.07 10.86 79.1 2.45 2.38 [8]

M-13 63.0 43.13 19.85 68.5 2.55 2.46 [8]

A-13 52.8 38.48 14.35 72.8 2.19 2.16 [8]

M-13 52.3 35.92 16.43 68.6 1.38 1.46 [8]

J-13 136.5 50.94 85.59 37.3 3.44 3.17 [8]

J-13 62.0 34.23 27.75 55.2 2.46 2.38 [8]

A-13 76.0 31.32 44.73 41.2 1.91 1.92 [8]

S-13 35.5 25.00 10.46 70.5 1.03 1.14 [8]

O-13 41.1 23.77 17.29 57.9 0.82 0.94 [8]

N-13 52.4 21.47 30.91 41.0 0.77 0.89 [8]

D-13 74.3 41.25 33.08 55.5 2013 1.96 2.79 1.94 2.65 [8]

J-14 79.4 53.52 25.90 67.4 3.54 3.25 [8]

F-14 112.3 52.10 60.22 46.4 3.26 3.03 [8]

M-14 109.1 56.41 52.72 51.7 3.68 3.36 [8]

A-14 148.4 64.51 83.93 43.5 3.86 3.50 [8]

M-14 138.0 66.91 71.04 48.5 3.18 2.97 [8]

J-14 69.2 43.86 25.35 63.4 1.07 1.17 [8]

J-14 47.6 28.49 19.13 59.8 0.45 0.56 [8]

A-14 29.3 21.34 7.95 72.9 0.23 0.32 [8]

S-14 19.6 17.12 2.48 87.3 0.21 0.29 [8]

O-14 27.9 18.28 9.61 65.5 0.40 0.51 [8]

N-14 41.1 21.51 19.61 52.3 1.07 1.17 [8]

D-14 47.7 26.40 21.34 55.3 2014 1.90 1.80 1.83 1.83 [8]

J-15 71.5 32.33 39.17 45.2 1.87 1.89 [8]

F-15 55.6 40.54 15.07 72.9 2.24 2.20 [8]

M-15

A-15

121.7

70.2

60.96

43.57

60.76

26.60

50.1

62.1

4.41

2.09

3.92

2.08

[8]

[8]

M-15 34.6 32.32 2.30 93.4 1.10 1.20 [8]

J-15 50.9 24.92 25.94 49.0 0.48 0.59 [8]

J-15 35.3 19.66 15.69 55.6 0.38 0.49 [8]

A-15 37.9 13.89 24.06 36.6 0.57 0.69 [8]

S-15 35.3 14.35 20.95 40.7 1.15 1.25 [8]

O-15 35.9 17.91 17.95 49.9 0.57 [4] 0.69 [8]

N-15

D-15

32.2

60.5

16.57

30.99

15.60

29.47

51.5

51.3 2015 1.41

0.47

1.56

[4]

[4] 1.43

0.58

1.61

[8]

[8]

88

1 See table 3.4 for Equations [1] and [4] relating base flow at gage with missing values (x) to base flow at Pennypack Creek at Pine Road. 

Equation [8] estimating monthly base flow at Iron Works gage (y) using monthly base flow at Pennypack Creek at Horsham (x): 

 y = 1.1084x0.8509 , where equation based on regression of instantaneous flow at Iron Works (IW) with Horsham daily values (dv) in CFS;

for v = dv Horsham, w = IW instantaneous flow, w = 1.1084v0.8509 , R2 = 0.5863; see table 1.5 and figure 1.2.
89 

Figure 1.1. Regressions with (A) 2010 monthly base flow at U.S. Geological Survey gage Pennypack 

Creek at Pine Road and three other gages on Pennypack Creek: Horsham, Hatboro, Willow Grove; and 

(B) 2015 monthly base flow at U.S. Geological Survey gage Pennypack Creek at Pine Road and four 

other gages on Pennypack Creek: Horsham, Hatboro, Willow Grove, Bethayres. 
90 

Figure 1.2. Regression between daily mean streamflow at gage 01467031 Pennypack Creek at Horsham and 

instantaneous streamflow at gage 01465460 Iron Works Creek at Richboro, 2010–15.

Figure 1.3. Instantaneous streamflow at gage 01465460 Iron Works Creek at Richboro and continuous 15-minute streamflow at gage 01467031 

Pennypack Creek at Horsham, 2010–15.

y = 1.1084x0.8509

R2 = 0.5863

0

2

4

6

8

10

0 1 2 3 4 5 6 7 8 9

Iron Works instantaneous flow, 

in cubic feet per second

Horsham mean daily flow, in cubic feet per second
91

Appendix 2. Model Calibration Results

PEST++ (Welter and others, 2015), which is based on the methods and algorithms of 

PEST (Doherty, 2015), was used for model calibration. Python scripts for PEST++ calibration of 

the model using FloPy (Bakker and others, 2016, 2018) and other Python libraries are provided 

in the model data release (Goode and Senior, 2020), along with all MODFLOW 6 input and 

output files. This appendix presents additional details on the model calibration results. 

Spatial Distribution of Water-Level Residuals

The water-level residuals, or differences between observed and modeled values, are 

shown for the study area for 1999, 2010, 2013, and 2016 are shown in figures 2.1–2.4. The 

modeled values are from the final calibration simulation. 

Parameter Sensitivities

Composite scaled sensitivities (fig. 2.5) measure the extent to which each parameter is 

‘sensitive’ to the measurements used for calibration (see White and others, 2019, p. 83). 

Parameters with high composite scaled sensitivity are relatively well constrained by the 

measurements, whereas parameters with low sensitivity are not. 
92 

Figure 2.1. Water-level (WL) residuals, differences between observed (Obs) and modeled (Model) values, in meters (m), for simulation of 1999 

conditions. 
93 

a

Figure 2.2. Water-level (WL) residuals, differences between observed (Obs) and modeled (Model) values, in meters (m), for simulation of 2010 

conditions. 
94 

Figure 2.3. Water-level (WL) residuals, differences between observed (Obs) and modeled (Model) values, in meters (m), for simulation of 2013 

conditions. 
95 

Figure 2.4. Water-level (WL) residuals, differences between observed (Obs) and modeled (Model) values, in meters (m), for simulation of 2016 

conditions. 
96 

Figure 2.5. Composite scaled sensitivities for model parameters. Yes parameters were estimated in final calibrations; Tied parameters were adjusted 

proportionally to a related Yes parameter (table 10, main text); No parameters were fixed based on manual adjustments and test calibrations. 
97

Recharge Rate Estimates

All recharge values were scaled by a single model parameter. This parameter had the 

highest sensitivity to observations used for calibration (fig. 2.5). Recharge rates were highest in 

the Carbonate hydrogeologic zone and lowest in the Lockatong and Other zones (table 2.1, also 

see fig. 15). 

Table 2.1. Optimum and approximate, individual, 95-percent confidence-interval values for recharge rates 

from calibration of the groundwater-flow model near Willow Grove and Warminster, Pennsylvania.

[2010 and 2013 rates; rates for 1999 are 0.9 times the values shown; rates for 2016 and 2017 are 0.7 times the values 

shown; m/day, meters per day; in/yr, inches per year]

Hydrogeologic zone Optimum Lower 

bound

Upper 

bound Optimum Lower bound

Upper 

bound

m/day in/yr

Brunswick 0.000587 0.000573 0.000601 8.43 8.23 8.63

Carbonate 0.00117 0.001142 0.001198 16.81 16.41 17.21

Upper Stockton 0.000685 0.000669 0.000701 9.85 9.61 10.08

Lockatong 0.000587 0.000573 0.000601 8.43 8.23 8.63

Other 0.000587 0.000573 0.000601 8.43 8.23 8.63

Middle and lower Stockton 0.00088 0.000859 0.000901 12.65 12.35 12.95
98

Hydraulic Conductivity Estimates

The final parameter values reflect regional-scale estimates of bulk hydraulic conductivity 

(fig. 2.6, table 2.2) based on water-level measurement, stream base flow, and the conceptual 

model. The Carbonate zone has the highest estimated K (42 ft/d). For the Newark Basin strata, 

the middle and lower members of the Stockton Formation, which underlies both Willow Grove 

and Warminster has the highest Kx (2.5 ft/d), followed by the Lockatong Formation (2.3 ft/d) 

and lastly the upper member of the Stockton Formation (1.2 ft/d). Ky/Kx in the dip direction is 

essentially equal for the Lockatong Formation (0.1) and the upper member of the Stockton 

Formation (0.2). The Brunswick Group strata underlie only a small portion of the model domain 

and K was fixed based the low sensitivity to the calibration data. 

The upper and lower confidence bounds reflect the relative uncertainty in the parameter 

values, which is smaller for the more sensitive parameters. The recharge parameter had the 

highest sensitivity, based on the weight used for the streamflow data, and has a narrow 

confidence interval (table 2.2). The streambed conductance has a relatively wide confidence 

interval. 

Figure 2.6. Estimated hydraulic conductivity in the strike direction (Kx) at depths between 125 and 400 

feet below land surface in model area in layers 4 and 5. Bases shown in magenta outline.
99 

Table 2.2. Optimum and approximate, individual 95-percent confidence-interval values for streambed and aquifer hydraulic conductivity estimated 

through calibration of the groundwater-flow model near Willow Grove and Warminster, Pennsylvania.

[m/d, meters per day; ft/d, feet per day; K, hydraulic conductivity; Kx, K in strike direction; Ky, K in dip direction; Kz, K perpendicular to bedding; --, not 

applicable or not available]

Hydrogeologic zone Layers Method Metric 

units

Optimum 

value

Confidence bounds English 

units

Optimum 

value

Confidence bounds

Lower Upper Lower Upper

Streambed hydraulic conductivity

All model streams 1 Calibrate m/d 1.61 0.23 10.7 ft/d 5.28 0.75 35.1

Hydraulic conductivity (Kx, max or strike value)

Weathered rock and overburden 1,2 Fixed m/d 0.12 -- -- ft/d 0.4 -- --

Brunswick Group 3-7 Fixed m/d 0.80 -- -- ft/d 2.6 -- --

Lockatong Formation 3-7 Calibrate m/d 0.70 0.53 0.88 ft/d 2.3 1.7 2.9

Stockton Formation, upper 3-7 Tied1 m/d 0.38 0.36 0.4 ft/d 1.24 1.18 1.31

Stockton Formation, middle and lower 3-7 Calibrate m/d 0.758 0.719 0.80 ft/d 2.5 2.4 2.6

Carbonate 3-7 Calibrate m/d 12.8 9.8 16.7 ft/d 42 32 55

Other 3-7 Fixed m/d 0.106 -- -- ft/d 0.35 -- --

1

Upper Stockton Kx = Middle and lower Stockton Kx/2.

Horizontal anisotropy Vertical anisotropy

Hydrogeologic zone Layers Method

Optimum 

value 

(Ky/Kx)

Optimum 

value 

(Kx/Ky)

Layers Method

Optimum 

value 

(Ky/Kx)

Optimum 

value 

(Kx/Ky)

Weathered rock and overburden 1,2 Fixed 1 1 1,2 Fixed 0.33 3

Brunswick Group 3-7 Fixed 0.33 3 3-7 Fixed 0.1 10

Lockatong Formation 3-7 Tied2 0.1 10 3-7 Tied2 0.1 10

Stockton Formation, upper 3-7 Tied3 0.18 5.4 3-7 Tied2 0.1 10

Stockton Formation, middle and lower 3-7 Calibrate 0.18 5.4 3-7 Tied2 0.1 10

Carbonate 3-7 Tied2 1 1 3-7 Tied2 0.33 3

Other 3-7 Fixed 1 1 3-7 Fixed 0.33 3

2

Tied to Kx for the respective hydrogeologic zone.

3

Upper Stockton Ky/Kx = Middle and lower Stockton Ky/Kx.
100 

References 

Bakker, M., Post, V., Langevin, C.D., Hughes, J.D., White, J.T., Starn, J.J. and Fienen, M.N., 

2016, Scripting MODFLOW model development using Python and FloPy: Groundwater, v. 54, 

p. 733–739, accessed July 9, 2017, at https://doi.org/10.1111/gwat.12413.

Bakker, M., Post, V., Langevin, C.D., Hughes, J.D., White, J.T., Starn, J.J., and Fienen, M.N., 

2018, FloPy—Python package for creating, running, and post-processing MODFLOW-based 

models (ver. 3.2.10): U.S. Geological Survey software release, accessed October 19, 2018, at

https://doi.org/10.5066/F7BK19FH. 

Doherty, J., 2015, Calibration and uncertainty analysis for complex environmental models: 

Brisbane, Australia, Watermark Numerical Computing, 227 p. 

Goode, D.J., and Senior, L.A., 2020, MODFLOW 6 and MODPATH 7 model data sets used to 

evaluate groundwater flow in the vicinity of Horsham and Warminster, Bucks and 

Montgomery Counties, Pennsylvania—Preliminary simulations for conditions in 1999, 2010, 

2013, 2016, and 2017: U.S. Geological Survey data release, 

https://doi.org/10.5066/P9K36P5S. 

Welter, D.E., White, J.T., Hunt, R.J., and Doherty, J.E., 2015, Approaches in highly 

parameterized inversion—PEST++ Version 3, a Parameter ESTimation and uncertainty 

analysis software suite optimized for large environmental models: U.S. Geological Survey 

Techniques and Methods, book 7, section C12, 54 p. 
101 

Appendix 3. Simulated Water Levels and 

Groundwater-Flow Paths 

MODPATH 7 (Pollock, 2016, 2017) was used to calculate groundwater-flow paths to 

discharge areas under assumed steady-state flow for conditions in 1999, 2010, 2013, 2016, and 

2017. The overall water budget of the model aquifer consists of a balance between the recharge 

inflow, and some stream loss, and total outflows, which include base-flow discharge to streams 

and well withdrawals. The ultimate disposition of recharge in each model cell can be visualized 

from MODPATH particle tracking to show areas of recharge that discharges either to wells or to 

streams. All model input and output files, including MODFLOW 6 files, are provided in the 

model archive data release (Goode and Senior, 2020). 

Per- and polyfluoroalkyl substances (PFAS) contamination has been identified at and 

near the Willow Grove and Warminster bases. To help understand groundwater migration of 

contaminants, regional flow paths are produced from identified possible locations of PFAS 

contamination in groundwater to their ultimate discharges, to either streams or wells for the five 

time periods. 

Maps of simulated water levels, areas contributing recharge to wells and to streams, and 

regional groundwater-flow paths from possible sources are shown for the four calibration periods 

(1999, 2010, 2013, and 2016) and for the simulation of conditions in 2017. 
102 

1999 Conditions

Simulated water levels for the calibrated model in layer 1 for 1999 conditions are shown 

in figure 3.1. The 1999 conditions included the largest well withdrawals, so a higher proportion 

of recharge was captured by wells in this simulation (fig. 3.2), compared to the other time 

periods. These conditions used data from other dates in the 1990s. For example, withdrawal rates 

were from 1996 records (Delaware River Basin Commission, 2019). Regional groundwater-flow 

paths for 1999 conditions are shown in figures 3.3–3.5. 

Figure 3.1. ModelMuse screenshot of simulated potentiometric surface in model layer 1 for 1999 

conditions. The contour interval is 20 feet. 
103

Figure 3.2. ModelMuse screenshot of simulated areas that contribute recharge to discharging wells 

(green) and streams (white) for 1999 conditions. 
104 

Figure 3.3. ModelMuse screenshot of simulated regional groundwater-flow paths (black) from locations of 

possible per- and polyfluoroalkyl substances (PFAS) contamination in groundwater for 1999 conditions. 
105 

Figure 3.4. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 1999 

conditions near Willow Grove, with withdrawal wells in black. Pink squares represent locations of possible 

per- and polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from 

green for relatively short travel times to magenta for relatively long travel times. 
106 

Figure 3.5. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 1999 

conditions near Warminster, with withdrawal wells in black. Pink squares represent locations of possible 

per- and polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from 

green for relatively short travel times to magenta for relatively long travel times.
107

2010 Conditions

Simulated water levels for the calibrated model in layer 1 for 2010 conditions are shown 

in figure 3.6 and recharge simulated to be discharged to wells and streams is shown in figure 3.7. 

Willow Grove base production wells were in operation for this simulation. Recharge from 

precipitation was the highest in this period and in 2013. Regional groundwater-flow paths for 

2010 conditions are shown in figures 3.8–3.10.

Figure 3.6. ModelMuse screenshot of simulated potentiometric surface in model layer 1 for 2010 

conditions. The contour interval is 20 feet. 
108 

Figure 3.7. ModelMuse screenshot of simulated areas that contribute recharge to discharging wells 

(green) and streams (white) for 2010 conditions. 
109 

Figure 3.8. ModelMuse screenshot of simulated regional groundwater-flow paths from possible per- and 

polyfluoroalkyl substances (PFAS) contamination sources for 2010 conditions. Colors along paths range 

from green for relatively short travel times to magenta for relatively long travel times.
110 

Figure 3.9. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2010 

conditions near Willow Grove with withdrawal wells in black. Pink squares represent areas of possible perand polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from green 

for relatively short travel times to magenta for relatively long travel times. 
111 

Figure 3.10. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2010 

conditions near Warminster, with withdrawal wells in black. Pink squares represent areas of possible perand polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from green 

for relatively short travel times to magenta for relatively long travel times.
112 

2013 Conditions

Simulated water levels for the calibrated model in layer 1 for 2013 conditions are shown 

in figure 3.11 and recharge simulated to be discharged to wells and streams is shown in figure 

3.12. Willow Grove base production wells and many public supply wells in the area were not in 

operation. Recharge from precipitation was the highest in this period and in 2010. Regional 

groundwater-flow paths for 2013 conditions are shown in figures 3.13–3.15. 

Figure 3.11. ModelMuse screenshot of simulated potentiometric surface in model layer 1 for 2013 

conditions. The contour interval is 20 feet. 
113

Figure 3.12. ModelMuse screenshot of simulated areas that contribute recharge to discharging wells 

(green) and streams (white) for 2013 conditions. 
114 

Figure 3.13. ModelMuse screenshot of simulated regional groundwater-flow paths from potential per- and 

polyfluoroalkyl substances (PFAS) contamination sources for 2013 conditions. Colors along paths range 

from green for relatively short travel times to magenta for relatively long travel times.
115

Figure 3.14. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2013 

conditions near Willow Grove, with withdrawal wells in black. Pink squares represent areas of possible perand polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from green 

for relatively short travel times to magenta for relatively long travel times.
116 

Figure 3.15. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2013 

conditions near Warminster, with withdrawal wells in black. Pink squares represent areas of possible perand polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from green 

for relatively short travel times to magenta for relatively long travel times.
117

2016 Conditions

Simulated water levels for the calibrated model in layer 1 for 2016 conditions are shown 

in figure 3.16 and recharge simulated to be discharged to wells and streams is shown in figure 

3.17. Willow Grove base production wells and some public supply wells in the area were 

operating at lower rates than in 2010. Recharge from precipitation was the lowest during this 

period and in 2017. Regional groundwater-flow paths for 2016 conditions are shown in figures 

3.18–3.20.

Figure 3.16. ModelMuse screenshot of simulated potentiometric surface in model layer 1 for 2016 

conditions. Model streams are shown in blue gray. The contour interval is 20 feet.
118 

Figure 3.17. ModelMuse screenshot of simulated areas that contribute recharge to discharging wells 

(green) and streams (white) for 2016 conditions. 
119 

Figure 3.18. ModelMuse screenshot of simulated regional groundwater-flow paths from potential per- and 

polyfluoroalkyl substances (PFAS) contamination sources for 2016 conditions. Model streams are shown in 

blue gray. Colors along paths range from green for relatively short travel times to magenta for relatively 

long travel times.
120 

Figure 3.19. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2016 

conditions near Willow Grove, with withdrawal wells in black. Pink squares represent areas of possible perand polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from green 

for relatively short travel times to magenta for relatively long travel times. 
121 

Figure 3.20. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2016 

conditions near Warminster, with withdrawal wells in black. Pink squares represent areas of possible perand polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from green 

for relatively short travel times to magenta for relatively long travel times.
122 

2017 Conditions 

Simulated water levels for the calibrated model in layer 1 for 2017 conditions are shown 

in figure 3.21 and recharge simulated to be discharged to wells and streams is shown in figure 

3.22. This time period had the least groundwater withdrawals and the lowest recharge rate, which 

was equal to that in 2016. This time period was not used for model calibration. Regional 

groundwater-flow paths for 2017 conditions are shown in figures 3.23–3.25. 

Figure 3.21. ModelMuse screenshot of simulated potentiometric surface in model layer 1 for 2017 

conditions. Model streams are shown in blue gray. The contour interval is 20 feet. 
123 

Figure 3.22. ModelMuse screenshot of simulated areas that contribute recharge to discharging wells 

(green) and streams (white) for 2017 conditions. 
124 

Figure 3.23. ModelMuse screenshot of simulated regional groundwater-flow paths from potential per- and 

polyfluoroalkyl substances (PFAS) contamination sources for 2017 conditions. Model streams are shown in 

blue gray. Colors along paths range from green for relatively short travel times to magenta for relatively 

long travel times. 
125 

Figure 3.24. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2017 

conditions near Willow Grove, with withdrawal wells in black. Pink squares represent areas of possible perand polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from green 

for relatively short travel times to magenta for relatively long travel times. 
126 

Figure 3.25. Detail of ModelMuse screenshot of simulated regional groundwater-flow paths for 2017 

conditions near Warminster, with withdrawal wells in black. Pink squares represent areas of possible perand polyfluoroalkyl substances (PFAS) contamination in groundwater. Colors along paths range from green 

for relatively short travel times to magenta for relatively long travel times.
127 

References 

Delaware River Basin Commission, 2019, Southeastern Pennsylvania ground water protected 

area (GWPA) GIS data sets: Delaware River Basin Commission web page, accessed 

November 1, 2019, at https://www.nj.gov/drbc/programs/project/pr/gwpa-data.html. 

Goode, D.J., and Senior, L.A., 2020, MODFLOW 6 and MODPATH 7 model data sets used to 

evaluate groundwater flow in the vicinity of Horsham and Warminster, Bucks and 

Montgomery Counties, Pennsylvania—Preliminary simulations for conditions in 1999, 2010, 

2013, 2016, and 2017: U.S. Geological Survey data release, 

https://doi.org/10.5066/P9K36P5S. 

Pollock, D.W., 2016, User guide for MODPATH version 7—A particle-tracking model for 

MODFLOW: U.S. Geological Survey Open-File Report 2016–1086, 35 p., accessed July 9, 

2018, at https://doi.org/10.3133/ofr20161086.

Pollock, D.W., 2017, MODPATH v7.2.01—A particle-tracking model for MODFLOW: U.S. 

Geological Survey software release, 15 December 2017, https://doi.org/10.5066/F70P0X5X. 
For additional information, contact:

Director, Pennsylvania Water Science Center

U.S. Geological Survey

215 Limekiln Road

New Cumberland, PA 17070-2424

Or visit our website at:

https://www.usgs.gov/centers/pa-water

ISSN 2331-1258 (online)

https://doi.org/10.3133/ofr20191137