U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigations Report 2009–5039

Prepared in cooperation with the Hampton Roads Planning District Commission

Simulation of Groundwater Flow in the Coastal Plain

Aquifer System of Virginia
Cover. Landsat image of the Chesapeake Bay and adjacent coastal plain in Virginia and southern Maryland (photograph was 

contributed by the U.S. Geological Survey/National Aeronautics and Space Administration Landsat Program).
Simulation of Groundwater Flow in the 

Coastal Plain Aquifer System of Virginia

By Charles E. Heywood and Jason P. Pope

Prepared in cooperation with the Hampton Roads Planning District Commission

Scientific Investigations Report 2009–5039

U.S. Department of the Interior

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

KEN SALAZAR, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

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

For more information on the USGS—the Federal source for science about the Earth, its natural and living resources, 

natural hazards, and the environment, visit http://www.usgs.gov or call 1-888-ASK-USGS

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

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

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

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

U.S. Government.

Although this report is in the public domain, permission must be secured from the individual copyright owners to 

reproduce any copyrighted materials contained within this report.

Suggested citation:

Heywood, C.E., and Pope, J.P., 2009, Simulation of groundwater flow in the Coastal Plain aquifer system of Virginia: 

U.S. Geological Survey Scientific Investigations Report 2009–5039, 115 p.

ISBN 978–1–4113–2418–3
iii

Contents

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

Introduction.....................................................................................................................................................1

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

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

Description of the Virginia Coastal Plain ...................................................................................................6

Hydrogeologic Framework..................................................................................................................6

Hydrologic Conditions........................................................................................................................10

Precipitation................................................................................................................................11

Evapotranspiration ....................................................................................................................11

Recharge.....................................................................................................................................11

Groundwater Underflow from the Piedmont Province........................................................14

Hydraulic Properties .................................................................................................................14

Hydraulic Conductivity.....................................................................................................14

Specific Storage ...............................................................................................................15

Specific Yield.....................................................................................................................16

Discharge....................................................................................................................................16

Saltwater Transition Zone ........................................................................................................17

Simulation of Transient Groundwater Flow.............................................................................................19

Projection System...............................................................................................................................19

Numerical Method..............................................................................................................................19

Spatial Discretization .........................................................................................................................20

Internal Flow Package .......................................................................................................................21

Temporal Discretization.....................................................................................................................22

Boundary Conditions..........................................................................................................................22

Specified Flux Boundaries .......................................................................................................22

Transient Areal Recharge................................................................................................22

Transient Underflow to Maryland ..................................................................................24

Constant Underflow from Piedmont ..............................................................................24

Transient Groundwater Pumpage..................................................................................24

Reported Withdrawals............................................................................................24

Domestic Withdrawals ...........................................................................................26

Head-Dependent Flux Boundaries..........................................................................................26

Evapotranspiration ...........................................................................................................26

Rivers ................................................................................................................................26

Chesapeake Bay and Ocean Floor Leakage ................................................................26

Hydraulic-Property Distribution .......................................................................................................27

Horizontal Interpolation of Hydraulic Conductivity ..............................................................27

Depth-Dependence of Potomac Aquifer Horizontal Hydraulic Conductivity...................27
iv

Model Calibration.........................................................................................................................................29

Calibration Data Set ...........................................................................................................................30

Water-Level Observation Errors..............................................................................................33

Weights........................................................................................................................................33

Modification of First Stress Period Time-Discretization for Parameter Estimation.................34

Hydraulic Property Parameterization..............................................................................................34

Boundary Flux Parameterization......................................................................................................36

Simulated Conditions...................................................................................................................................36

Water Levels........................................................................................................................................36

Surficial Aquifer.........................................................................................................................37

Yorktown-Eastover Aquifer ......................................................................................................41

St. Marys Aquifer.......................................................................................................................45

Piney Point Aquifer....................................................................................................................49

Aquia Aquifer..............................................................................................................................53

Peedee Aquifer ..........................................................................................................................57

Virginia Beach Aquifer..............................................................................................................60

Potomac Aquifer........................................................................................................................63

Water Budget ......................................................................................................................................73

Steady-State Water Budget.....................................................................................................73

Transient Water Budget............................................................................................................74

Calibration Assessment..............................................................................................................................74

Residual Plots......................................................................................................................................74

Weighted Residual Maps ..................................................................................................................77

Parameter Correlation .......................................................................................................................84

Parameter Sensitivities......................................................................................................................85

Estimated Parameters........................................................................................................................85

Aquifer and Confining-Unit Spatial Variability ...............................................................................87

Model Limitations.........................................................................................................................................87

Summary......................................................................................................................................................105

Acknowledgments.....................................................................................................................................106

References Cited........................................................................................................................................106

Appendix......................................................................................................................................................110

Transient Saltwater-Transport Simulation....................................................................................110

Boundary Condition Modifications .......................................................................................111

Initial Conditions ......................................................................................................................111

Simulation of Mass Transport................................................................................................111

Appendix References Cited ............................................................................................................115

Figures

1. Map showing locations of counties and independent cities, the Chesapeake Bay 

impact crater, and other important geographic and physiographic features of the 

Virginia Coastal Plain ...................................................................................................................2

2. Graph showing Virginia reported and domestic groundwater withdrawals ......................3
v

3. Map showing finite-difference model grid with row and column numbers 

of cross sections...........................................................................................................................5

4. Vertical model cross sections showing distribution of hydrogeologic units from 

east to west along row 80 (A–A ́), east to west along row 95 (B–B ́), east to west 

along row 107 (C–C ́), east to west along row 120 (D–D ́), north to south along 

column 50 (E–E ́), north to south along column 60 (F–F ́), north to south along 

column 70 (G–G ́), and north to south along column 87 (H–H ́).............................................7

5. Stratigraphic correlations of hydrogeologic units of the Virginia Coastal Plain................9

6. Diagram showing generalized predevelopment groundwater flow patterns in the 

Coastal Plain aquifer system ....................................................................................................10

7. Map showing locations of precipitation-measurement stations, control points, 

and rivers simulated with the River package, other streams, and sea-level areas 

simulated with the General Head Boundary Package .........................................................12

8. Graph of historical measured precipitation at the four precipitation-measurement 

stations in the Virginia Coastal Plain.......................................................................................13

9. Pie chart showing estimated groundwater withdrawal-use categories in 2000 .............17

10. Map showing location of Coastal Plain reported withdrawal wells, model domain, 

and self-supplied withdrawal rate per square mile in 2000.................................................18

11. Graph showing simulated recharge in the Virginia Coastal Plain by year........................23

12. Vertical sections of finite-difference model grid, showing cells designated for timevarying specified flux to Maryland in row 1 and constant specified flux from the 

Piedmont Province in column 1................................................................................................25

13. Map showing pilot-point locations and simulated horizontal hydraulic conductivity 

at the top of the Potomac aquifer ............................................................................................28

14. Graph showing decrease with depth of the simulated horizontal hydraulic 

conductivity of the Potomac aquifer .......................................................................................29

15. Map showing locations and screened hydrogeologic unit of water-level 

observation wells........................................................................................................................31

16. Histogram showing time distribution and spatial distribution by State of water-level 

observations used for model calibration ................................................................................32

17. Map of estimated vertical hydraulic conductivity of the Calvert confining unit..............35

18. Contour maps of the simulated water table in the Virginia Coastal Plain:

a. before development ..........................................................................................................38

b. drawdown in 2001..............................................................................................................39

19. Map showing locations of observation wells screened in the surficial, 

Yorktown-Eastover, St. Marys, Piney Point, Aquia, and Peedee aquifers for 

hydrographs depicted in figures 20, 22, 24, 26, 28, and 30 ...................................................40

20. Hydrograph of observed and simulated water levels in the surficial aquifer 

well 57D 23 ...................................................................................................................................41

21. Contour maps of simulated water levels in the Yorktown-Eastover aquifer:

a. before development ..........................................................................................................42

b. in 2003 ..................................................................................................................................43

c. drawdown from predevelopment to 2003 ......................................................................44

22. Hydrograph of observed and simulated water levels in wells 57H 14, 61B 2, 62A 2, 

62C 9, and 57B 8 in the Yorktown-Eastover aquifer...............................................................45
vi

23. Contour maps of simulated water levels in the St. Marys aquifer:

a. before development ..........................................................................................................46

b. in 2003 ..................................................................................................................................47

c. drawdown from predevelopment to 2003 ......................................................................48

24. Hydrograph of observed and simulated water levels in well 57C 21 in the St. Marys 

aquifer...........................................................................................................................................49

25. Contour maps of simulated water levels in the Piney Point aquifer:

a. before development ..........................................................................................................50

b. in 2003 ..................................................................................................................................51

c. drawdown from predevelopment to 2003 ......................................................................52

26. Hydrograph of observed and simulated water levels in wells 56H 29, 54Q 73, and 

56G 38 in the Piney Point aquifer..............................................................................................53

27. Contour maps of simulated water levels in the Aquia aquifer:

a. before development ..........................................................................................................54

b. in 2003 ..................................................................................................................................55

c. drawdown from predevelopment to 2003 ......................................................................56

28. Hydrograph of observed and simulated water levels in wells 51M 18, 57G 2, 

and 55B 67 in the Aquia aquifer................................................................................................57

29. Contour maps of simulated water levels in the Peedee aquifer:

a. before development ..........................................................................................................58

b. in 2003 ..................................................................................................................................59

30. Hydrograph of observed and simulated water levels in well 61B 6 

in the Peedee aquifer.................................................................................................................60

31. Contour maps of simulated water levels in the Virginia Beach aquifer:

a. before development ..........................................................................................................61

b. in 2003 ..................................................................................................................................62

32. Contour maps of simulated water levels in the Potomac aquifer:

a. before development ..........................................................................................................64

b. in 2003 ..................................................................................................................................65

c. drawdown from predevelopment to 2003 ......................................................................66

33. Map showing locations of observation wells screened in the Potomac aquifer 

for hydrographs depicted in figures 34–39.............................................................................67

34. Hydrograph of observed and simulated water levels of selected wells screened 

in the Potomac aquifer in the southwest model area...........................................................68

35. Hydrograph of observed and simulated water levels of selected wells screened 

in the Potomac aquifer in the southeast model area............................................................68

36. Hydrograph of observed and simulated water levels of selected wells screened 

in the Potomac aquifer in the western model area...............................................................69

37. Hydrograph of observed and simulated water levels of selected wells screened 

in the Potomac aquifer in the eastern model area................................................................70

38. Hydrograph of observed and simulated water levels of selected wells screened 

in the Potomac aquifer in the central model area.................................................................71

39. Hydrograph of observed and simulated water levels of selected wells screened 

in the Potomac aquifer in the northern model area..............................................................72

40. Graphs of:

a. Relation between observed and simulated water levels............................................75

b. Relation between weighted observed and weighted simulated water levels ........76
vii

41. Graph of relation between weighted residuals and simulated water levels....................77

42. Maps showing maximum residual magnitudes for observation wells screened in the:

a. surficial aquifer..................................................................................................................78

b. Yorktown-Eastover aquifer...............................................................................................79

c. Piney Point aquifer ............................................................................................................80

d. Aquia aquifer......................................................................................................................81

e. Peedee and St. Marys aquifers.......................................................................................82

f. Potomac aquifer.................................................................................................................83

43. Graph of composite scaled sensitivities for model parameters .........................................85

44. Graph of estimated values, confidence intervals, and reasonable ranges for 

parameters representing horizontal and vertical hydraulic conductivity, specific 

storage, and other fluxes...........................................................................................................86

45. Map of estimated transmissivity of the:

a. Yorktown-Eastover aquifer...............................................................................................88

b. St. Marys aquifer ...............................................................................................................89

c. Piney Point aquifer ............................................................................................................90

d. Aquia aquifer......................................................................................................................91

e. Peedee aquifer...................................................................................................................92

f. Virginia Beach aquifer......................................................................................................93

46. Map of estimated leakance of the:

a. Yorktown confining zone ..................................................................................................94

b. St. Marys confining unit....................................................................................................95

c. Calvert confining unit........................................................................................................96

d. Chickahominy confining unit............................................................................................97

e. Exmore matrix confining unit...........................................................................................98

f. Exmore clast confining unit..............................................................................................99

g. Nanjemoy-Marlboro confining unit ..............................................................................100

h. Peedee confining zone ...................................................................................................101

i. Virginia Beach confining zone.......................................................................................102

j. Upper Cenomanian confining unit ................................................................................103

k. Potomac confining zone .................................................................................................104

A1. Graph of sea levels during the Pleistocene and Holocene epochs .................................110

A2. Map showing simulated water density near the saltwater transition zone 

of the Virginia Coastal Plain....................................................................................................112

A3. Cross section showing simulated water density near the saltwater transition 

zone of the Virginia Coastal Plain ..........................................................................................113

A4. Graph of relation between simulated and observed normalized chloride 

concentrations ..........................................................................................................................114

Tables

1. Summary of Potomac aquifer-performance test results......................................................15

2. Magnitudes of 2003 reported and 2000 self-supplied domestic groundwater 

withdrawals from the Virginia Coastal Plain by hydrogeologic unit..................................17

3. Model grid-layer thicknesses...................................................................................................21
viii

4. Time discretization for 113-year historical transient simulation.........................................22

5. Number of control points for simulated rivers.......................................................................27

6. Number of observations and wells primarily associated with each hydrogeologic 

unit simulated in the Virginia Coastal Plain model................................................................32

7. Steady-state water-budget components for model domain ...............................................73

8. Steady-state and transient water-budget components for Virginia sub-domain............73

9. Parameter pairs with absolute correlation coefficients greater than 0.8.........................84

A1. Specified Pleistocene-Epoch sea levels. .............................................................................111

A2. Values used for transport simulation.....................................................................................111

Conversion Factors and Datums

Inch/Pound to SI

Multiply By To obtain

Length

inch (in.) 2.54 centimeter (cm)

inch (in.) 25.4 millimeter (mm)

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Area

square mile (mi2) 259.0 hectare (ha)

square mile (mi2) 2.590 square kilometer (km2) 

Volume

cubic foot (ft3) 0.02832 cubic meter (m3) 

gallon (gal) 3.785 liter (L) 

gallon (gal) 0.003785 cubic meter (m3) 

Flow rate

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

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

Flux

inch per year (in/yr) 2.54 centimeter per year (cm/yr)

Mass

pound, avoirdupois (lb) 0.4536 kilogram (kg) 

Density

pound per cubic foot (lb/ft3) 16.02 kilogram per cubic meter (kg/m3)

pound per cubic foot (lb/ft3) 0.01602 gram per cubic centimeter (g/cm3)

Hydraulic conductivity

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

Transmissivity*

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

Temperature in degrees Fahrenheit (°F) may be converted to degrees Celsius (°C) as follows:

°C = (°F – 32) / 1.8

Vertical coordinate information is referenced to the National Geodetic Vertical Datum of 1929 

(NGVD 29).

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

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

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

foot of aquifer thickness [(ft3/d)/ft2]ft. In this report, the mathematically reduced form, foot 

squared per day (ft2/d), is used for convenience.

Acronyms and Abbreviations:

BCF Block-Centered Flow package

CSS Composite scaled sensitivity

DEM Digital elevation model

DEQ Department of Environmental Quality

ET evapotranspiration

EVT Evapotranspiration package

FHB Flow and Head Boundary package

GHB General Head Boundary package

GIS Geographic information system

GPS Global Positioning System

GWMA Ground Water Management Area

GWSI Groundwater Site Inventory

HRPDC Hampton Roads Planning District Commission

HUF Hydrogeologic-Unit Flow package

IMT Integrated Mass Transport

KDEP hydraulic-conductivity depth dependence

LPF Layer-Property Flow package

mg/L milligram per liter

MNW Multi-Node Well package

MSL mean sea level

PET potential evapotranspiration

RASA Regional Aquifer System Analysis

RCH Recharge package

RIV River package

RMSE Root mean square error

SSE sum-of-squares error

TDS total dissolved solids

USGS U.S. Geological Survey

VDF Variable-Density Flow
x

Parameter Names and Definitions:

et_max_rate Maximum evapotranspiration flux.

evapotransp Evapotranspiration extinction depth.

k_aquia Horizontal hydraulic conductivity of the Aquia aquifer.

k_aquiclud Horizontal hydraulic conductivity of the Exmore matrix, Exmore clast, and Chickahominy confining units.

k_aquitard Horizontal hydraulic conductivity of other confining units.

k_charlesci Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in Charles City.

k_chesapeak Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in Chesapeake County.

k_chesterfi Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in Chesterfield County.

k_columbia Horizontal hydraulic conductivity of the surficial aquifer.

kdepth1 Lambda exponent in equation 8.

k_franklin Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in the city of Franklin.

k_henrico Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in Henrico County.

k_islewight Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in Isle of Wight County. 

k_jamescity Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in James City County.

k_kinggeorg Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in King George County.

k_kingwilli Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in King William County.

k_lancaster Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in Lancaster County.

k_maryland Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in southern Maryland.

k_megabloc Horizontal hydraulic conductivity of the Potomac aquifer within the Chesapeake Bay Impact Crater.

k_newportne Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in the city of Newport News.

k_peedee_a Horizontal hydraulic conductivity of the Peedee aquifer.

k_pineyp_jc Horizontal hydraulic conductivity of the Piney Point aquifer.

k_potom_cz Horizontal hydraulic conductivity of the Potomac confing zone.

k_princegeo Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in Prince George County.

k_s_mary_a Horizontal hydraulic conductivity of the Saint Marys aquifer.

k_sohampt_n Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in northern Southhampton 

County.

k_sohampton Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in Southhampton County.

k_sussex Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in Sussex County.

k_virgin_a Horizontal hydraulic conductivity of the Virginia Beach aquifer.

k_westmore Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in Westmoreland County.

k_westwestm Horizontal hydraulic conductivity of the Potomac aquifer at the pilot point located in western Westmoreland 

County.

k_york_e_a Horizontal hydraulic conductivity of the Yorktown-Eastover aquifer.

mdflow Specified flux from finite-difference cells in row 1 shown in figure 12a.

piedflow Specified flux to finite-difference cells in columns 4–20 shown in figure 12b.

recharge Specified recharge flux.

river_cond Hydraulic conductance of simulated riverbeds.

rw Effective well radius in Multi-Node Well (MNW) package. 

spec_yield Specific yield.

ss_aquitard Specific storage of other confining units.

ss_other Specific storage of other aquifers.

ss_potomac Specific storage of the Potomac aquifer.

ss_upper_ce Specific storage of the Upper Cenomanian confining unit.

vk_aquifer Vertical hydraulic conductivity of other hydrogeologic units simulated as aquifers.

vk_calve_c Vertical hydraulic conductivity of the Calvert confining unit at other pilot point locations.

vk_chickah Vertical hydraulic conductivity of the Exmore matrix, Exmore clast, and Chickahominy confining units.

vk_james_c Vertical hydraulic conductivity of the Calvert confining unit at the pilot point located in James City county.

vk_nanjemo Vertical hydraulic conductivity of the Nanjemoy-Marlboro confining unit.

vk_peede_c Vertical hydraulic conductivity of the Peedee confining unit.

vk_poto_cz Vertical hydraulic conductivity of the Potomac confining zone.

vk_potomac Vertical hydraulic conductivity of the Potomac aquifer.

vk_smary_c Vertical hydraulic conductivity of the Saint Marys confining zone.

vk_surface Vertical hydraulic conductivity of the surficial aquifer.

vk_upper_c Vertical hydraulic conductivity of the Upper Cenomanian confining unit.

vk_virgb_c Vertical hydraulic conductivity of the Virginia Beach confining unit.

vk_york_cz Vertical hydraulic conductivity of the Yorktown-Eastover confining zone.
Simulation of Groundwater Flow in the Coastal Plain 

Aquifer System of Virginia

By Charles E. Heywood and Jason P. Pope

Abstract

The groundwater model documented in this report 

simulates the transient evolution of water levels in the 

aquifers and confining units of the Virginia Coastal Plain 

and adjacent portions of Maryland and North Carolina since 

1890. Groundwater withdrawals have lowered water levels in 

Virginia Coastal Plain aquifers and have resulted in drawdown 

in the Potomac aquifer exceeding 200 feet in some areas. The 

discovery of the Chesapeake Bay impact crater and a revised 

conceptualization of the Potomac aquifer are two major 

changes to the hydrogeologic framework that have 

been incorporated into the groundwater model. The spatial 

scale of the model was selected on the basis of the primary 

function of the model of assessing the regional water-level 

responses of the confined aquifers beneath the Coastal 

Plain. The local horizontal groundwater flow through the 

surficial aquifer is not intended to be accurately simulated. 

Representation of recharge, evapotranspiration, and interaction 

with surface-water features, such as major rivers, lakes, the 

Chesapeake Bay, and the Atlantic Ocean, enable simulation 

of shallow flow-system details that influence locations of 

recharge to and discharge from the deeper confined flow 

system. The increased density of groundwater associated 

with the transition from fresh to salty groundwater near the 

Atlantic Ocean affects regional groundwater flow and was 

simulated with the Variable Density Flow Process of SEAWAT 

(a U.S. Geological Survey program for simulation of threedimensional variable-density groundwater flow and transport). 

The groundwater density distribution was generated by a 

separate 108,000-year simulation of Pleistocene freshwater 

flushing around the Chesapeake Bay impact crater during 

transient sea-level changes. Specified-flux boundaries simulate 

increasing groundwater underflow out of the model domain 

into Maryland and minor underflow from the Piedmont 

Province into the model domain. Reported withdrawals 

accounted for approximately 75 percent of the total groundwater withdrawn from Coastal Plain aquifers during the year 

2000. Unreported self-supplied withdrawals were simulated in 

the groundwater model by specifying their probable locations, 

magnitudes, and aquifer assignments on the basis of a separate 

study of domestic-well characteristics in Virginia. 

The groundwater flow model was calibrated to 7,183 

historic water-level observations from 497 observation wells 

with the parameter-estimation codes UCODE–2005 and 

PEST. Most water-level observations were from the Potomac 

aquifer system, which permitted a more complex spatial 

distribution of simulated hydraulic conductivity within the 

Potomac aquifer than was possible for other aquifers. Zone, 

function, and pilot-point approaches were used to distribute 

assigned hydraulic properties within the aquifer system. The 

good fit (root mean square error = 3.6 feet) of simulated to 

observed water levels and reasonableness of the estimated 

parameter values indicate the model is a good representation 

of the physical groundwater flow system. The magnitudes 

and temporal and spatial distributions of residuals indicate no 

appreciable model bias.

The model is intended to be useful for predicting changes 

in regional groundwater levels in the confined aquifer system 

in response to future pumping. Because the transient release of 

water stored in low-permeability confining units is simulated, 

drawdowns resulting from simulated pumping stresses may 

change substantially through time before reaching steady state. 

Consequently, transient simulations of water levels at different 

future times will be more accurate than a steady-state simulation for evaluating probable future aquifer-system responses to 

proposed pumping. 

Introduction

During the 20th century, groundwater withdrawals 

increasingly supplied the domestic, municipal, and industrial 

water demands of the growing population of the Virginia 

Coastal Plain (fig. 1). Groundwater withdrawals from the 

Virginia Coastal Plain aquifers increased substantially after 

the onset of U.S. involvement in World War II (fig. 2), 

which resulted in regional cones of depression in which 

drawdown exceeded 200 feet (ft) in places by 2003. By 

the early 1970s, concern over the availability, quality, and 

sustainability of groundwater resources resulted in the 

Virginia Groundwater Act of 1973. During subsequent years, 

continued concern over the sustainability of fresh groundwater 

resources and the increased potential for seawater intrusion 
2 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

resulted in management of groundwater withdrawals through 

a State-controlled regulatory process. The Virginia Ground 

Water Management Act of 1992 empowered a State Water 

Control Board to create Ground Water Management Areas 

(GWMA), within which groundwater withdrawals greater 

than 300,000 gallons (gal) per month must be permitted by the 

Department of Environmental Quality (DEQ). The southeastern two-thirds of the Virginia Coastal Plain are currently a 

declared GWMA. Among the criteria for issuing groundwater 

withdrawal permits are the evaluation of the “area of impact” 

of a withdrawal and an assessment of the probable additional 

groundwater drawdown resulting from the proposed 

withdrawal. A regional groundwater flow model (Harsh and 

Laczniak, 1990) developed as part of the U.S. Geological 

Survey (USGS) Regional Aquifer System Analysis (RASA) 

Program was modified (McFarland, 1998) and used for these 

Figure 1. Locations of counties and independent cities, the Chesapeake Bay impact crater, and other 

important geographic and physiographic features of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

39 ̊

38 ̊

37 ̊

78 ̊ 77 ̊ 76 ̊

PRINCE

WILLIAM

STAFFORD

CAROLINE

SPOTSYLVANIA

KING GEORGE

Fredericksburg

FAIRFAX

KING AND QUEEN

KING WILLIAM

MIDDLESEX

LANCASTER

RICHMOND

NORTHUMBERLAND

ESSEX

WESTMORELAND

HANOVER

HENRICO NEW KENT

Richmond

MATHEWS

GLOUCESTER

ACCOMACK

NORTHAMPTON

DINWIDDIE

SUSSEX

SURRY

ISLE OF WIGHT

CHARLES CITY

CHESTERFIELD

SOUTHAMPTON

Chesapeake

Virginia

Emporia Beach

GREENSVILLE

Suffolk

Franklin

PRINCE GEORGE

Petersburg

Hopewell

Colonial Heights

Norfolk

Portsmouth

Poquoson

Hampton Newport News

YORK

Williamsburg JAMES CITY

NORTH CAROLINA

MARYLAND

VIRGINIA

VIRGINIA

Delaware

Fall Line

Piedmont

Province

Outer limit of Chesapeake

Bay impact crater

River

James

River

Eastern Virginia Ground Water

 Management Area

EXPLANATION

Washington, D.C.

ARLINGTON

Alexandria

CHARLES

ST. MARYS

NORFOLK ARCH

Point Lookout

CHESAPEAKE

BAY

ATLANTIC

OCEAN

VIRGINIA

NORTH CAROLINA

DELAWARE

MARYLAND

Washington, D.C.

PENNSYLVANIA

VIRGINIA COASTAL PLAIN

CHESAPEAKE BAY

ATLANTIC OCEAN

Potomac River

River

Rappahannock

York

Eastern Shore

of Virginia
Introduction 3

evaluations. Water levels simulated with this modified RASA 

model differed from corresponding observed water levels 

by up to 100 ft in some areas of the Coastal Plain (Quinlan, 

2004, 2005; Wright, 2006). This magnitude of inaccuracy 

has decreased the utility of the modified RASA model as a 

water-resource management tool.

The discovery of a buried impact crater beneath the 

Chesapeake Bay (Poag and others, 1994; Poag, 1999) revolutionized the conceptualization of the Virginia Coastal Plain 

hydrogeologic framework (Powars and Bruce, 1999; Powars, 

2000). The low permeability of the tsunami-generated Exmore 

breccia within the Chesapeake Bay impact crater influences 

regional groundwater flow and the distribution of salty 

groundwater near the Atlantic Coast. Additional advancements 

in understanding the Coastal Plain hydrogeologic framework 

have been made in areas of Virginia outside of the crater.

The changes in conceptualization of the hydrogeologic 

framework and advances in hydrologic-modeling techniques 

motivated the development of the new groundwater flow 

model documented in this report. The DEQ and the Hampton 

Roads Planning District Commission (HRPDC) intend to use 

this model in place of the modified RASA model for future 

predictions of aquifer-system response to proposed groundwater withdrawals. The Virginia Coastal Plain groundwater 

regulatory scheme currently allocates withdrawals from the 

Potomac aquifer system as pumpage from either a lower, 

middle, or upper Potomac aquifer (or some combination) 

according to the RASA-era hydrogeologic framework. In 

the latest Virginia Coastal Plain hydrogeologic framework 

(McFarland and Bruce, 2006), the Potomac aquifer system is 

considered to be a single aquifer. This new conceptualization 

of the Potomac aquifer, which is incorporated in this model, 

may require a re-evaluation of the existing drawdown criteria 

used for the process of issuing withdrawal permits.

Purpose and Scope

This report describes a new computer simulation of 

groundwater flow in the aquifer system of the Atlantic Coastal 

Plain in Virginia. The simulation incorporates the recently 

revised hydrogeologic framework (McFarland and Bruce, 

2006), which includes a redefined Potomac aquifer and 

representation of the units associated with the Chesapeake 

Bay impact crater. Development of the model included the 

identification of pertinent physical properties and processes 

to be simulated, design of the numerical discretization 

scheme, specification of initial and boundary conditions, and 

calibration of simulated to observed water levels by non-linear 

regression. The primary purpose of the model is to assist in 

the evaluation of the responses of the confined aquifers to 

groundwater withdrawals. The changes to aquifer-system 

water levels since the “predevelopment conditions” prevalent 

in 1890 are simulated by a 113-year transient simulation 

consisting of 34 stress periods and 218 time steps. Watermanagement or planning authorities who wish to assess the 

possible future regional groundwater drawdown resulting from 

proposed groundwater withdrawals may modify appropriate 

model-input files to run predictive simulations. 

The scale of the model is appropriate for assessing the 

regional water-level responses in the confined aquifers of the 

Coastal Plain. Because groundwater flow in the unconfined 

Figure 2. Virginia reported and domestic groundwater withdrawals.

0

10

20

30

40

50

60

70

80

90

100

GROUNDWATER WITHDRAWALS, IN MILLION GALLONS PER DAY

Domestic

Reported

1940 1950 1960 1970 1980 1990 2000

YEAR
4 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

surficial aquifer is largely influenced by topographic variations 

that cannot be represented with the 1-mile (mi) finitedifference cells of the model, local horizontal groundwater 

flow through the surficial aquifer is not simulated to the degree 

of accuracy necessary for its assessment as a groundwater 

resource. However, simulated vertical flow through the 

surficial aquifer from and to the water table in areas of 

groundwater recharge and discharge does allow an accurate 

assessment of groundwater transmission to and from underlying confined aquifers. 

The modeled area encompasses all the Virginia Coastal 

Plain south of Stafford, VA, and the adjacent Coastal Plain 

areas of North Carolina and Maryland (fig. 3). The small 

Coastal Plain area in Virginia west of the Potomac River and 

north of Stafford, VA, is not included in the model domain. 

Although the Eastern Shore of Virginia (also known as 

the Delmarva Peninsula) is included in the model domain, 

accurate simulation of groundwater flow in the shallow 

unconfined and confined aquifers located there is not within 

the scope of this study. A more detailed groundwater flow 

model of the shallow aquifers on the Eastern Shore of Virginia 

was developed separately by Sanford and others (2009).

The regional scale of the model precludes its use for predicting the fate and transport of anthropogenic contaminants 

in surficial aquifers. The model does not accurately simulate 

horizontal flow in the surficial aquifers. However, the model 

could be used to provide boundary water levels and fluxes for 

more detailed smaller-scale models assembled to address local 

groundwater supply or contamination problems. 

The increase in groundwater density associated with the 

transition from fresh to salty groundwater near the Atlantic 

Ocean affects regional groundwater flow. An approximation 

of the distribution of modern groundwater density was 

generated with an ancillary 108,000-year transport simulation 

of the evolution of the saltwater transition zone during the 

Pleistocene epoch. The details of this simulation are described 

in the appendix.

The 113-year transient model described in this report 

does not simulate saltwater intrusion. Unlike the ancillary 

108,000-year transport simulation, the variable-density 

groundwater flow equation was not coupled to the masstransport equation for the 113-year transient simulation. 

Therefore, the three-dimensional distribution of groundwater 

density does not change with time in the 113-year simulation. 

Although the numerical dispersion associated with 1-mi-wide 

finite-difference cells is acceptable for generating the regional 

groundwater density distribution with the 108,000-year 

transport simulation, it is not adequate when the simulation 

of temporal changes in solute concentration require a more 

accurate solution of the mass-transport equation. Smaller 

finite-difference cells would be required to accurately simulate 

solute transport. Appropriate representation of aquifer 

dispersivity is required for simulation and prediction of 

solute transport, which requires that numerical dispersion be 

minimized. However, the more-refined spatial and temporal 

discretization and computationally intensive solution algorithm required to accomplish this are not compatible with the 

scale and scope of the 113-year transient model, nor are the 

long simulation-execution times required by such algorithms 

compatible with the need for the large number of simulation 

runs required during parameter estimation.

Previous Studies

An extensive review of previous geologic and hydrologic 

studies of the Virginia Coastal Plain was recently provided by 

McFarland and Bruce (2006); a short synopsis from selected 

studies follows.

Sanford (1913) surveyed the groundwater resources 

of the Virginia Coastal Plain and provided a preliminary 

overview of the hydrogeology of the area. During and after 

the acceleration of groundwater development coinciding with 

the start of U.S. involvement in World War II, Cederstrom 

(1943, 1945, 1957, 1969) reported on more detailed studies of 

Coastal Plain hydrogeology. His groundwater quality analyses 

recognized the existence of an anomalous “inland saltwater 

wedge,” now known to coincide with the location of the 

Chesapeake Bay impact crater.

As a part of the USGS RASA Program, Meng and Harsh 

(1988) conducted the first comprehensive hydrogeologic 

analysis of the entire Virginia Coastal Plain aquifer system, 

in which they defined a hydrogeologic framework consisting 

of nine aquifers with eight intervening confining units. In 

this framework, the Potomac aquifer system was subdivided 

into the upper, middle, and lower Potomac aquifers, which 

were separated by laterally extensive confining units. Harsh 

and Laczniak (1990) used this framework to construct the 

first digital model of Coastal Plain groundwater flow in 

Virginia and adjacent parts of Maryland and North Carolina. 

This RASA model was subsequently modified (McFarland, 

1998) for use as a management tool by personnel involved 

in Virginia’s regulatory and groundwater withdrawal permit 

process. Leahy and Martin (1993) noted the differences in 

simulated water levels between constant- and variable-density 

formulations of groundwater flow under the northern Atlantic 

Coastal Plain, including Virginia. Richardson (1994a) 

simulated a sharp transition from freshwater to saltwater with 

a dual-density model of groundwater flow in the shallow 

aquifers on the Eastern Shore of Virginia.

The discovery of the Chesapeake Bay impact crater 

enabled a major advancement in understanding the stratigraphy of the eastern part of the Coastal Plain, which was 

documented by Powars and Bruce (1999) and Powars (2000). 

McFarland and Bruce (2006) have provided a thorough 

revision to the earlier RASA hydrogeologic framework, 

including incorporation of hydrogeologic units associated with 

the Chesapeake Bay impact crater. 
Introduction 5

Figure 3. Finite-difference model grid of the Virginia Coastal Plain with row and column numbers of cross sections.

Delaware

Maryland

North Carolina

Virginia

Cross section

0

0

25 50 MILES

25 50 KILOMETERS

80

95

107

120

50 60 70 87

A A’

EXPLANATION

E’

E F G H

F’ G’ H’

A

B

C

D

A’

B’

C’

D’

38°

37°

36°

77° 76° 75°
6 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Description of the Virginia 

Coastal Plain

The Coastal Plain of Virginia is underlain by Early 

Cretaceous and younger unconsolidated sedimentary deposits 

that overlie older consolidated sedimentary, igneous, and 

metamorphic bedrock. The “Fall Line” separates the Coastal 

Plain from the Piedmont Physiographic Province to the west 

(fig. 1), which is composed of consolidated lower Mesozoic 

rift-basin and older crystalline bedrock that also underlie the 

Coastal Plain sedimentary strata. The “Fall Line” is named 

for the occurrence of rapids and small waterfalls in rivers that 

cross the contact from the consolidated rocks of the Piedmont 

to the more easily eroded Coastal Plain sediments.

Most of the topography of the Virginia Coastal Plain is 

flat to gently rolling and interspersed with numerous small 

gaining streams. Some deeply incised stream valleys accentuate the rolling topography in the northwest part of the Coastal 

Plain. The land-surface altitude ranges from zero to about 

350 ft in the northwest part of the Coastal Plain, near the Fall 

Line in Stafford County. Three large rivers—the Rappahannock, York, and James—separate peninsular regions known as 

the Northern Neck, Middle, and York-James Peninsulas. The 

Northern Neck is bounded on the north by the Potomac River. 

The area south of the James River is known as southeastern 

Virginia. The Chesapeake Bay separates the mainland Coastal 

Plain from the Delmarva Peninsula, the Virginia portion of 

which is also known as the Eastern Shore (fig. 1). 

The Coastal Plain climate is temperate; summers are 

hot and humid, and winters are moderate. The mean annual 

temperature is about 58 degrees Fahrenheit (°F). Precipitation 

throughout the year is predominantly rainfall, with occasional 

light snowfall in winter. The mean annual precipitation for the 

Coastal Plain is 45.9 inches (in.; Virginia State Climatology 

Office, 2007).

Hydrogeologic Framework

The sediments that underlie the Virginia Coastal Plain 

reflect changing depositional conditions during the past 

120 million years, resulting in a sedimentary assemblage of 

gravels, sands, silts, and clays with a complex distribution 

of hydraulic properties. The Norfolk Arch (fig. 1) is an 

east-west trending structural high that separated post-Jurassic 

sedimentation into northern and southern depocenters, termed 

the Salisbury and Albemarle embayments, respectively. 

Cretaceous fluvial-deltaic and early Tertiary marine deposition 

formed a generally eastward-thickening unconsolidated sedimentary assemblage that attains a thickness of over 4,000 ft 

near the Atlantic coast. A mid-Eocene bolide impact excavated 

the post-Jurassic sedimentary strata and some older underlying 

rocks within a 56-mi-diameter area around Cape Charles, on 

the present-day Eastern Shore of Virginia. At the time of the 

impact, the sea level was approximately 200 ft higher than the 

present day, and the shore line was approximately coincident 

with the Fall Line. Tsunamis generated by this bolide impact 

on the submerged continental shelf re-deposited a chaotic mix 

of the excavated strata as clast- and matrix-supported breccias within the excavated crater (Powars and Bruce, 1999). 

Subsequent late-Eocene through Pliocene low-energy marine 

deposition buried the impact structure and the surrounding 

undisrupted Coastal Plain sediments with more clay, silt, and 

sand. The lower half of the crater contains rotated blocks and 

submarine-landslide deposits resulting from the impact event 

but not associated with the tsunami-resurge deposits.

The Coastal Plain strata have been classified into a 

layered sequence of aquifers and confining units, each of 

which generally thicken and deepen to the east toward the 

Atlantic Ocean. McFarland and Bruce (2006) describe the 

structural configuration of 19 hydrogeologic units within the 

Virginia Coastal Plain. Borehole electric logs were used along 

with drill-cutting and undisturbed core samples to define the 

altitude of the top of each hydrogeologic unit at borehole 

control points. Cross-borehole correlation enabled generation 

of two-dimensional surfaces that define the top altitude 

over the horizontal extent of each unit. Subtraction of the 

altitudes of adjacent surfaces generated raster representations 

of the thickness of each hydrogeologic unit. The composite 

assemblage of the top surfaces and associated thicknesses 

of all 19 hydrogeologic units comprise a three-dimensional 

model of the Coastal Plain hydrogeologic framework. This 

framework provides the primary basis for the distribution 

of model parameters representing the hydraulic properties 

that control the simulation of groundwater flow. None of 

the 19 hydrogeologic units are present throughout the entire 

domain of the groundwater model; they pinch out at lateral 

depositional boundaries, some are truncated by the Chesapeake Bay impact crater, and shallow units are discontinuous 

across incised river drainages. The complex three-dimensional 

distribution of these hydrogeologic units can be appreciated 

with computer software for viewing the solid model of the 

hydrogeologic framework. A series of east-west and northsouth two-dimensional cross sections (the locations of which 

are shown in figure 3) through this solid model that illustrate 

the vertical relations between the various hydrogeologic units 

are shown in figure 4. The horizontal extents of each aquifer 

and confining unit or zone are presented in figures later in 

the report.

McFarland and Bruce (2006) classified Cretaceous and 

younger unconsolidated Virginia Coastal Plain sediments into 

19 hydrogeologic units, of which 8 are identified as aquifers, 

7 are confining units, and 4 are considered “confining zones.” 

The age, stratigraphic sequence, and relation to geologic 

formation names in Virginia are summarized in figure 5. 

The reader is referred to McFarland and Bruce (2006) for an 

explanation of hydrologic associations that cross stratigraphic 

boundaries. McFarland and Bruce (2006) distinguish 

“homogenous aquifers,” which “can have locally gradational 

changes in texture and(or) composition” from “heterogeneous 

aquifers,” which “exhibit sharp contrasts in sediment texture 
Description of the Virginia Coastal Plain 

7

Figure 4. Vertical model cross sections showing distribution of hydrogeologic units from west to east along row 80 (A–A ́), west to east along row 95 (B–B ́), west to east along row 107 (C–C ́), west to east along row 120 (D–D ́), north to south along column 50 (E–E ́), north to south along column 60 (F–F ́), north to south along column 70 (G–G ́), and north to south along column 87 (H–H ́). (Locations of cross sections are shown in figure 3.)

West

Basement Basement

Basement Basement

East West East

A A’ B B’

West West East

C D C’ D’

East

–1,670

–2,470

210

–870

0

–1,670

–2,470

210

–870

0

ALTITUDE, IN FEET ABOVE OR BELOW NGVD 29

EXPLANATION

Surficial aquiferYorktown-Eastover confining zoneYorktown-Eastover aquiferSt. Marys confining unitSt. Marys aquifer

Calvert confining unitPiney Point aquiferChickahominy confining unitExmore matrix confining unitExmore clast aquifer

Nanjemoy-Marlboro confining unitAquia aquiferPeedee confining zonePeedee aquiferVirginia Beach confining zone

Virginia Beach aquiferUpper Cenomanian confining unitPotomac confining zonePotomac aquifer

0 40 MILES 10 20 300 40 KILOMETERS 10 20 30

VERITCAL SCALE GREATLY

EXAGGERATED
8 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 4. Vertical model cross sections showing distribution of hydrogeologic units from west to east along row 80 (A–A ́), west to east along row 95 (B–B ́), west to east along row 107 (C–C ́), west to east along row 120 (D–D ́), north to south along column 50 (E–E ́), north to south along column 60 (F–F ́), north to south along column 70 (G–G ́), and north to south along column 87 (H–H ́). (Locations of cross sections are shown in figure 3.)—Continued

–1,670

–2,470

210

–870

0

–1,670

–2,470

210

–870

0

ALTITUDE, IN FEET ABOVE OR BELOW NGVD 29

North South

E E’

North South

F F’

North South

G G’

North South

H H’

Basement

Basement

Basement

EXPLANATION

Surficial aquiferYorktown-Eastover confining zoneYorktown-Eastover aquiferSt. Marys confining unitSt. Marys aquifer

Calvert confining unitPiney Point aquiferChickahominy confining unitExmore matrix confining unitExmore clast aquifer

Nanjemoy-Marlboro confining unitAquia aquiferPeedee confining zonePeedee aquiferVirginia Beach confining zone

Virginia Beach aquiferUpper Cenomanian confining unitPotomac confining zonePotomac aquifer

0 100 MILES 20 40 60 80

0 100 KILOMETERS 20 40 60 80

VERITCAL SCALE GREATLY

EXAGGERATED
Description of the Virginia Coastal Plain 9

Figure 5. Stratigraphic correlations of hydrogeologic units of the Virginia Coastal Plain (modified 

from McFarland and Bruce, 2006.) Vertical arrows indicate major hydrologic associations that cross 

stratigraphic boundaries; minor overlaps of hydrogeologic units among adjacent geologic formations 

are not depicted.

GEOLOGIC UNIT

Powars and Bruce, 1999;

Powars, 2000

HYDROGEOLOGIC UNIT

McFarland and Bruce, 2006

ERATHEM/

ERA PERIOD SERIES/

EPOCH

Quaternary

Holocene undifferentiated

Pleistocene surficial aquifer

Bacons Castle Yorktown

confining zone late Chowan River Pliocene

Yorktown Yorktown confining zone

early Yorktown-Eastover

aquifer Eastover

late St. Marys

confining unit 

St. Marys Miocene St. Marys aquifer

middle Calvert confining unit

Calvert

Cenozoic

Tertiary

early Piney Point

aquifer Oligocene late Old Church

early Delmarva Beds

Chickahominy Chickahominy confining unit

late

Exmore

tsunami-breccia

Exmore matrix confining unit

Exmore clast confining unit

Eocene megablock beds Potomac confining zone

Potomac aquifer

middle Piney Point Piney Point

aquifer Nanjemoy early

Marlboro

Clay

Nanjemoy-Marlboro

confining unit

late Aquia Paleocene Aquia aquifer

early Brightseat

Peedee

clayey silty sand confining zone

organic-rich clay

glauconite quartz sand Peedee

Late red beds aquifer

Virginia Beach confining zone

glauconitic sand Virginia Beach aquifer

upper Cenomanian

beds

Upper Cenomanian

Me confining unit sozoic

Cretaceous

Potomac confining zone

Early Potomac Potomac aquifer

Jurassic

Triassic Undifferentiated Basement Paleozoic

Proterozoic
10 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

across small distances in the form of discontinuous and locally 

variable fine-grained sediments that are interbedded with 

coarse-grained sediments.” The surficial, Yorktown-Eastover, 

Peedee, and Potomac aquifers are considered heterogeneous, 

and the St. Marys, Piney Point, Aquia, and Virginia Beach 

aquifers are considered to be homogeneous.

McFarland and Bruce (2006) also identify “...indistinguishable relations that can exist adjacent to heterogeneous 

aquifers. Confining zones bound major parts of heterogeneous 

aquifers that extend across distances as great as several tens 

of miles and that are separated from adjacent aquifers solely 

by interbeds because no confining unit is present. Moreover, 

neither an interbed nor a confining unit may be present at 

many locations across the confining zones, thereby allowing 

direct contact of the coarse-grained sediments of the adjacent 

aquifers.” Note that areas that are not considered “confining” 

by common hydrologic meaning may be included in units designated as a “confining zone” according to this stratigraphic 

usage of this term by McFarland and Bruce (2006).

In contrast to the earlier three-aquifer conceptualization 

by Meng and Harsh (1988), the Potomac aquifer system is 

currently considered to be a three-dimensional fluvial-deltaic 

network of anastomosing sandy river channels and finegrained overbank deposits composing a single hydrogeologic 

unit. McFarland and Bruce (2006) determined that the lowresistivity signals identifying individual fine-grained layers 

within the Potomac Formation generally cannot be correlated 

between boreholes separated by more than several thousand 

feet. If more laterally extensive fine-grained layers exist 

within the Potomac Formation, they have not been mapped 

or documented as of 2009. The maximum lateral extent of 

most of the numerous low-permeability clay interbeds within 

the Potomac Formation is, therefore, probably on the order of 

several hundred feet.

Hydrologic Conditions

Climatic conditions and surface-water features influence 

the recharge to and flow through the aquifer system. Before 

the onset of extensive groundwater extraction in the 20th century, shallow groundwater flowed toward discharge areas 

at rivers, the Chesapeake Bay, and the Atlantic Ocean and 

also discharged by evapotranspiration in low-altitude areas. 

Deeper recharged groundwater flowed east toward a mixing 

zone with higher-density seawater, where it was deflected 

upward toward discharge areas near the coast above the deeper 

higher-density salty groundwater. These generalized patterns 

of predevelopment groundwater flow are depicted in figure 6. 

The spatial distributions of the hydraulic properties associated 

with the various hydrogeologic units, as well as the recharge 

fluxes and their geometric relation to discharge areas at the 

land or sea-floor surface, all controlled the directions and 

Figure 6. Generalized predevelopment groundwater flow patterns in the Coastal Plain aquifer system 

(modified from Leahy and Martin, 1993).

Confined

aquifer

Unconfined

aquifer

Freshwater

Transpiration

Water table

Recharge

Evaporation

Stream

Major bay

Atlantic

Ocean

Fall

line

Sea level

Confining unit

Confining unit

Confined

aquifer

Not to scale

Freshwater - Saltwater

interface

Basement bedrock

Precipitation

Saltwater

Groundwater flow paths–

 shows general direction of

 groundwater flow

EXPLANATION
Description of the Virginia Coastal Plain 11

rates of predevelopment groundwater flow. After the installation of pumping wells, the commencement of groundwater 

withdrawals introduced stresses within the aquifer system that 

modified the groundwater flow rates and directions from their 

predevelopment conditions.

Precipitation

For the period from 1895 through 2001, the mean annual 

precipitation for Virginia was 42.5 in. (National Climate 

Data Center, 2002). The Coastal Plain receives slightly more 

precipitation than the remainder of Virginia; the average of 

annual precipitations recorded at Richmond and Norfolk, VA, 

for the period 1895–1998 was 44.6 in. (Virginia State Climatology Office, 2007). July and August are the wettest months 

for the Coastal Plain; monthly precipitation for those months 

averages slightly over 5 in., and other months average between 

3 and 4 in. of precipitation. 

Precipitation is a large component of groundwater 

recharge, which is specified as a temporally variable boundary 

flux in the groundwater model. A simple analysis of the spatial 

and temporal variations of Coastal Plain precipitation was 

necessary to quantify groundwater recharge for the model. The 

locations of precipitation-measurement stations at Williamsburg, Norfolk, Hopewell, and Fredericksburg (Virginia 

State Climatology Office, 2007) are shown in figure 7, and 

the annual totals of monthly precipitation measured at these 

stations over the period from 1891 through 2003 are depicted 

in figure 8. The mean annual precipitation for these measurement sites for the period from 1895 through 2003 is 45.9 in. 

The precipitation records illustrate good correlation of total 

precipitation between these four recording stations. A change 

of mean annual precipitation measured at Williamsburg is 

apparent in figure 8. Although precipitation measured at Williamsburg is consistently higher than precipitation measured at 

Fredericksburg, the difference between these stations becomes 

less pronounced during the period after 1978. The mean 

annual precipitation at Williamsburg from 1905 through 1953 

was 59.1 in., and from 1958 through 2003 the mean annual 

precipitation was 48.2 in.

Evapotranspiration

Water evaporates from the surfaces of land and water 

bodies, the water table, the intervening unsaturated zone, 

and through the plant stomata (transpiration). It is difficult to 

separate the amounts of transpiration and evaporation over 

vegetated land areas, and these two processes commonly are 

combined and represented by an evapotranspiration (ET) term 

in hydrologic models. Annual potential ET (PET) is the maximum ET flux that could occur without a shortage of water. The 

Virginia State Climatology Office (2007) has calculated PET 

for sites in Virginia by the Thornthwaite (1948) method, which 

is based on an empirical relation between PET and mean 

air temperature. The mean annual PET for five sites in the 

Coastal Plain (Richmond, Williamsburg, Suffolk, Norfolk, and 

Fredericksburg) is 32 in.

ET from the water table is considered separately from 

vadose-zone ET because the groundwater model documented 

in this report only simulates saturated-zone processes. In 

the absence of plants, evaporative flux decreases with depth 

to the water table. Gardner and Fireman (1958) studied the 

evaporation rate as a function of water-table depth for two 

soil types (sandy loam and clay); evaporation rate decreased 

to less than 10 percent of the maximum rate for both soil 

types at depths below about 8 ft. Transpiration losses can 

occur to the maximum depth to which plant roots are able to 

penetrate. Canadell and others (1996) reported rooting depths 

by vegetation and climate types. In temperate climates similar 

to that of the Virginia Coastal Plain, they reported average 

maximum-root depths of 9.5 ft (2.9 meters [m]) for deciduous 

trees and 12.8 ft (3.9 m) for coniferous trees. In the Virginia 

Coastal Plain, ET losses from the saturated groundwater 

system are probably negligible when the water table is deeper 

than about 13 ft below land surface.

Recharge

The fraction of precipitation infiltrating Coastal Plain 

sediments that is not evaporated or transpired by plants may 

recharge aquifers and subsequently flow toward areas of 

groundwater discharge. In areas not covered by a surfacewater feature, such as a lake or stream, recharge to the 

saturated groundwater flow system primarily occurs from deep 

percolation of precipitation through the unsaturated zone to 

the water table. Some precipitation that falls on the ground 

surface (P) runs off to surface streams and rivers without 

infiltrating the ground (Qrunoff). Some of the water that does 

infiltrate below ground surface is evaporated or transpired 

from the vadose zone (ETvadose). If the moisture content within 

the unsaturated zone is sufficiently high, the remaining 

infiltrated water percolates to the water table and is considered 

recharge to the saturated groundwater system. Secondary (and 

volumetrically minor compared to infiltrated precipitation) 

contributions to recharge are effluent from septic systems, 

return percolation from irrigation, and leakage from water and 

septic pipelines (Qother). These components of recharge (R) are 

summarized by the mass-balance equation:

R= P + Qother – Qrunoff – ETvadose.

The quantity R was required for the simulation and was 

estimated as a fraction of P; quantification of the components 

in the remaining fraction (Qother – Qrunoff – ETvadose) was not 

necessary. Groundwater at or below the water table may also 

be lost to evaporation or transpiration by plants or trees with 

roots that penetrate below the depth of the water table (ETsat). 

The remaining net recharge to the saturated groundwater flow 

system is:

 R – ETsat = P + Qother – Qrunoff – ETvadose – ETsat. 

(1)

(2)
12 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 7. Locations of precipitation-measurement stations, control points, and rivers simulated with the River 

package, other streams, and sea-level areas simulated with the General Head Boundary Package in the Virginia 

Coastal Plain. (River-name abbreviations reference table 5.)

0

0

25 50 MILES

25 50 KILOMETERS

DelMarva Peninsula

rapp

poni

Fredericksburg

marr

catp

pisc

pian

wico

jame

chic

toto

soua

nora

appo

ston assa

Hopewell

swif

blac rowa

Williamsburg

MARYLAND

VIRGINIA

NORTH CAROLINA

seac

rivb

some

nans

thre

mehe

font

kirb

pote

icwat

norw

Norfolk

Figure 7. Location of precipitation-measurement stations, control points, and rivers simulated with the River

Package, other streams, and sea-level areas simulated with the General Head Boundary Package in the

Virginia Coastal Plain. (River-name abbreviations reference table 6.)

mill

CHESAPEAKE

BAY

ATLANTIC

OCEAN

General head boundary

Non-simulated stream

Simulated stream

Precipitation-measurement

station

River control point

EXPLANATION

37 ̊

38 ̊

77 ̊ 76 ̊
Description of the Virginia Coastal Plain 13 Figure 8. Historical measured precipitation at the four precipitation-measurement stations in the Virginia Coastal Plain. (Site locations are shown in figure 7. Data from Virginia State Climatology Office, 2007.) 2030

405060708090

1890 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000

ANNUAL PRECIPITATION, IN INCHES

Norfolk Williamsburg Hopewell Fredericksburg

YEAR
14 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Richardson (1994b) used hydrograph-separation 

(Rutledge, 1993) and hydrogeologic-area approaches to 

estimate an average discharge flux (base flow) of about 

10 in/yr from the Coastal Plain of Virginia. The discharge 

flux from that study is an estimate of the effective recharge 

quantified in equation 2.

Groundwater Underflow from the 

Piedmont Province

Groundwater flow in the Piedmont Province of Virginia 

generally is thought to be compartmentalized between 

topographic drainages, and “water is not gained from or lost 

to deeper regional ground-water flow systems...” (Powell 

and Abe, 1985). In the Piedmont, recharge over interfluves 

contributes to storage in the shallow regolith and the flow 

though the deeper crystalline bedrock fractures that discharges 

to local rivers and streams. On the basis of this conceptual 

model, previous Coastal Plain modeling studies have assumed 

that there is no continuity of groundwater flow across the Fall 

Line from the Piedmont Province into the Coastal Plain. Harsh 

and Laczniak (1990) state:

The westernmost boundary of the model coincides 

with the Fall Line and is considered impermeable 

to flow. This assumption is supported by the large 

difference in hydraulic conductivity between the 

rocks of the Piedmont physiographic province and 

the unconsolidated sediments of the Coastal Plain 

physiographic province.

A statistical analysis by Daniel (1989) of North Carolina 

Piedmont well records indicated that well yields increase 

with depth to a maximum at 700 to 800 ft below land surface, 

a greater depth than previously thought. It is possible that 

some flow through deeper Piedmont bedrock fractures may 

be regional and discharges through Coastal Plain sediments to 

the east. The magnitude of this underflow is unknown but is 

almost certainly a minor component of the Piedmont groundwater budget and a minor source of water to the Coastal Plain 

groundwater flow system. 

Tensional stress in the Norfolk Arch may have fractured 

Piedmont rocks, thereby enhancing fracture permeability and 

the flow of groundwater toward the Coastal Plain over the 

Norfolk Arch (fig. 1). The larger area of the Piedmont Province in southern Virginia may also contribute to more deep 

groundwater infiltration and underflow to Coastal Plain sediments in southern Virginia than north of the Norfolk Arch. It 

is speculative if either or both of these effects cause a greater 

flux of Piedmont groundwater to the Coastal Plain in southern 

Virginia than further north. In this modeling study, the effect 

of possible groundwater underflow from the Piedmont to the 

Coastal Plain over the Norfolk Arch was simulated.

Hydraulic Properties

In steady-state groundwater flow, the water-level 

distribution and groundwater flux is controlled by hydraulic 

conductivity (K), which depends on the intrinsic permeability 

of the porous media and the density and viscosity of the 

fluid. The variability of these fluid properties (which depend 

on temperature and solute concentration) in Coastal Plain 

groundwater is small and has a relatively insignificant effect 

on changes of hydraulic conductivity. By contrast, the difference in hydraulic conductivity between a gravel and marine 

clay, which represent end-members on the range of intrinsic 

permeability expected for the sedimentary deposits 

comprising Coastal Plain aquifers and confining units, may 

be 10 orders of magnitude (Domenico and Schwartz, 1990). 

Within the Coastal Plain aquifer system, changes in hydraulic 

conductivity associated with the different sedimentary 

materials occur with different magnitudes at various scales. 

Hydraulic conductivity contrasts within the fluvial-deltaic 

Potomac aquifer system may be greater than those in 

hydrogeologic units formed under more uniform marine 

depositional environments.

The response of the aquifer or confining unit to a change 

in hydraulic stress under transient groundwater flow condition 

is also influenced by a storage property, which is represented 

in the groundwater flow equation by either specific storage (Ss) 

or specific yield (Sy).

Hydraulic Conductivity

Hydraulic conductivities determined from laboratory 

measurements of core-size samples and analyses of aquiferperformance tests may differ greatly because of scale effects, 

sampling, and testing and analytical methodologies. Review 

of the range of hydraulic conductivities determined by these 

methods is useful for comparison to those appropriate at the 

scale of the regional groundwater flow model, which were 

estimated during model calibration. These laboratory or field 

measurements of hydraulic conductivity also help establish the 

expected “reasonable range” for model parameters representing the horizontal or vertical hydraulic conductivity of various 

hydrogeologic units.

McFarland and Bruce (2006) reported measured vertical 

hydraulic conductivities from 33 sediment core samples 

collected from 3 boreholes within the area of the Chesapeake 

Bay impact crater; these hydraulic conductivities are 

considered representative only within an order of magnitude. 

The largest vertical hydraulic conductivity of 2.8 feet per day 

(ft/d) was measured from a sand collected from the Potomac 

aquifer. Permeameter measurements of samples representing 

the Yorktown-Eastover and Piney Point aquifers ranged 

from 6x10–4 to 2x10–2 ft/d. Permeameter measurements of 

samples from Potomac sands ranged from 6x10–4 to 2.8 ft/d. 

Permeameter measurements of samples from Potomac clay 

interbeds ranged from 3x10–5 to 2x10–3 ft/d. Permeameter 

measurements of samples representing the St. Marys, Calvert, 
Description of the Virginia Coastal Plain 15

Chickahominy, and Exmore matrix and clast confining units 

ranged from 3x10–5 to 6x10–2 ft/d. Separate measurements 

representing the same confining unit differ by up to three 

orders of magnitude. Measured hydraulic conductivities are 

not correlated to sample burial depth in any of the sampled 

hydrogeologic units. McFarland and Bruce (2006) state that 

non-cohesive coarse-grained sediments could not be collected, 

and consequently, the higher hydraulic conductivities expected 

for those sediments may be underrepresented in their 

sample set. 

Transmissivity (T) values from analyses of Potomac 

aquifer-performance tests are summarized in table 1. These 

data were provided by the Virginia DEQ (Robin Patton, 

written commun., 2004). In most cases, either different 

tests of the same well, different analyses of the same data, 

or analyses using different observation wells resulted in a 

range of transmissivity values associated with each pumping 

well (table 1). The transmissivity values were divided by the 

screened interval of each well to estimate hydraulic conductivities, which range from 15 to 337 ft/d. Multiple analyses of 

some wells give a range of transmissivity (and subsequently, 

hydraulic conductivity) values. For example, hydraulic 

conductivities from three analyses of aquifer-performance data 

(by two consulting firms and the Virginia DEQ) from well 

161-00371 (located in the City of Suffolk and operated by the 

Western Tidewater Authority), range from 83 to 337 ft/d. Low 

hydraulic conductivities from each aquifer performance test 

range from 15 to 175 ft/d, with a mean of 59 and a standard 

deviation of 42 ft/d.

Specific Storage

When water levels change under confined conditions, 

the amount of water stored or released from a representative volume of the porous media is due to the volumetric 

compressibility of water and the porous medium. These are 

quantified by the elastic specific storage, which is approximated by the expression: 

Sse=rf g(a+nbw),

Table 1. Summary of Potomac aquifer-performance test results.

[-, no analysis; T, transmissivity; ft2/d, feet squared per day; S, storage coefficient; K, hydraulic conductivity; Ss, specific storage; ft/d, 

feet per day; ft–1, per foot. Data provided by the Virginia Department of Environmental Quality] 

Well ID T low

(ft2/d)

T high

(ft2/d) S low S high

Interval

length 

(ft)

K low

(ft/d)

K high

(ft/d)

Ss low

(ft–1)

Ss high

(ft–1)

161-00380 11,963 - - - 202 59 - - -

146-00304 5,665 6,564 - - 76 75 86 - -

216-00086 6,112 - 0.00022 - 35 175 - 6.3E-06 -

216-00036 3,789 6,234 - - 66 57 94 - -

161-00397 8,328 15,774 0.00082 0.0011 150 56 105 5.5E-06 7.3E-06

187-00026 2,693 3,456 0.000052 0.000298 16 168 216 3.3E-06 1.9E-05

147-00262 8,082 8,823 0.000385 - 220 37 40 1.8E-06 -

147-00260 2,874 5,080 0.00054 - 190 15 27 2.8E-06 -

147-00294 12,346 13,700 0.000692 - 305 40 45 2.3E-06 -

147-00232 3,413 3,719 0.0007 0.0012 125 27 30 5.6E-06 9.6E-06

161-00371 31,345 127,395 - - 378 83 337 - -

161-00372 14,791 23,323 - - 380 39 61 - -

146-00250 10,999 11,782 - - 220 54 54 - -

146-00251 17,902 98,160 - - 310 58 317 - -

234-00220 2,779 2,779 0.0003 - 188 15 15 1.6E-06 -

234-00221 4,186 7,183 0.000071 0.00044 79 53 91 9.0E-07 5.6E-06

234-00222 6,961 13,334 0.00019 0.00034 94 74 142 2.0E-06 3.6E-06

234-00223 4,805 8,171 0.00016 0.00033 104 46 79 1.5E-06 3.2E-06

216-00041 7,198 - - - 455 16 - - -

161-00396 12,352 14,611 0.00067 0.001 160 77 91 4.2E-06 6.3E-06

147-00287 10,994 13,700 0.000792 - 440 25 31 1.8E-06 -

(3)
16 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

where 

	 rf is the density of water, 

g is the acceleration of gravity, 

	 a is the uniaxial compressibility of the porous 

media, 

n is porosity, and 

	 bw is the isothermal compressibility of water. 

For porous media with a porosity of 0.3, the component of 

specific storage resulting from the compressibility of freshwater (rf gnbw) at 5 °C is 4.3 x10–7 ft–1. Because bw is essentially constant over the range of pressures and temperatures 

in groundwater flow systems and n varies by less than an 

order of magnitude, reasonable values of specific storage 

are constrained by the plausible range of the porous media 

compressibility (a). For alluvial sediment, the lower limit of a

is approximately the same as bw, so the lower limit of reasonable values of elastic specific storage is about 1x10–6 ft–1. 

Laboratory consolidation tests of the virgin (non-recoverable) 

compressibility of normally consolidated alluvial materials 

at effective stress levels substantially higher than the in situ 

stress are not representative of the elastic compressibility 

applicable to in situ effective stress changes caused by groundwater level declines that do not exceed the preconsolidation 

stress in over-consolidated sediments. Widely cited laboratory compressibility measurements (Domenico and Mifflin, 

1965) from such tests may be useful for estimating inelastic 

specific storage (Ssv) but are too large and not applicable for 

estimating reasonable values of elastic specific storage (Sse). 

(For a further discussion of this issue see also Riley, 1998.) 

Extensometric measurements of elastic specific storage in 

alluvial aquifer systems containing interbedded aquitards 

from a variety of field sites are consistently in the range from 

1.8x10–6 to 4x10–6 ft–1 (Ireland and others, 1984; Heywood, 

1995, 1998). Using the one-dimensional vertical compaction 

model of Helm (1974), Pope and Burbey (2004) determined 

elastic specific storage values between 1.4x10–6 and 1.8x10–6

ft–1 for Coastal Plain sediments penetrated by vertical extensometers at Franklin and Suffolk, VA. Pope and Burbey (2004) 

also estimated that the inelastic specific storage values of finegrained interbeds within the Potomac aquifer system range 

from 4.6x10–6 to 1.6x10–5 ft–1 and that inelastic specific storage 

values of post-Cretaceous fine-grained sediments range from 

3x10–5 to 4.6x10–5 ft–1 at these sites. Other studies employing 

compaction modeling have similarly found that, for individual aquitards, elastic specific storage is typically < 2x10–5

ft–1 for water-level declines resulting in elastic compression, 

and inelastic specific storage is < 7x10–4 ft–1 when water-level 

declines result in inelastic compaction (Ireland and others, 

1984; Epstein, 1987).

Some of the Potomac aquifer-performance test analyses 

summarized in table 1 reported a storage-coefficient value. 

Estimates of specific storage in table 1 were obtained by 

dividing these storage-coefficient values by the screened 

interval of each well; estimates ranged from 9 x10–7 to 

1.9x10–5 ft–1. In the absence of detailed stratigraphic logs, this 

procedure provides a useful test for reasonableness of aquifer 

storage-coefficient values. However, it must be recognized that 

screened intervals are not necessarily truly representative of 

the aggregate thickness of aquifer material that contributes to 

the initial transient response to pumping. The cited maximum 

value of 1.9x10–5 ft–1 is about an order of magnitude too 

large for typical sands and may reflect a short screen in a 

thick aquifer as well as interpretative problems arising from 

borehole-storage effects, partial penetration, aquitard leakage, 

or other circumstances. 

Specific Yield

At the water table where the aquifer is unconfined, 

changes in storage occur primarily as drainage or filling of 

pore space. The fraction of aquifer volume that drains with 

a unit decline of the water table is quantified by the specific 

yield, which is a fraction of the total porosity. Johnson (1967) 

compiled specific yields for various alluvial materials from 

clays to coarse gravels and determined an average value for 

fine sands of 0.21 and an average value for coarse sands 

of 0.27. 

Discharge

Before groundwater pumping from the Coastal Plain 

began, groundwater flowed exclusively to natural discharge 

areas, some of which continue to be discharge areas today. 

Much of the recharge discharges to the numerous small 

streams and rivers (fig. 7) in the Coastal Plain. Water that 

has infiltrated deeper in the aquifer system may flow toward 

discharge areas along larger river systems, as well as to the 

Chesapeake Bay and Atlantic Ocean. Groundwater also 

discharges by evapotranspiration from the water table in 

low-altitude areas.

Pumping has increasingly resulted in the capture of 

groundwater that otherwise would have ultimately flowed to 

the natural discharge areas under predevelopment conditions. 

By the year 2003, total groundwater withdrawals from the 

Coastal Plain study area were about 117 million gallons per 

day (Mgal/d; fig. 2). Commercial, industrial, and public-supply 

withdrawals accounted for approximately 74 percent of the 

total groundwater withdrawals during the year 2000 (fig. 9). 

The well locations for these reported withdrawals are shown 

in figure 10. The majority of the groundwater withdrawal flux 

was concentrated in well fields supplying the industrial centers 

near West Point and Franklin, VA, where the largest declines 

in groundwater levels have also occurred.

Domestic groundwater withdrawals accounted for 

approximately 25 percent of the total groundwater withdrawals during 2000 (fig. 9). Although the individual locations 

and associated withdrawal magnitudes from wells supplying 

groundwater for domestic use are unknown, the spatial 

distribution of self-supplied domestic withdrawals has been 

estimated by Pope and others (2007). The estimated horizontal 

distribution of the magnitude of self-supplied withdrawals 
Description of the Virginia Coastal Plain 17

for the year 2000 is also shown in figure 10. Groundwater 

withdrawals for municipal, industrial, and commercial uses 

are sufficiently large to require reporting and are summarized 

by hydrogeologic unit under the “reported” column in table 2. 

The estimated total domestic withdrawal from each hydrogeologic unit in the Virginia Coastal Plain, excluding the Eastern 

Shore, is also summarized in table 2.

Saltwater Transition Zone

The increase in groundwater density associated with the 

transition from fresh to salty groundwater near the Atlantic 

Ocean has an effect on regional groundwater flow. Mixing 

of freshwater and saltwater has created a “transition zone,” 

through which salinity generally increases toward the ocean 

and with depth. This transition zone vertically crosses the 

layered sequence of Coastal Plain aquifers and confining units 

of contrasting hydraulic conductivities. Preferential flow and 

flushing through the relatively conductive aquifer strata has 

resulted in local vertical salinity inversions (McFarland and 

Bruce, 2005) within the transition zone.

Cederstrom (1943) mapped the distribution of dissolved 

solids in groundwater in the Atlantic Coastal Plain of Virginia 

and noted anomalously high concentrations in an “inland 

wedge,” now known to correspond with the extent of the 

Chesapeake Bay impact structure (Powers and Bruce, 1999). 

Geochemical and other geologic and hydrologic evidence 

suggest a Pliocene or earlier seawater source of the salty 

groundwater within the impact structure (McFarland and 

Bruce, 2005). The persistence of saltwater within the crater 

through the Pleistocene, and possibly since the Eocene, 

suggests that the low hydraulic conductivity of crater-fill 

sediments inhibited advective flushing by fresh groundwater. 

The occurrence of saline groundwater within the Chesapeake 

Bay impact structure and vertical salinity inversions in other 

locations indicate that geologic controls, in addition to the 

hydrodynamic control, affect the salinity distribution near the 

Atlantic coast.

The transition from fresh groundwater of meteoric origin 

to seawater occurs across a mixing zone whose horizontal and 

vertical widths are determined, in part, by the heterogeneity 

and permeability of the porous media and the periodic tidal or 

seasonal variations or longer non-periodic variations 

of both groundwater flux and sea level. Meisler (1989) 

defined the transition zone beneath the northern 

Atlantic Coastal Plain as the zone with chloride 

concentrations between 250 and 18,000 milligrams per 

liter (mg/L). By this criterion, the horizontal width of 

the transition zone in Virginia varies from about 25 mi 

in southeast Virginia to 45 mi across the Chesapeake 

Bay impact structure.

The equilibrium location of the transition zone 

is hydrodynamically controlled by the magnitude and 

distribution of recharge flux, aquifer-system hydraulic 

conductivities, and sea level. Because climatic conditions and sea level typically are changing through 

time, the transition zone may not be in equilibrium for 

the current climate and sea level. The Coastal Plain 

aquifers were entirely inundated and likely saturated 

with saltwater most recently during the Pliocene Epoch 

2 million years ago, and numerous sea-level changes 

have occurred since that time. The magnitude of 

relative sea-level lowering during even the most recent 

glacial episode, however, is uncertain. The ICE–4G 

Figure 9. Estimated groundwater withdrawals by 

category of use in the Virginia Coastal Plain during 2000 

(derived from U.S. Geological Survey, 2002).

Table 2. Magnitudes of 2003 reported and 2000 self-supplied 

domestic groundwater withdrawals from the Virginia Coastal Plain 

by hydrogeologic unit.

[Mgal/d, million gallons per day; <, less than; <<, much less than] 

Aquifer Domestic1

(Mgal/d)

Reported2

(Mgal/d)

Total

(Mgal/d)

Percent 

of total

Surficial 5.2 0.09 5.3 5

Yorktown-Eastover 10.8 .51 11.3 10

Saint Marys .05 .0 .05 <<1

Piney Point 1.9 4.6 6.5 5

Aquia 2.8 .48 3.3 3

Peedee 0 0 0 0

Virginia Beach .09 .08 .17 <1

Potomac 7.5 83 90.5 77

Total 28.3 88 117 100

1 Estimated from 2000 census data (Pope and others, 2007). Yorktown-Eastover 

includes Yorktown confining zone. Potomac includes Potomac confining zone. 

2

 Groundwater withdrawals greater than 300,000 gallons per month must be 

permitted. Reported withdrawal magnitudes provided by Virginia Department of 

Environmental Quality, Groundwater Withdrawal Permit Program.

Irrigation, 0.58 million

gallons per day, 0.5 percent

Domestic, 28.8 million gallons

per day, 24.9 percent

Public supply, 27.5 million

gallons per day, 23.8 percent

Commercial/Industrial, 58.6 million

gallons per day, 50.8 percent
18 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 10. Location of Coastal Plain reported withdrawal wells, model domain, and self-supplied withdrawal rate per 

square mile in 2000.

Figure 10. Location of Coastal Plain reported withdrawal wells, model domain, and self-supplied withdrawal rate per

square mile in 2000.

0

0

25 50 MILES

25 50 KILOMETERS

CHESAPEAKE

BAY

ATLANTIC

OCEAN

VIRGINIA

NORTH CAROLINA

! MARYLAND !

! !

! !

!!!!!!

!!

!

!

! !

!

!

!

!

!!

!

!

!!

!

!

!

!

!

!!

!!

!

!

!!!

! !

!

!

!

!!

!

!

!!!

!

!

!

!!!

!

! !

! ! !

!

!

!

!

!

!

!

!

! !

!

! !

!!!!! ! !

!!

!!

!!

!

! ! ! !

!

!

! !

! !!!

!!

! !

!

! !

!!!!!

!!!

!!!!

!!!!

!!!

!

!!!

!

! !

!

! ! !

!!

!!

!

!

!

!!

!

!

!

!

!

!

!

!

! ! !

!

!

!!!

!

!

!

!

!!

!!!

!

!

!

!

!

!

!

!!!

!

!!

! !! !

!

!

!

!

!

!!

!

! !

!

!

!

!

!

!

!

! !

!!

!

! !

!!

!

!!

!!!

!

!

!

!

!

!

!

!

!

!

!

!

!!

!!!!

!

!

!!

!

!

!

!

!

!

!

!

! !

!

! !

!

!

!

!

!

!

!

!

!

!

!

!

!

! !

!

!!

!

!

!

! ! !

!

! !

! ! ! !! !

!

!

!

!

!

!! !

! !

! !

!

!

!

!

!

!

!

! !

! !

!

!!

! !

!!!

!!

! !

!!

!

! ! !

!

! !

!

!

!

! !

!

!

!!

!

!

!

!

!!

!

!

!

!

! !!

!

!

!

!

!

!

!

!!

! ! !

!

!!

!

!!!!

!

!!

!!!

!

!

!

!!

!

!

!

!

!

!

!

!

!

!

!!!!

!

!!

! !

!

!! !

! ! !

!!!

!

!!

!

!

!

!

!

!!!!

!!

!

!

!!!!

!!!!!!

!!!

!

!

!

!!

!!

!

!

!!!

!!

!

!

!!!!

!

!

!!!

!

!!!!!

!!

! !!

!

!

! ! ! ! ! !! !!!!! !!

! ! ! ! ! ! !

!

!

! ! !! !! ! ! ! ! ! !! ! ! ! ! ! !!! ! ! ! !!! ! !! ! !! ! ! ! !

!! ! !! !! ! ! !! ! ! ! ! ! ! ! ! ! !

!! ! ! ! !!!!! ! ! !

!

!

!!!!! !! ! !!! !! ! !! !!!!! ! ! ! ! ! !! ! ! ! ! ! !! ! ! !

! ! !! !

! !!!!!! ! !!!!

!!!!!!!!! !! ! !!!!! ! ! !!! !

!!!!! ! ! !! ! !!

!

! ! ! !!! !!!!

! ! ! ! !! ! !!!!! ! ! !

! !! ! !!

!! !

! ! ! !

!!

!! ! ! ! ! !!! ! ! !!! ! !

!! ! !!! !

! !! ! ! ! !!

! !

! !! ! !! ! !!! ! ! !!! ! ! ! !! ! ! ! !!! !!! !!!!!

! !!

!! ! !! ! ! !

!!! !!! ! !! ! !! ! ! !! ! !!!! ! !!!! ! !

!!!! !!!! ! !! !! ! !

! ! !!!!!! !! !! ! ! ! !! !! !! ! ! !! ! ! !! ! !! ! ! ! ! !!! !!! !! ! ! ! ! ! ! !!!!

!

!!

! ! !! !! !! ! ! !!! ! !! !! ! !! ! ! !! !!! !! ! ! !

!!!!! ! !! !!! ! ! ! ! !! ! ! ! !! ! !! ! ! ! ! !! ! ! !!! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !!! !!! !!!!! ! ! ! !!!! ! !! ! ! ! ! ! !!! ! ! ! ! ! ! !! !!!!!!!!! !!

!!

!!!!

!!!!!! !! ! !!

!!!!!!!!! ! ! !

!! ! ! !!!! !!!!!

!!!!!!!!!!!!! ! !!!! ! ! ! ! ! ! ! !! ! !!!! ! ! ! !! ! ! ! !!!!!!! !!!!

!!!!!!!! !!!! !!!!!!!! ! ! ! !!!! !!! ! ! ! ! ! ! !!!!!!!!!!!!!!!!!

!!!!!!!!!! ! !!!! !!!!!! !! ! ! ! ! ! !! ! ! ! !!!!!!!!! !! ! ! ! ! ! !!!! !!! !!!! !! ! !! !! !! !!

! ! !! ! ! ! ! ! !!! !!!!! !!! ! ! !!! !! !! !!!! ! ! !!!!! ! !! ! ! !! ! !

!! !! !!!! ! !!! ! !! !!! ! ! !! !! !!!!!!! !!! !!! ! !! !!!!!! !! ! ! !!! !!!!! ! ! !!! !!! ! !!! !!!!!! !! ! !! ! !! !! !

!!!!!!!!!!!!!!! ! ! !!!!! !!!! !! !!! !!!!!!!!!!!!!!!!!!!!!!

!!

!!!!!!! !!! !!! !! !! ! !!!!!!!

!!!!!!! !!!!!!!!!!!!!!!! ! ! !! !! ! !!!!!!! !!!!!!!!!!!!! !

!!!

!!!! ! !!! ! !!! ! !!!!

!!!!!!!!!!!! !!!!!!! ! ! ! !!! ! !!!!!!!!!!!!

!!!!! !!!!!!!!!! ! !! !!!!!!!!! !!!!! ! ! !!! ! ! !

!

! !!!!!!!!!! ! ! !! !!!!!

! !

! ! ! ! ! ! ! ! ! ! !

! ! !

! ! ! ! ! ! ! ! !

West Point

0 to 200

201 to 400

401 to 800

801 to 1,600

1,601 to 8,948

Model domain

Reported withdrawal well

SELF-SUPPLIED WITHDRAWAL FLUX, IN

CUBIC FEET PER DAY PER SQUARE MILE

!

EXPLANATION

Franklin

38°

37°

36°

77° 76° 75°
Simulation of Transient Groundwater Flow 19

model (Peltier, 1994), which incorporates mantle viscosity and 

lithospheric flexure, computes sea level during the last glacial 

maximum (21,000 years before present) 470 ft lower than 

present sea level for the Virginia Coastal Plain. Age dating 

of basal peat below salt marshes and estuarine sediments 

in the Chesapeake Bay indicates sea level was 25 ft lower 

6,000 years ago (Larsen and Clark, 2003), which is about 

one-third of the ICE–4G model value for that time. During 

glacial episodes when the sea receded to the continental-shelfslope break, the average aquifer-system hydraulic gradient 

was much greater than present under the current sea-level 

conditions, especially east of the coastline, causing enhanced 

flushing of saltwater from aquifers and a seaward migration 

of the freshwater-saltwater transition zone. Repeated sea-level 

fluctuations during the late Tertiary and Quaternary have 

mixed freshwater with seawater, widening the freshwatersaltwater transition zone (Meisler and others, 1984). The 

general trend of sea-level rise over the past 18,000 years 

suggests that hydraulic forces affecting the location of the 

freshwater-saltwater interface may be out of equilibrium, 

causing landward migration of the transition zone.

The spatial configuration of the transition zone may 

be represented, either conceptually or from simulated or 

measured concentrations, by one or more surfaces of equal 

or iso-concentration (iso-conc). An iso-conc representing 

50 percent of the seawater concentration of either chloride or 

total dissolved solids (TDS) is convenient for approximate 

depiction of the position of a seawater transition zone with 

respect to a coastline or other reference. Such a representative 

iso-conc position can be calculated for simple groundwater 

systems under static conditions, or simulated for more 

complex natural systems under conditions of hydrodynamic 

equilibrium between the groundwater flow and the denser 

seawater. Because climatic and sea-level conditions change 

through time, the actual position of the seawater transition 

zone may differ from the equilibrium condition simulated for 

the current sea level and climate. The configuration of the 

simulated transition zone is shown in the appendix (fig. A2).

Simulation of Transient 

Groundwater Flow

Development of the groundwater model required assembly and organization of a considerable quantity of spatially 

referenced hydrologic, geologic, and topographic data. A 

geographic information system (GIS) was used to construct 

coverages of these diverse data sets, which facilitated 

preliminary data analysis, construction of model data sets, 

and visualization and analysis of model output and results. 

The altitudes of the land surface, tops of the 19 hydrogeologic 

units, and the pre-Cretaceous basement rocks are represented 

by raster data structures. Vector GIS coverages describe 

the location and attributes for water-level and water-quality 

observation wells, precipitation observations, pumping wells, 

county and independent city boundaries, state boundaries, 

surface-water bodies, geologic faults, the impact crater, etc. 

GIS coverages of the model finite-difference cells and nodes 

created with the program MODELGRID (Kernodle and Philip, 

1988) contain attributes for derived quantities required for 

model-input definition, such as domestic-withdrawal fluxes 

and hydrogeologic-unit top altitudes and thicknesses. Specification of MODFLOW (McDonald and Harbaugh, 1988; 

Harbaugh and others, 2000) River Package input from GIS 

coverages of rivers and attribute control points was facilitated 

by use of the program RIVGRID (Leake and Claar, 1999). 

Construction of coverages of model-simulated water levels 

was facilitated with the program MODTOOLS (Orzol, 1997).

Projection System

The GIS coverages utilized a Lambert Conformal Conic 

Projection system, which was chosen to minimize area and 

angular distortions in the Virginia Coastal Plain. The specific 

projection parameters include a central longitudinal meridian at –77°W and standard parallels of latitude at 37° and 

38°N. Coordinates in this projection space have units of feet 

measured positive east of the central meridian and north of 

latitude 36°N. 

Numerical Method

Spatial variations of dissolved-solid concentrations 

and consequent groundwater density can affect the patterns 

and rates of groundwater flow. For many groundwater flow 

systems, spatial variations of the density and temperature of 

the groundwater are sufficiently small that the fluid density 

may be considered constant, in which case the appropriate 

form of the groundwater flow equation may be conveniently 

solved with a simulation code such as MODFLOW. In coastal 

aquifers, the 2.5-percent increase in density across a transition zone from fresh groundwater (62.43 pounds per cubic 

foot [lb/ft3]) to seawater (63.99 lb/ft3) substantially affects 

groundwater flow. In these cases, more accurate simulations may be obtained with a simulation code that solves 

a variable-density form of the groundwater flow equation. 

If groundwater concentrations change with time, it is most 

appropriate to employ a simulation code that solves a system 

of coupled equations describing variable-density groundwater 

flow and mass transport. The major components of seawater 

are mass-conservative solutes that are not affected by chemical reactions. For the purposes of this model, the effects of 

groundwater temperature gradients on groundwater flow in the 

Virginia Coastal Plain are negligible. Accordingly, the governing equations for variable-density flow and mass transport can 

be simplified to the forms listed in equations 4 and 5:

( )

t

C

t C

P

p S s q q ∂

∂

∂

∂

+

∂

∂ − ∇ ⋅ + = 

     (4)
20 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

where 

∇ is the gradient operator, 

r is fluid density [ML–3], 

q

 is specific discharge [LT–1], 

P is fluid pressure [ML–1T–2], 

S

p is specific storage in terms of pressure 

[M–1LT2], 

t is time [T], 

θ is effective porosity, 

C is solute concentration [ML–3], 

D is the dispersion coefficient [L2T–1], 

qs is flow rate per unit volume representing 

sources or sinks [T–1], 

r is source or sink fluid density [ML–3], and

Cs is source or sink solute concentration [ML–3]. 

The dependence of fluid density on concentration can be 

approximated by a linear equation of state:

where f  is the density of freshwater (about 62.43 lb/ft3). The 

groundwater flow (eq. 4) and transport (eq. 5) equations are 

coupled by the specific discharge ( q

 ) and concentration (C) 

terms. The specific discharge ( q

 ) is expressed by the pressure 

form of Darcy’s law:

where

 is the tensor of permeability [L2], 

μ is the dynamic viscosity [ML–1T–1], 

g is the acceleration of gravity [LT–2], and 

ezˆ is a unit vector pointing downward. 

Guo and Langevin (2002) reformulated the flow equation 

(eq. 4) using Darcy’s law (eq. 7) to express the specific 

discharge ( q

 ) in terms of equivalent freshwater head. The 

resulting flow and transport equations can be iteratively solved 

by the Variable-Density Flow (VDF) and Integrated Mass 

Transport (IMT) processes, respectively, of the simulation 

code SEAWAT–2000 (Langevin and others, 2003). For some 

groundwater flow simulations, the temporal change of solute 

concentration may be relatively small and may not appreciably 

affect the regional spatial distribution of groundwater density. 

If the objective is the simulation of groundwater levels and not 

concentrations, as it is in this study, the spatial groundwater 

density distribution may be considered non-transient in these 

cases, and the governing equations for variable-density flow 

and mass transport are de-coupled. The primary objective of 

this study is to simulate groundwater levels; therefore, the 

governing equations were de-coupled. Because the iterative 

solution of the coupled equations is computationally intensive, 

solution of the variable-density flow equation for a temporally 

constant-density case using only the VDF process in 

SEAWAT–2000 saves considerable computation time.

The SEAWAT code used to simulate the variable-density 

groundwater flow is documented in both Guo and Langevin 

(2002) and Langevin and others (2003). For this study, the 

SEAWAT–2000 (Langevin and others, 2003) version of the 

code was used and is referred to as SEAWAT throughout this 

report. The specified groundwater density distribution was 

generated by a separate 108,000-year simulation of Pleistocene freshwater flushing around the Chesapeake Bay impact 

crater during transient sea-level changes and is described in 

the appendix. 

Spatial Discretization

The finite-difference groundwater model grid was 

designed to provide good spatial resolution on a regional scale 

while maintaining acceptable limits on computer execution 

times. The 24,486-mi2 area (159 mi by 154 mi) encompassed 

by the model grid extends in an east-west direction from the 

Fall Line to the continental shelf and from 31 mi south of the 

Virginia-North Carolina State line to about 25 mi north of 

Point Lookout in Maryland (fig. 1). The grid extent is almost 

twice the area of the Coastal Plain in Virginia (13,000 mi2), 

although areas on the continental shelf and in North Carolina 

have wide column and row spacing, respectively, and therefore 

did not substantially increase execution times.

The finite-difference model grid is composed of 

134 rows, 96 columns, and 60 layers; row numbers increase 

from north to south, column numbers increase from west 

to east, and layer numbers increase from the top down. 

The north-west corner of the grid (row 1 and column 1) is 

876,480 ft north and 153,120 ft west of the previously defined 

original latitude and Lambert central meridian, respectively. 

The Coastal Plain of Virginia (exclusive of the Eastern Shore) 

and part of Maryland is discretized with 1-mi2

 grid-cell 

spacing. The grid spacing along rows (north-south direction) 

within the southernmost six rows of finite-difference cells in 

the model, which are south of the Virginia-North Carolina 

State line, increases through 2, 3, and 5 mi to a maximum 

of 7 mi. East of the mainland Virginia Beach area, the grid 

spacing along columns (east-west direction) increases through 

2, 3, 5, and 7 mi to a maximum of 10 mi.

Median land-surface altitudes and bathymetric depths 

for each model cell (1 mi2 or more) derived from 30-m digital 

elevation models (DEMs) were used to define the uppermost 

active layer in a particular row and column. The vertical thickness of the upper 48 layers is uniformly 35 ft. Layers 49–52 

are 50 ft thick, and layers 53–60 are 100 ft thick (table 3).

The altitude of the base of the lowest model layer is –2,470 ft. 

This relatively fine vertical discretization was needed to 

accurately represent the hydrogeologic framework and to 

simulate the slope of the saltwater transition zones within 

( ) s s

C q q D C C C

t  

∂   = ∇⋅ ∇ − ∇⋅ −   ∂  



f C

C

  ∂

= +

∂

(5)

(6)

(7) ( ) ˆ

z q P pg e

μ

− = ∇ − 
Simulation of Transient Groundwater Flow 21

individual aquifers. A model cell was simulated as active if the 

GIS raster representation of a hydrogeologic unit indicated at 

least 1 ft of thickness at the center of the finite-difference cell. 

The GIS coverage of basement altitude was used to define the 

lowermost active layer in a particular row and column.

Internal Flow Package

The numerical representation of a hydrogeologic framework depends on the choice of method for simulating internal 

flow between finite-difference cells. Three MODFLOW 

packages exist that can be used to simulate internal flow 

with SEAWAT. The Block Centered Flow (BCF) and Layer 

Property Flow (LPF) packages are conceptually identical 

in their formulation of inter-cell conductance equations; 

their difference lies in details of model input, with the LPF 

package facilitating the use of parameters. These packages are 

well-suited to the approach of associating model layers with 

hydrogeologic units. The Hydrogeologic-Unit Flow (HUF) 

package (Anderman and Hill, 2000) is an alternative internalflow package for formulating the inter-cell conductance equations. Although the horizontal discretization of hydrogeologic 

units and finite difference cells is identical, the HUF package 

allows the flexibility of defining the layering and thickness of 

all hydrogeologic units independently from that of the model 

finite-difference-grid-cell layers. The HUF package internally 

calculates storage and inter-cell conductance terms for each 

cell from the hydraulic-property distributions assigned to 

individual hydrogeologic units present within that cell. Four 

considerations favored the use of the HUF package over either 

the BCF or LPF packages. 

1. The complexity of the hydrogeologic framework in the 

Virginia Coastal Plain rendered the “traditional” approach 

of associating model finite-difference-cell layers with 

individual hydrogeologic units impractical. As discussed 

under the hydrogeologic framework section of the report, 

no hydrogeologic unit is present throughout the entire 

groundwater model domain. All 19 units are present in 

only a portion of the model domain; they may be discontinuous, pinch out laterally, or be incised by river drainages. Some hydrogeologic units contain a “hole” in areas 

excavated during the Chesapeake Bay bolide impact, 

which is filled by different hydrogeologic units. Because 

of the limited spatial extent of the hydrogeologic units, 

their strict association with individual finite-differencecell layers was not possible. Approximating that association would have required a very complex finite-difference-cell layering scheme and was deemed impracticable.

2. In the Virginia Coastal Plain, many hydrogeologic units 

deepen and thicken toward the Atlantic coast and the 

zone of transition from freshwater to saltwater. Because 

the greater density of the saltwater affects the direction 

of groundwater flow, the transition zones are commonly 

delineated by equal concentration (iso-conc) surfaces that 

dip landward within individual aquifers. Under natural 

(steady-state) flow conditions, fresh groundwater discharges upward parallel to iso-conc surfaces. To simulate this variable-density flow effect in the groundwater 

model, the thicker, seaward portions of individual aquifers 

had to be represented with multiple layers of finitedifference cells. Using the HUF package conveniently 

allows improved discretization for numerical accuracy 

in these areas.

3. The ability to update the groundwater model if and 

when new hydrogeologic information is obtained was an 

important design goal. Model features such as observation 

wells and boundary conditions are referenced to the finitedifference grid and, therefore, require extensive modification if the model grid is modified. Because the extent 

and thickness of hydrogeologic units as well as their 

associated hydraulic-property distributions are defined 

in the HUF package independently from the model grid, 

they can be updated without changing observation data or 

boundary specifications. Consequently, the stated design 

goal is more effectively accomplished with the HUF package than the alternative internal flow packages. 

4. Because the hydrogeologic framework (McFarland and 

Bruce, 2006) was being developed contemporaneously 

with the groundwater model documented in this report, 

specification of the hydrogeologic-framework geometry 

independently from the other model features was crucial 

to timely completion of the model. During the course 

of model development, several interim versions of the 

hydrogeologic framework were used before the final 

version was incorporated in the model. Because these 

separate and distinct frameworks could be incorporated 

and tested without modifying the model grid, the required 

effort at construction of model input for observation wells 

or boundary conditions (pumpage, recharge, evapotranspiration, rivers or bay leakage, etc.) was economized.

Table 3. Model grid-layer 

thicknesses.

Layer 

numbers1

Cell 

thickness

(feet)

1 – 48 35 

49 – 52 50 

53 – 60 100 

1

Layers are numbered from 

the top to the bottom of the 

model.
22 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Temporal Discretization

The 113-year historical transient simulation was timediscretized into 34 stress periods of lengths indicated in 

table 4. An initial steady-state stress period was used to 

simulate predevelopment conditions. The time-step length 

within each transient stress period increased as a geometric 

progression with ratio 1.4. The number of time steps within 

each stress period was chosen such that the initial time-step 

length of each stress period was approximately 33 days. 

The hydraulic stresses in the second and third stress periods 

(table 4), which together represent the 30-year interval from 

1891 through 1920, are identical; this interval was divided 

to facilitate calculation of water-level observation times 

referenced from January 1, 1900. As explained later in the 

model-calibration section, the initial steady-state stress period 

was modified to transient to facilitate convergence during 

parameter-estimation iterations.

Boundary Conditions

Specified Flux Boundaries

The fraction of precipitation resulting in recharge to 

the water table and groundwater pumpage within the model 

domain are represented as specified flux boundaries. With 

two exceptions, no groundwater flow is simulated across the 

lateral and bottom boundaries of the model domain. The two 

exceptions where lateral boundary flow is non-zero simulate 

groundwater underflow: (1) to the north out of the model 

domain toward Maryland, and (2) from the Piedmont Province 

(west of the model domain) over the Norfolk Arch (fig. 1). 

Transient Areal Recharge

The correlation and similar magnitudes of annual 

precipitation measured at four stations in the Virginia Coastal 

Plain (fig. 8) indicate that the spatial variability of annual 

precipitation is relatively minor and likely less than the spatial 

variability of either runoff or vadose-zone evapotranspiration. 

Quantification of the individual magnitudes of either runoff 

or vadose-zone evapotranspiration or their spatial variability 

was beyond the scope of this study. Annual precipitation 

was, therefore, considered spatially homogeneous over the 

Coastal Plain, but variable in time. Recharge was simulated 

to be the fraction of this precipitation not lost to the combined 

magnitude of runoff and vadose-zone evapotranspiration. 

Recharge to the water table in the upper hydrogeologic 

unit (usually the surficial aquifer) was specified for the uppermost active model finite-difference cell using the Recharge 

(RCH) package of SEAWAT. This recharge flux simulates 

the fraction of the annual precipitation not lost to runoff or 

vadose-zone evapotranspiration and includes volumetrically 

minor amounts of septic-system effluents, return percolation 

from irrigation, and leakage from water and septic pipelines. 

Although this specified recharge flux is simulated as spatially 

homogeneous, it varies in time as a constant fraction of the 

measured annual Coastal Plain precipitation. This fraction was 

parameterized for model calibration and found to be correlated 

with parameters controlling the simulated evapotranspiration 

flux. Because this fraction could not be uniquely determined 

through model calibration, it was set to equal 50 percent of the 

measured annual precipitation. The variation of this simulated 

flux between stress periods is shown in figure 11. In low-lying 

areas, most of this recharge simulated within a model time step 

is lost from the groundwater system as simulated ET within 

the same model time step.

Table 4. Time discretization for 113-year historical transient 

simulation.

Stress

period

Year(s)

represented

Length 

in years

Length

in days

Time steps

in period

1 predevelopment - - 1 

2 1891–1899 9 3,287 11 

3 1900–1920 21 7,670 4 

4 1921–1939 19 6,939 13 

5 1940–1945 6 2,192 10 

6 1946–1952 7 2,557 11 

7 1953–1957 5 1,826 10 

8 1958–1964 7 2,557 11 

9 1965–1967 3 1,095 8 

10 1968–1972 5 1,827 10 

11 1973–1977 5 1,826 10 

12 1978–1981 4 1,461 9 

13 1982 1 365 5 

14 1983 1 365 5 

15 1984 1 366 5 

16 1985 1 365 5 

17 1986 1 365 5 

18 1987 1 365 5 

19 1988 1 366 5 

20 1989 1 365 5 

21 1990 1 365 5 

22 1991 1 365 5 

23 1992 1 366 5 

24 1993 1 365 5 

25 1994 1 365 5 

26 1995 1 365 5 

27 1996 1 366 5 

28 1997 1 365 5 

29 1998 1 365 5 

30 1999 1 365 5 

31 2000 1 366 5 

32 2001 1 365 5 

33 2002 1 365 5 

34 2003 1 365 5 
Simulation of Transient Groundwater Flow 23 Figure 11. Simulated recharge in the Virginia Coastal Plain by year.

0

0.008

0.007

0.0060.005

0.004

0.0030.002

0.001

steady

state

1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000

SIMULATED RECHARGE, IN FEET PER DAY

YEAR
24 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Transient Underflow to Maryland

Groundwater withdrawals from the Aquia aquifer in 

Maryland have resulted in water-level declines of over 200 ft 

in part of St. Marys County, MD. This area is included within 

the model domain, and withdrawals from the Aquia aquifer 

in that area were specified in the model. Withdrawals from 

the Upper and Lower Patapsco aquifers in Charles County, 

MD (which are correlative and laterally contiguous with the 

Potomac aquifer system in Virginia), have resulted in lower 

water levels in those aquifers. These withdrawals occur 

outside of the spatial domain of the Virginia Coastal Plain 

model (fig. 3) but affect water levels in simulated areas of 

Virginia and result in northward flow out of Virginia under the 

Potomac River into Maryland. To simulate declining water 

levels in Virginia resulting from this flow, a flux boundary was 

specified from the north side (row 1) of the model domain, 

from columns 5 to 37, layers 6 to 41, corresponding to the 

depths of the Potomac aquifer equivalent to the Upper and 

Lower Patapsco aquifers in Maryland (fig. 12a). This transient 

specified-flow was simulated with the Flow and Head Boundary (FHB) package (Leake and Lilly, 1997) and increased 

in a segmented linear fashion from zero before 1950 to a 

maximum of 8.86 x 105 ft3/d (6.63 Mgal/d) in 2003.

Constant Underflow from Piedmont

The magnitude of groundwater underflow from preCretaceous bedrock of the Piedmont Province into the Coastal 

Plain unconsolidated sediments is thought to be volumetrically 

minor. Under this assumption, the western model boundary 

was initially designated as no-flow. Analysis of water-level 

residuals during the course of model calibration suggested 

a component of inflow through the western boundary of the 

model. As discussed in the subsequent section entitled Boundary Flux Parameterization, the magnitude of flow through this 

area was parameterized and estimated in subsequent model 

calibration runs. In addition to the initial no-flow boundary 

specification, the hydraulic effects of possible underflow from 

the Piedmont were evaluated by testing three other specifiedflux boundary configurations, each representing different 

areas of underflow. In the calibrated model, the constant-flux 

boundary is specified along the west side in the western-most 

active cell from row 80 to 112 and model layer 7 to 16, 

corresponding to altitudes from zero to –350 ft (fig. 12b). The 

total simulated flow through this boundary is 2.16 x 106 ft3/d 

(16.2 Mgal/d).

Transient Groundwater Pumpage

Virginia Coastal Plain groundwater withdrawals are for 

municipal public-supply, industrial, commercial, agricultural, 

or domestic self-supply uses. Because the locations of municipal and industrial withdrawal wells are known and withdrawals from those wells are reported, these withdrawals can be 

simulated in space and time by flux assignment to appropriate 

finite-difference cells and stress periods. In contrast, the 

location and magnitude of self-supplied domestic groundwater 

withdrawals are uncertain. Reported and domestic withdrawals 

were simulated with different packages in SEAWAT and 

are discussed in the following subsections. Reported values 

for agricultural withdrawals are included in the withdrawal 

database and are simulated in the groundwater model. Because 

not all groundwater withdrawals for agricultural use are 

reported, not all agricultural withdrawals are simulated in the 

groundwater model. 

Reported Withdrawals

Groundwater withdrawal records were obtained from 

various sources for different areas and periods of time. In 

general, specificity with respect to individual wells improves 

for later simulated time periods. During earlier simulated 

times, withdrawals from multiple withdrawal wells in a well 

field may be combined in a single withdrawal record. Because 

the flux from all wells within the area of a finite-difference 

cell are simulated as being withdrawn from the center of 

the cell, this lack of well-specificity for earlier simulated 

times usually did not introduce any substantial additional 

numerical inaccuracy.

Site-specific groundwater withdrawal information for 

the years 1891 through 1981 were adapted from a previous 

compilation by Harsh and Laczniak (1990). Well-specific 

annual groundwater withdrawal information for the years 

from 1982 through 2003 was obtained from the Ground-Water 

Withdrawal Permitting and the Water Use/Water Supply 

Planning Programs of the Virginia DEQ (Robin Patton, written 

commun., 2002, 2003). Historical groundwater withdrawal 

data for wells in Maryland were obtained from the USGS 

Water Use Program (Judy Wheeler, written commun., 2005) 

and the Maryland Geological Survey (D.D. Drummond, written commun., 2005). Historical groundwater withdrawal data 

for wells in North Carolina were obtained from a compilation 

for the southern Atlantic Coastal Plain aquifer system (J. Fine, 

U.S. Geological Survey, written commun., 2005).

The 866 wells for which groundwater withdrawal records 

exist were simulated with the Multi-Node Well (MNW) package (Halford and Hanson, 2002). The locations of these wells 

are shown in figure 10. For reported withdrawals from wells 

with screens that penetrate more than one finite-difference-cell 

layer, the total reported annual withdrawal from each well 

is specified, but the layer-by-layer distribution of that flux 

depends, in part, on the simulated water-level differences 

between the withdrawal well and each of the screened 

finite-difference cells. Because of this head-dependency, these 

reported withdrawals may be considered a head-dependent 

flux for each layer. Water levels in withdrawal wells may be 

affected by head losses resulting from turbulent flow near 

the well or flow through drilling-damaged formation, the 

gravel pack, or the well screen; these effects were not directly 

simulated. The hydraulic conductance between simulated 

withdrawal wells and the finite-difference cells containing 

the well screen are calculated by the MNW package and 
Simulation of Transient Groundwater Flow 25

Figure 12. Finite-difference model grid showing cells designated for (a) time-varying specified flux to Maryland in row 1 and (b) constant 

specified flux from the Piedmont Province in columns 4–20.

(a)

(b)

MODEL ROW 1—NORTHERN BOUNDARY

MODEL COLUMN 1—WESTERN BOUNDARY
26 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

depend, in part, on the simulated well radius (rw) for each 

well. Because detailed well-construction information was not 

available for most reported withdrawal wells, the well radius 

was assumed to be similar for all wells and was estimated to 

be 0.55 ft during model calibration. Simulated water levels 

were insensitive to the parameter representing well radius.

Domestic Withdrawals

Unreported self-supplied domestic withdrawals comprise 

about 25 percent of the groundwater use in the Virginia 

Coastal Plain. Pope and others (2007) sampled private 

domestic well-completion records from Virginia Coastal 

Plain counties and independent cities to estimate the relative 

magnitudes of domestic groundwater withdrawals from 

hydrogeologic units within each county or independent city. In 

order to estimate the horizontal distribution of the magnitude 

of domestic withdrawals, road-density calculations were 

used to distribute the total census-block-group populations 

among the 1-mi2 model-grid cells within each census-block 

group. Using a combination of city boundaries and census 

spatial data on urban areas, grid cells thought to be publicly 

supplied were excluded to generate a horizontal coverage of 

model cells thought to be self supplied (fig. 10). Prior to 1950, 

all domestic-use groundwater was simulated as withdrawals 

from the surficial aquifer. For simulated times after 1950, the 

county- or independent-city-specific domestic-flux distribution 

from Pope and others (2007) was applied to apportion the total 

estimated domestic flux among the hydrogeologic units present within the area of each model cell. Within each cell area, 

the appropriate flux for each hydrogeologic unit was specified 

from the uppermost finite-difference-grid cell containing the 

hydrogeologic unit. The FHB package (Leake and Lilly, 1997) 

was used to simulate these domestic groundwater withdrawals 

and to temporally interpolate the change in flux from every 

specified finite-difference cell according to changes in decadal 

population-census data.

Head-Dependent Flux Boundaries

In contrast to the specified-flux boundaries, the simulated 

flow into or out of the groundwater model through headdependent-flux boundaries depends on the simulated water 

level adjacent to the boundary. 

Evapotranspiration

Evaporation and transpiration from the water table were 

simulated with the Evapotranspiration (EVT) package of 

SEAWAT. Because simulated stress periods were 1 year or 

longer in length, the maximum possible evapotranspiration 

flux was simulated to be constant at the potential evapotranspiration rate of 32.4 in/yr (0.0074 ft/d) for Virginia. This 

simulated maximum flux may occur when the water table 

is at or above the median land-surface altitude within each 

finite-difference cell (the specified evapotranspiration surface). 

The simulated evapotranspiration flux linearly decreases with 

depth of the simulated water table below this surface and 

becomes zero at the simulated evapotranspiration extinction 

depth. This depth, which was simulated to be uniform 

throughout the model domain, was estimated during model 

calibration to be 8.5 ft.

Rivers

Reaches of 32 rivers in the Coastal Plain for which the 

stage is above sea level were simulated with the River (RIV) 

package of SEAWAT. The total length of these reaches is 

792 mi. Reaches of these rivers with stage altitude less than 

1 ft were considered tidal and simulated with the General 

Head Boundary (GHB) package of SEAWAT. Utilizing the 

best available GIS coverage, the total length of mapped rivers 

within the model domain is approximately 19,000 mi. Thus, 

only 4.2 percent of the total river length is explicitly simulated 

in the model. (Because the larger river segments are simulated, 

the simulated riverbed seepage area is probably greater than 

4.2 percent of the actual riverbed area.) The discrepancy 

between simulated and actual seepage area available for interaction with the groundwater flow system may be compensated 

for numerically in two ways: (1) more groundwater discharge 

to model-represented river segments may be simulated than 

actually occurs for those segments, partially compensating for 

flux to river segments that are not represented in the model; 

and (2) in areas where the water table is sufficiently close to 

the land surface, discharge to many small streams may be 

accounted for by increased flux through the EVT package. 

These issues are discussed further in the water-budget section.

Both the RIV and GHB packages require specification 

of the hydraulic conductance and water level in the boundary 

feature being simulated for each model cell. Riverbed leakance (T –1) is the vertical hydraulic conductivity of the riverbed 

material divided by its thickness (Kz/b). The hydraulic 

conductance (dimension: [L2/T]) is the product of the leakance 

with the area of the boundary feature within the model cell. 

The water level in the river (river stage), river width, and the 

bottom altitude of the riverbed sediments were interpolated 

between river-property control points using the program 

RIVGRID (Leake and Claar, 1999), which also calculates the 

length of the river segment within each model cell. The number of control points used for each simulated river depicted in 

figure 7 are summarized in table 5. Riverbed leakance, which 

was assumed to be spatially homogeneous, was parameterized 

and estimated during model calibration to be 1 d–1.

Chesapeake Bay and Ocean Floor Leakage

Surface-water features with water-level altitudes less than 

1 ft, such as the Chesapeake Bay, Atlantic Ocean, and river 

tidal areas were simulated with the GHB package of SEAWAT. 

To determine the boundary conductance for each model cell, 
Simulation of Transient Groundwater Flow 27

the leakance parameter, which represents the vertical hydraulic 

conductivity of the bay or ocean floor material divided by 

its thickness (Kz/b), was multiplied by the area occupied by 

these features within each cell. Spatially uniform leakance 

parameters representing riverbed, bay, or ocean floor materials 

initially were parameterized separately for estimation. After 

boundary-leakance-parameter composite-scaled sensitivities 

indicated that these boundary leakance parameters were not 

sufficiently sensitive for independent estimation, they were 

combined into a single parameter representing riverbed, 

Chesapeake Bay, and ocean-floor materials.

Hydraulic-Property Distribution

Because hydraulic conductivity (K) and specific 

storage (Ss) depend on the intrinsic permeability and 

elasticity of the alluvial materials comprising the aquifer 

system, the spatial distribution of their magnitudes is 

expected to correlate with the lithologically defined 

hydrogeologic units. Parameters were, therefore, defined 

to represent the horizontal conductivity (Kh), vertical 

conductivity (Kv), and specific storage (Ss) within each 

hydrogeologic unit. Because the thickness (h) of every 

aquifer is spatially variable, the transmissivity (T = Khh) 

and storage coefficient (Ssh) of every aquifer varies spatially. Similarly, the thickness (b) of each confining unit 

and confining zone is spatially variable; consequently, 

the leakance (Kv/b) of these hydrogeologic units also 

varies spatially. Where observation data warranted a 

heterogeneous distribution of horizontal or vertical 

hydraulic conductivity within a hydrogeologic unit, this 

complexity was simulated using zonation, horizontal 

interpolation, or a depth-dependent function.

Horizontal Interpolation of 

Hydraulic Conductivity

To simulate horizontal heterogeneity of vertical or 

horizontal hydraulic conductivity within hydrogeologic 

units, 20 “pilot points” (Doherty, 2003) were placed 

within the model domain at locations depicted in 

figure 13. Hydraulic conductivities at these points were 

parameterized and estimated during model calibration. 

Hydraulic conductivities between pilot points were 

calculated by spatial interpolation of the native values 

of hydraulic conductivity estimated for the pilot points. 

The spatial interpolation was computed by ordinary 

kriging with an isotropic exponential variogram.

Depth-Dependence of Potomac Aquifer 

Horizontal Hydraulic Conductivity

The Potomac is a heterogeneous aquifer in 

which sharp changes in vertical and horizontal hydraulic 

conductivity (Kh) occur at interbed boundaries at a scale 

smaller than can be represented with the model discretization. 

At the model scale, a change of the average magnitude of Kh

with depth may also occur within the Potomac aquifer. This 

possibility was investigated using the hydraulic-conductivity 

depth-dependence capability (KDEP) of the HUF package 

(Anderman and Hill, 2003). The simulated horizontal 

Table 5. Number of control points for simulated rivers.

River, creek, reservoir, or swamp name Abbreviation

Number 

of control 

points

Appomattox River appo 10

Assamoosiak Swamp assa 2

Blackwater River blac 24

Cat Point Creek catp 3

Chickahominy River chic 4

Fontaine Creek font 2

North Landing River icwat 3

Kirby’s Creek kirb 2

Marracossic Creek marr 2

Mattaponi River matt 5

Meherrin River mehe 4

Mill Swamp mill 2

Nansemond River nans 7

North Anna River nora 9

Northwest River norw 2

Nottoway River nott 8

Piankatank River pian 2

Piscataway Creek pisc 3

Poni River poni 2

Potecasi Creek pote 2

Rappahannock River rapp 3

Lake Burnt Mills – Western Branch Reservoir rivb 7

Seacock Swamp seac 2

Somerton Creek some 2

South Anna River soua 4

Stony Creek ston 2

Swift Creek swif 3

Three Creek thre 2

Totoponomoy Creek toto 4

Wicomico River wico 2

James River jame 10

Rowanty Creek rowa 5
28 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 13. Pilot-point locations and simulated horizontal hydraulic conductivity at the top of the Potomac aquifer of the Virginia 

Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

HORIZONTAL HYDRAULIC

CONDUCTIVITY,

IN FEET PER DAY

EXPLANATION

0.0001

0.8 to 5.00

5.01 to 10.00

10.01 to 20.00

20.01 to 35.00

35.01 to 50.00

50.01 to 65.00

65.01 to 74.20

Aquifer not present

Pilot-point location

CHESAPEAKE

BAY

ATLANTIC

OCEAN

VIRGINIA

NORTH CAROLINA

38 ̊

37 ̊

77 ̊ 76 ̊
Model Calibration 29

hydraulic conductivity (Kh) was calculated for each model 

cell containing Potomac aquifer material by the exponential 

equation: ,

where

Kdepth is the horizontal hydraulic conductivity at 

depth d, 

Ksurface is the horizontal hydraulic conductivity at a 

reference altitude of 210 ft, 

λ is the depth-dependence coefficient, and 

d is the depth below the reference altitude of 

210 ft. 

The depth-dependence coefficient λ, which was parameterized 

and estimated during model calibration, was initially allowed 

to be either negative or positive to allow for the possibility of 

either a general increase or decrease of horizontal hydraulic 

conductivity with depth. The large composite scaled sensitivity of λ is due, in part, to the large number of Potomac aquifer 

water-level observations, which were used to estimate its 

value during the model-calibration regression. The simulated 

decrease of horizontal hydraulic conductivity with depth for 

the calibrated value of λ is depicted in figure 14. The highest altitude of the Potomac aquifer (along the Fall Line) 

is 174 ft, where the magnitude of the horizontal hydraulic 

conductivity is about 99 percent of the reference depth 

value (calculated from the horizontal interpolation between 

pilot points). The magnitude of the horizontal hydraulic 

conductivity simulated for the Potomac aquifer decreases to 

about 53 percent of the reference depth value in the bottom 

finite-difference-cell layer of the model. This decrease may 

represent either a general increase of sediment compaction 

with depth or a general change in sediment facies associated 

with the prograding fluvial-deltaic system that deposited the 

Potomac Formation.

Model Calibration

The model-calibration goal was to achieve the best fit 

of simulated to observed water levels (measured at discrete 

points in space and time) with as few parameters as practical. The parameter-estimation codes UCODE–2005 (Poeter 

and others, 2005) and PEST (Doherty, 2004) were used to 

calibrate horizontal and vertical hydraulic conductivity, 

specific storage, and boundary-flux model parameters to 

7,183 historical groundwater level observations in Virginia, 

Maryland, and North Carolina. Input to UCODE–2005 and 

PEST is similar, which facilitated utilization of both codes 

for comparative purposes and maximal advantage of their 

respective best features. Both PEST and UCODE–2005 

utilize non-linear regression to determine an optimal 

parameter vector b that minimizes an objective function S(b) 

of the form:

where 

hi is observed head, 

h ́

i is the simulated head, and 

wi is the weight for the observation. 

As suggested by Hill (1998), the weights wi were computed 

to be proportional to the inverse of the total error variance for 

each observation. The calculation of observation weights is 

described further in the subsequent section entitled “weights.”

Although only hydraulic heads were used as observations in the model-calibration regression, groundwater flux 

and age dates also constrained the calibration. Specification 

of groundwater withdrawals constrained estimated values 

of hydraulic conductivity in the confined aquifers, where 

pumpage and water-level measurements co-exist. The 

simulated river leakage was compared with discharge fluxes 

calculated by Richardson (1994b). Backward particle tracking 

with MODPATH (Pollock, 1994) was utilized to calculate 

simulated residence times for wells where carbon-14 (C-14) 

( ) 7183 2

1

( ) ( ) i i i iS b w h h b

=

= − ∑ ′

Figure 14. Decrease with depth of the simulated horizontal 

hydraulic conductivity of the Potomac aquifer of the Virginia 

Coastal Plain.

10 d

depth surface K K − = (8) (9)

–2,500

–2,000

–1,500

–1,000

–500

500

0

HYDRAULIC CONDUCTIVITY MULTIPLIER

ALTITUDE, IN FEET ABOVE NGVD 29

0.5 0.6 0.7 0.8 0.9 1.0
30 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

groundwater age measurements exist (Nelms and others, 

2003). Although river-leakage and C-14 data were not used in 

the model-calibration regression, similarity of the simulated 

equivalents to these data supported the general reasonableness 

of the conceptual flow model and values of simulated hydraulic conductivity.

Calibration Data Set

Within the area encompassed by the groundwater model 

domain, 497 observation wells in Virginia, Maryland, and 

North Carolina (fig. 15) were identified from the USGS 

Groundwater Site Inventory (GWSI). Historical measurements 

of groundwater levels in the identified wells, together with 

the associated accuracies of the observation-well location and 

the water-level measurement techniques, were extracted from 

GWSI and incorporated into a spreadsheet database. This 

water-level data set was subsequently modified to obtain a 

maximum of one measurement per year for each observation 

well. For wells with multiple observations in any year, the 

last observation in a calendar year was selected. The number 

of observation wells and associated water-level observations 

in each hydrogeologic unit that were used in the calibration 

process is summarized in table 6. The temporal distribution 

of water-level observations used in the calibration process is 

depicted in figure 16.

For the purpose of calibrating the groundwater model, 

observation wells ideally should have short screen intervals 

that enable water-level measurement in a single hydrogeologic 

unit. Water levels measured in wells with long screen intervals 

spanning multiple hydrogeologic units represent a composite 

of a potentially vertically varying head distribution of those 

hydrogeologic units. Comparison of simulated water levels 

to those observed in multi-aquifer wells is complicated by 

uncertainty of the influence of the head conditions in the 

different screened aquifer intervals to the water level in 

the observation well. The relative influence of the different 

screened intervals on the water-level measurement depends, 

in part, on hydraulic conductivity and water-level differences 

between the different screened hydrogeologic units. Aquifers 

containing sands or other high-permeability lithologies are 

expected to transmit fluid-pressure differentials more effectively than confining units composed predominantly of clays 

or other low-permeability lithologies. Consequently, the water 

level measured in a well with a screened interval penetrating 

both aquifer and adjacent confining units is generally biased 

toward the water level in the more highly permeable aquifer. 

To avoid the complications of multi-aquifer water-level 

observations and improve consistency among water-level 

observations used for model calibration, the Virginia Coastal 

Plain observation-well data set was culled to retain water-level 

observations from wells with screen intervals in a single 

aquifer. For those observation wells for which screen-interval 

information exists, the screened intervals of each observation 

well were compared with the hydrogeologic-unit top altitude 

from the digital representation of the hydrogeologic framework. Observation wells identified with screened intervals 

in more than one aquifer were removed from the observationwell data set. Observation wells with screened intervals in 

more than one hydrogeologic unit were retained if the additional hydrogeologic unit was classified as a confining unit.

The Head-Observation Process (Hill and others, 2000) of 

SEAWAT interpolates between finite-difference-cell centers to 

compute simulated equivalents at user-specified observationwell locations. The interpolation can fail for observation-well 

locations in finite-difference cells adjacent to no-flow boundaries or inactive finite-difference cells. Although SEAWAT 

continues to execute when this occurs, the correspondence 

of simulated water levels to observations is corrupted, which 

creates difficulties for parameter estimation with either 

UCODE–2005 or PEST. It was, therefore, necessary to 

remove these observation wells from the water-level observation data set.

Visual inspection of the time series of the observationwell water levels (hydrographs) revealed some problematic 

observation wells, in which hydrographs were either excessively noisy or contained measurements that were outliers 

to an overall water-level trend. The alpha-numeric identifier 

of each of these wells was catalogued according to the type 

of problem encountered. The problematic observations were 

either removed from the water-level observation data set or 

assigned a weight (so as to not unduly influence the parameterestimation regressions).

After editing, the calibration data set contained a total of 

7,183 observations of water levels from 497 observation wells, 

of which 5,087 observations are from 382 wells in Virginia, 

1,826 observations are from 89 wells in Maryland, and 

270 observations are from 26 wells in North Carolina.

During final model-calibration regression runs, optimization to the entire water-level data set, which included 1,826 

Maryland water-level observations, caused some degradation 

of model-fit to Virginia observations. An optimal fit to water 

levels in Virginia was considered more important than accurate 

simulation of conditions in Maryland, and the final calibration 

was performed using only the subset of Virginia water-level 

observations.
Model Calibration 31

Figure 15. Locations and screened hydrogeologic unit of water-level observation wells in the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

Figure 15. Locations and screened hydrogeologic unit of water-level observation wells in the

Virginia Coastal Plain.

Delmarva Peninsula

!

!

!

!!

!

! !

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!! !

!

!

! !

! !

!! !

! ! !

!!

!! !!

!

!

!

!!! !!!!

!!

!

! !

!! !!

"

"

"

""" "

"

" "

"

""

"

"

" "

!

!

!

!

! !

!

!!!

!

!

! !

!

!

!

!

!

!

!!

!

!! !!

!

!!!!

!

!

!! !

!

!

!

!!!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

! !

!

!

!

!

"

"

""

"

!

!

!!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!!

!!!!

!

!

!

!

!

!

!

!

!

!!

! !

!

! !

!

!

!

!

!

!

!

"

"

"

"

"

"

!

"

"

!

"

"

"

!

!

!

!

!

!

!

!

!

!

!

!

! !

! !

!

!

!

!

!

!

!

!

!

!

!

!

!!

!

!

!

!

!

!

!

!

! !!

!

!

!

!

!!

!

!

!

!

!

!

!!! !

!

! !

!!

!! !

!

!

!

!

! ! !

!

!

!

!

!

!

!

!

!!

!

!

!

!

!

!!

!

!!!!!

!

!

! !

!

!

! !

!

! !

!!

!!

!

!

!!

!

!

!

!

!

!

!

!

!

! !

! ! !!

!

!

!

!

! ! ! !

!

!

!

!

! !

!

!

!

!

!

!

!

!

! !

!

!

!

!

!

! !

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

! !

!!

! !

!

"

"

"

"

"

"

"

"

(

( (

(

(

(

(

(

(

(

(

(

(

(

(

( (

(

VIRGINIA

NORTH CAROLINA

MARYLAND

CHESAPEAKE

BAY

ATLANTIC

OCEAN

!

!

!

!

!!

!

!

! ! !

! !!! !!

! !

!

!

!!!

!

!

! !

!!

!

!

!

!

!

!

!

! !

!

!!

!

!

!

!

!

!

"

"

"

"

" "

"

Surficial aquifer

Yorktown confining zone

Yorktown-Eastover aquifer

St. Marys confining unit

St. Marys aquifer

Piney Point aquifer

Chickahominy confining unit

Exmore Matrix confining unit

Nanjemoy-Marlboro confining unit

Aquia aquifer

Peedee confining zone

Peedee aquifer

Virginia Beach confining zone

Virginia Beach aquifer

Upper Cenomanian confining unit

Potomac confining zone

Potomac aquifer

HYDROGEOLOGIC UNIT

EXPLANATION

!

"

!

!

"

!

!

"

"

"

!

"

!

"

"

!

37 ̊

38 ̊

77 ̊ 76 ̊
32 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 16. Time distribution and spatial distribution by State of water-level observations used for model calibration.

Table 6. Number of observations and wells primarily associated with each hydrogeologic 

unit simulated in the Virginia Coastal Plain model. 

Hydrogeologic 

unit1 Hydrogeologic unit name Abbreviation Wells2 Observations3

1 Surficial aquifer SURF 62 873

2 Yorktown confining zone YTCZ 6 110

3 Yorktown-Eastover aquifer YEAQ 59 833

4 Saint Marys confining unit SMCU 5 85

5 Saint Marys aquifer SMAQ 2 39

6 Calvert confining unit CACU 0 0

7 Piney Point aquifer PPAQ 49 523

8 Chickahominy confining unit CHCU 4 72

9 Exmore Matrix confining unit XMCU 0 0

10 Exmore Clast confining unit XCCU 2 28

11 Nanjemoy-Marlboro confining unit NMCU 17 136

12 Aquia aquifer AQAQ 63 969

13 Peedee confining zone PDCZ 1 26

14 Peedee aquifer PDAQ 1 25

15 Virginia Beach confining zone VBCZ 2 37

16 Virginia Beach aquifer VBAQ 1 16

17 Upper Cenomanian confining unit UCCU 3 53

18 Potomac confining zone POCZ 8 116

19 Potomac aquifer POAQ 212 3,242

1

Units are numbered from the top of the model to the bottom. 

2

Number of observation wells in unit.

3

Number of water-level measurements used for model calibration.

0

50

100

150

200

250

300

1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000

YEAR

NUMBER OF OBSERVATIONS

Maryland

North Carolina

Virginia
Model Calibration 33

Water-Level Observation Errors

According to Hill and Tiedeman (2007), “Observation 

error is error related to any aspect of the observation not 

accounted for by the model considered, for which the expected 

value is zero.” Quantification of the different sources of error 

in observations of water levels is necessary to calculate the 

weights used for model calibration. The sources of these errors 

are categorized under the following three types:

1. Observation-well location errors:

a) Altitude error: Water-level measurements typically 

are referenced to a measuring point on the well head 

near the land surface. The uncertainty associated 

with the measuring-point altitude varies among wells 

depending on the method used to determine the wellhead altitudes. Well-head altitudes in Virginia are 

determined by either a level survey or interpolation 

between the contours of a topographic map. Altitudes 

interpolated from topographic maps are considered to 

be accurate to within one-half of the contour interval 

of the map. For the 1:24,000 scale USGS topographic 

maps used in the Coastal Plain, the associated altitude 

uncertainties are either 2.5 ft or 5 ft. 

b) Horizontal-position error: Because water levels vary 

spatially, the uncertainty in the location (horizontal 

position) of an observation well results in further 

uncertainty in the water-level altitude associated 

with a specified horizontal location. Coastal Plain 

observation-well locations were determined either by 

a survey, such as with the Global Positioning System 

(GPS) or a theodolite, or by interpolation on a map. 

The expected accuracies of these methods range from 

0.01 second of arc (for some GPS measurements) to 

1 minute of arc (for some locations interpolated on 

a map). For this study, wells for which the positiondetermination method is unknown are assumed accurate to 1 minute of arc. A typical simulated waterlevel gradient of 5x10–4

 was used to translate these 

horizontal-position errors to errors in water-level 

altitude, which range from insignificant (for GPS 

measurements) to about 9 ft (for horizontal-position 

uncertainty of about 1 mi).

2. Water-level measurement errors: Most measurements of 

Virginia Coastal Plain water levels used for model calibration were made with steel measuring tapes (3,993 measurements). The method of measurement was unknown 

for 294 reported observations used in the calibration. 

Other types of measurements and the number of associated observations used in the calibration include graphic 

recorders (277 measurements), “electric tapes” (249 measurements), “calibrated electric tapes” (215 measurements), manometers (51 measurements), airlines (4 measurements), “pressure gauges” (2 measurements), and 

“calibrated pressure gauges” (1 measurement), and 1 measurement by a presumably known “other” method. Water 

levels measured with steel measuring tapes are assumed 

to be accurate to within 0.05 ft. Water levels measured by 

unknown methods are assumed to be accurate to within 

1 ft. Water levels measured by other methods are assumed 

to be accurate to within 0.1 ft.

3. Non-simulated transient-stress errors: Some observation 

wells respond to transient stresses not explicitly simulated 

in the groundwater model. While it can be argued that 

these observed responses are not errors in the observation, 

they are included as such in accordance with the opening 

definition of this section for the purpose of observation 

weighting. Examples of these stresses include (a) watertable fluctuations resulting from temporal variations in 

recharge and(or) evapotranspiration rates, (b) episodic 

pumping in nearby wells, (c) seasonal or episodic changes 

in rivers or stream stage, and (d) tides. The amplitude of 

water-level variation resulting from these stresses was 

estimated by comparing observed and simulated hydrographs for each observation well. Observed hydrographs 

were separated into “trend” and “noise” components, 

corresponding to water-level variations resulting from 

simulated and non-simulated processes. The amplitude of 

the interval containing 95 percent of the “noise” component was estimated for each observation hydrograph.

The magnitude of the individual error components for 

observation-well location and water-level measurements 

represent one-half of the width of their associated 95-percent 

confidence intervals. The error-component confidence 

intervals have a significance level of 5 percent, for which 

the critical value is 1.96. The component-error variances for 

each observation well were computed from their confidence 

intervals and the critical value, which were summed to obtain 

a total variance associated with each water-level observation. 

Weights

In general, the weights (wi) assigned to water-level 

observations were calculated as the inverse of the total 

observation-error variance. The total water-level observation 

error variance is the sum of the component-error variances 

described under the sub-heading “Water-Level Observation 

Errors.” Because UCODE–2005 and PEST were used to 

perform non-linear regressions for parameter estimation, the 

square root of the weight (square root of the inverse of the 

total observation error variance) was specified for both program input files to facilitate comparison of objective function 

values between the two different programs. The magnitudes of 

observation weights ranged from essentially zero to about 10. 

Most observations (91.5 percent) had weights between 0.1 and 

0.4; and 6.5 percent had weights between 0.4 and 1.1.
34 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Modification of First Stress Period TimeDiscretization for Parameter Estimation

Some non-linear effects of head-dependent boundary 

conditions, particularly the simulation of evapotranspiration, 

can cause numerical instability resulting in nonconvergence 

for the initial steady-state solution used to simulate predevelopment conditions. For this model, the convergence problem 

typically occurred when parameters associated with the 

boundary condition (such as the ET extinction depth) were 

adjusted in the regression process used for model calibration. 

To enable the regression to proceed in these instances, a 

100-year transient stress period with no hydraulic stress was 

substituted for the initial steady-state stress period. This 

surrogate predevelopment transient stress period consisted of 

60 time steps that increased in length in a geometric progression of ratio 1.2, such that the initial time step was approximately 3 hours in length. The storage term in the transient flow 

equation and small initial time-step length prevented the large 

head oscillations related to numerical instability and nonconvergence for a steady-state solution. The long (100‐year) 

period with no hydraulic stress allowed specified starting 

water levels sufficient time to adjust to the updated hydraulicparameter distribution and boundary conditions for each 

parameter-estimation iteration, and approached steady-state 

conditions in 100 years. The water-level distribution at the 

end of this 100-year predevelopment stress period was saved 

for subsequent use as the starting water-level distribution in 

later steady-state model solutions computed between regression runs. These steady-state solutions were used to update 

regression-run starting water-level distributions in order to 

minimize the possible effect of transient residual starting water 

levels (calculated from a previous parameter distribution) in 

low-permeability hydrogeologic units.

Hydraulic Property Parameterization

A parsimonious approach was used for assignment of 

parameters representing hydraulic conductivity and specific 

storage in sub-domains of the model. Each hydrogeologic unit 

was initially assigned a unique parameter representing vertical 

hydraulic conductivity, horizontal hydraulic conductivity, 

and specific storage. Computation of the composite scaled 

sensitivities (CSS) of each of these model parameters identified parameters that could be effectively estimated through 

non-linear regression. The magnitude of composite scaled 

sensitivities for these initial parameters was largely dependent 

on the number of water-level observations within the corresponding hydrogeologic unit. Where possible, parameters 

with small composite scaled sensitivities were combined in 

a manner consistent with current knowledge of the hydraulic 

properties of the corresponding hydrogeologic units. For 

example, with the exception of the Potomac confining zone, 

the composite scaled sensitivities of parameters representing 

horizontal hydraulic conductivity in confining units and 

zones were low. On the basis of the initial composite scaled 

sensitivity and lithology of these units, they were subsequently 

sorted into four groups: (1) the Yorktown confining zone, 

(2) the Potomac confining zone, (3) the Chickahominy and 

Exmore confining units, and (4) all other confining units and 

zones. A single parameter was used to represent the horizontal 

hydraulic conductivity of each group. In later parameterestimation runs, the continued small composite scaled 

sensitivity of the parameter representing horizontal hydraulic 

conductivity of the Yorktown confining zone suggested further 

modification of this parameterization scheme, after which 

it was combined with the parameter representing horizontal 

hydraulic conductivity of the Yorktown-Eastover aquifer. A 

similar process of composite scaled sensitivity evaluation 

and re-parameterization was conducted for those parameters 

representing vertical hydraulic conductivity in each of the 

aquifers. With the exception of the Potomac aquifer, all 

parameters representing vertical hydraulic conductivity in 

individual aquifers had low composite scaled sensitivities. 

Following several iterations of re-parameterization and testing, 

these individual parameters were combined into three groups 

representing the vertical hydraulic conductivity in (1) the 

surficial and Yorktown-Eastover aquifers, (2) the Potomac 

aquifer, and (3) all other aquifers.

For each new parameterization scheme, the model fit was 

optimized by minimizing the objective function with either 

UCODE–2005 or PEST. Observation residual maps for each 

hydrogeologic unit were subsequently generated and inspected 

to identify (1) large residuals and (2) non-random spatial 

distribution of the positive and negative weighted residuals. 

Hydrographs of the observations used in the calibration 

along with their simulated equivalents were generated and 

inspected following each parameter-estimation run. Inspection 

of the hydrographs containing large residuals was useful for 

diagnosing probable causes of inadequate model fit in some 

areas. If parameters with large composite scaled sensitivities 

represented the hydraulic conductivity or specific storage for a 

hydrogeologic unit in an area of apparent non-random residual 

distribution, one or more additional parameters were added to 

represent the hydraulic conductivity or specific storage in that 

area. The area of influence of the additional parameter was 

defined either by a refined area zonation or association with an 

additional pilot point, which effectively modified the spatial 

interpolation of the hydraulic parameter. 

Because aquifer water levels are also affected by leakage through adjacent confining units, the vertical hydraulic 

conductivity (Kv) of confining units may also affect the pattern 

of residuals in aquifers. Because water levels within the Piney 

Point aquifer were sensitive to the vertical hydraulic conductivity of the Calvert confining unit, the spatial variability of 

vertical hydraulic conductivity within the Calvert confining 

unit (fig. 17) was simulated with pilot points.
Model Calibration 35

Figure 17. Estimated vertical hydraulic conductivity of the Calvert confining unit of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

1.0x10–6 to 1.2x10–5

1.3x10–5 to 2.0x10–5

2.1x10–5 to 2.7x10–5

2.8x10–5 to 3.2x10–5

3.3x10–5 to 3.5x10–5

VERTICAL HYDRAULIC

CONDUCTIVITY,

IN FEET PER DAY

Richmond

Williamsburg

Suffolk

Norfolk

Fredericksburg

EXPLANATION

VIRGINIA

NORTH CAROLINA

CHESAPEAKE

BAY

ATLANTIC

OCEAN

MARYLAND

37 ̊

38 ̊

77 ̊ 76 ̊
36 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Boundary Flux Parameterization

Early in the calibration process, water-level residual maps 

showed an area around the updip Norfolk Arch (fig. 1) of 

non-random positive residuals, indicating that observed water 

levels were greater than simulated water levels. This clustering 

suggested insufficient groundwater flux to the Potomac aquifer 

system in that area. Existing boundary-condition specifications 

(particularly specified river stages) were verified, and various 

parameterizations of overlying confining units were tested. 

Two alternative approaches were tested to achieve improvement of the model fit in this area: (1) decrease the horizontal 

hydraulic conductivity to the west (which increased the 

hydraulic gradient), or (2) allow some recharge as underflow 

from the Piedmont groundwater flow system over the Norfolk 

Arch. Following the first approach, parameter optimization 

estimated low (1–2 ft/d) hydraulic conductivities in the area. 

For the second approach, a zone of underflow from the west 

was defined (fig. 12b) over the Norfolk Arch through which 

specified flux was parameterized. Sensitivity analysis indicates 

that water-level observations are sensitive to the parameter 

representing this flux. Parameter optimization estimated an 

underflow equivalent to approximately 0.2 percent of precipitation flux over the area of the Piedmont west of the flux 

boundary. This (underflow) approach was considered more 

plausible than the alternative (low hydraulic conductivity) 

approach and was retained for the final model configuration.

Simulated Conditions

The goal of the calibration was to simultaneously meet 

the following conditions: (1) The non-linear regression converged with a minimum sum-of-squares error (SSE), (2) the 

regression-determined parameter values were within the range 

of pre-defined reasonable values, (3) parameters were not 

excessively correlated, (4) the simulated hydraulic-property 

distribution within hydrogeologic units was reasonable for the 

model scale, and (5) positive and negative model residuals 

were randomly distributed in space and time. The model was 

considered calibrated when these criteria, which are discussed 

further in the subsequent section entitled “Calibration Assessment” were met as closely as possible.

Water Levels

The 113-year transient groundwater model simulates the 

changing distribution of water levels throughout the aquifer 

system in space and time. Maps of simulated water levels 

along any altitude or surface of interest within the aquifer 

system may be generated from model output corresponding 

to the end of any model-simulated time step. Maps of water 

levels within hydrogeologic units designated as aquifers are 

useful for determining directions of flow and assessing aquifer 

impacts resulting from groundwater development. The maps 

of simulated water levels in this report illustrate the simulated 

altitudes of the water table, which are generally in the surficial 

aquifer and water levels in the middle of the YorktownEastover, St. Marys, Piney Point, Aquia, Peedee, and Virginia 

Beach aquifers, and at the top of the Potomac aquifer. Because 

of vertical gradients within the thick, vertically anisotropic 

Potomac aquifer, maps of water levels differ at various 

depths within the Potomac aquifer, particularly in areas near 

large withdrawal wells during transient stress periods. The 

potentiometric surface at the top of the Potomac aquifer was 

selected for depiction because conditions at the top of the 

Potomac aquifer are most likely to be evaluated for groundwater withdrawal-permit applications (Robin Patton, Virginia 

Department of Environmental Quality, oral commun., 2006). 

In subsequent sections, simulated water levels corresponding 

to predevelopment conditions are depicted and discussed for 

each aquifer, and water levels simulated for the year 2003 

are depicted and discussed for the seven aquifers beneath 

the surficial aquifer. Although simulated water levels within 

the surficial aquifer are generally similar to those simulated 

for predevelopment conditions, the drawdown simulated in 

the surficial aquifer for the drought year 2001 is depicted. 

The simulated water-level change between predevelopment 

and 2003 for the Yorktown-Eastover, St. Marys, Piney Point,

Aquia, and Potomac aquifers is also depicted as drawdown 

maps. No substantial regional drawdown within the Peedee 

and Virginia Beach aquifers was simulated and, therefore, is 

not depicted for those aquifers.

An estimate of the magnitude of water-level change 

(drawdown) due to domestic withdrawals within each aquifer 

was made by assuming model linearity and application of the 

principal of superposition. Although the model is somewhat 

non-linear, any error in this calculation from the assumption 

of model linearity is smaller than that from the uncertainty of 

domestic withdrawal magnitudes.

Because the model was calibrated to water levels 

measured at discrete points in space and time, comparison of 

hydrographs containing the observations used for calibration 

with their simulated equivalents provides a visual indication 

of the quality of the model fit in various aquifers and areas of 

the Coastal Plain. In addition to depicting the quality of the 

model fit, hydrographs of simulated and observed water levels 

illustrate the rate of change in water levels at various points in 

the aquifer system through time. Representative hydrographs 

of water-level measurements with their simulated equivalents 

are shown and discussed in subsequent sections for seven of 

the eight aquifers. Water-level observations for the Virginia 

Beach aquifer were not available.
Simulated Conditions 37

Surficial Aquifer

The simulated altitude of the water table in the Virginia 

Coastal Plain for the steady-state stress period representing 

predevelopment conditions is shown in figure 18a. Although 

the water table is generally located within the surficial aquifer, 

there are some areas where the simulated water table resides 

within one of the underlying aquifers or confining units. 

Because of the strong influence of topography, the watertable altitude varies with a higher spatial frequency than the 

potentiometric surface of any underlying confined aquifer. 

Simulated evapotranspiration largely controls the depth of 

the simulated water table, resulting in a water-table surface 

that is a smoothed and subdued replica of the land-surface 

altitude specified in the model EVT package. The simulated 

water-table altitudes depicted in figure 18a suggest that 

shallow groundwater flows downward away from higheraltitude areas, such as near the Fall Line and the mid-lines of 

major peninsular interfluves, toward lower altitudes adjacent 

to simulated rivers. However, most of the water recharged to 

the surficial aquifer likely discharges to one of the adjacent 

smaller streams (fig. 7) that were not explicitly simulated. As 

discussed in the subsequent section entitled “Water Budget,” 

detailed simulation of surficial-aquifer flow from interfluves to 

the non-simulated streams is not possible with the resolution 

of the finite-difference grid of this regional model. 

Differences in simulated water-table altitudes between 

stress periods, including the differences between steady-state 

and 2003, are predominantly because of simulated differences 

in recharge between stress periods. Simulated water-table 

altitudes at the end of the simulation in 2003 are within several 

feet of the predevelopment water-table altitudes and are, 

therefore, not depicted in a figure. The below-normal precipitation over the Virginia Coastal Plain during 2001 (fig. 8) was 

simulated with proportionately smaller recharge, resulting in 

lower simulated water-table altitudes for 2001. The difference 

between water-table altitudes simulated for predevelopment 

conditions and the drought year 2001 is depicted in figure 18b.

The locations of observation wells in the surficial, 

Yorktown-Eastover, St. Marys, Piney Point, Aquia, and 

Peedee aquifers for which hydrographs are depicted are shown 

in figure 19. In the example hydrograph from surficial-aquifer 

observation well 57D 23 (fig. 20), the simulated water 

level responds to simulated changes in annual recharge and 

reasonably mimics observed annual changes in water levels. A 

superposition analysis indicated that the drawdown resulting 

from domestic pumpage with a constant recharge rate was less 

than 1 ft throughout the surficial aquifer.
38 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 18a. Simulated predevelopment water table in the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

WATER-TABLE ALTITUDE,

IN FEET

EXPLANATION

–2.1 to 0

0.1 to 20.0

20.1 to 40.0

40.1 to 60.0

60.1 to 80.0

80.1 to 100.0

100.1 to 120.0

120.1 to 140.1

140.1 to 160.0

160.1 to 180.0

180.1 to 351.6

CHESAPEAKE

BAY

ATLANTIC

OCEAN

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
Simulated Conditions 39

Figure 18b. Simulated water-table drawdown in 2001 in the Virginia Coastal Plain.

DRAWDOWN, IN FEET

EXPLANATION

–0.18 to 0

0.01 to 0.50

0.51 to 1.00

1.01 to 1.50

1.51 to 2.00

2.01 to 2.50

2.51 to 3.96

0

0

25 50 MILES

25 50 KILOMETERS

CHESAPEAKE

BAY

ATLANTIC

OCEAN

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
40 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 19. Locations of observation wells screened in the surficial, Yorktown-Eastover, St. Marys, Piney Point, Aquia, 

and Peedee aquifers for hydrographs depicted in figures 20, 22, 24, 26, 28, and 30.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

Delmarva Peninsula

VIRGINIA

NORTH CAROLINA

CHESAPEAKE

BAY

ATLANTIC

OCEAN

Figure 19. Locations of observation-well hydrographs from the surficial, Yorktown-Eastover, St. Marys, 

Piney Point, Aquia, and Peedee aquifers depicted in figures 19, 21, 23, 25, 27, and 29.

57G 2

57H 14

56G 38

56H 29

57D 23

51M 18

54Q 73

62A 2

62C 9

61B 6

61B 2

55B 67 57B 8

57C 21

Surficial

Yorktown-Eastover

St. Marys

Piney Point

Aquia

Peedee

SCREENED AQUIFER

EXPLANATION

37 ̊

38 ̊

77 ̊ 76 ̊
Simulated Conditions 41

Figure 20. Observed and simulated water levels in the surficial aquifer well 57D 23 in the 

Virginia Coastal Plain. (Well location is shown in figure 19.)

Yorktown-Eastover Aquifer

The simulated steady-state water levels in the YorktownEastover aquifer (fig. 21a) indicate that predevelopment 

groundwater in the Yorktown-Eastover aquifer in Virginia 

flowed away from higher-altitude areas, near the midlines 

of the Northern Neck, Middle Peninsula, and York-James 

Peninsula, toward lower-altitude discharge areas, including 

major rivers and the Chesapeake Bay. The simulated water 

levels in the Yorktown-Eastover aquifer for 2003 (fig. 21b) 

suggest that regional groundwater flow during 2003 was 

similar to that simulated for predevelopment conditions. The 

drawdown of groundwater levels in the Yorktown-Eastover 

aquifer, derived by subtracting simulated 2003 water levels 

from predevelopment water levels, is shown in figure 21c. 

The drawdown under low-altitude areas, such as major rivers 

and the Chesapeake Bay, typically is less than 1 ft and likely 

results from decreased groundwater discharge through the 

Yorktown-Eastover aquifer due to groundwater withdrawals 

from the regional aquifer system. Simulated water-level 

drawdowns under land areas, resulting from localized 

groundwater withdrawals from the Yorktown-Eastover aquifer, 

generally are less than 2 ft, with a maximum of 2.77 ft in the 

southeastern part of the city of Chesapeake. Large land areas 

with up to 1 ft of simulated groundwater-level recovery, result 

from transient recharge greater than the long-term average 

value simulated during the steady-state stress period.

Wells 57H 14, 61B 2, 62A 2, 62C 9, and 57B 8 are 

screened in the Yorktown-Eastover aquifer at various locations 

in the Coastal Plain (fig. 19). The hydrograph in figure 22 

depicts water levels measured in these wells along with their 

corresponding simulated values.

60

70

72

74

76

78

62

64

66

68

80

1980 1985 1990 1995 2000

Simulated 57D 23

Observed 57D 23

YEAR

WATER-LEVEL ALTITUDE, IN FEET
42 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 21a. Simulated predevelopment water levels in the Yorktown-Eastover aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

Figure 21a. Simulated predevelopment water level in the Yorktown-Eastover aquifer of the Virginia Coastal Plain.

–2 to 0

1 to 10

11 to 20

21 to 40

41 to 60

61 to 80

81 to 100

101 to 120

121 to 158

WATER LEVEL,

IN FEET

EXPLANATION

VIRGINIA

NORTH CAROLINA

Northern Neck

Middle

Peninsula

York-James

Peninsula

CHESAPEAKE

BAY

ATLANTIC

OCEAN

MARYLAND

WATER LEVEL—Contour

interval 20 feet

20

80

60

80

80

100

100

60

60

20

20

80

80

80

80

60

60

60

60

40

40

40

20

20

37 ̊

38 ̊

77 ̊ 76 ̊
Simulated Conditions 43

Figure 21b. Simulated 2003 water levels in the Yorktown-Eastover aquifer of the Virginia Coastal Plain.

80

60

80

80

100

100

60

60

20

20

80

80

80

80

60

60

60

40

40

40

20

20

MARYLAND

–10 to 0

1 to 10

11 to 20

21 to 40

41 to 60

61 to 80

81 to 100

101 to 120

121 to 160

WATER LEVEL, IN FEET

EXPLANATION

VIRGINIA

NORTH CAROLINA

Northern Neck

Middle

Peninsula

York-James

Peninsula

0

0

25 50 MILES

25 50 KILOMETERS

Figure 21b. Simulated 2003 water level in the Yorktown-Eastover aquifer of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

WATER LEVEL—Contour

interval 20 feet

20

37 ̊

38 ̊

77 ̊ 76 ̊
44 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 21c. Simulated drawdown from predevelopment to 2003 in the Yorktown-Eastover aquifer of the Virginia Coastal Plain.

–1

MARYLAND

CHESAPEAKE

BAY

ATLANTIC

OCEAN

Northern Neck

Middle

Peninsula

York-James

Peninsula

0

0

25 50 MILES

25 50 KILOMETERS

Figure 21c. Simulated 2003 drawdown in the Yorktown-Eastover aquifer of the Virginia Coastal Plain.

VIRGINIA

NORTH CAROLINA

0

–1

–1

0

–1

–1

0

0

–1

–1

–1

–1

–1

0

–1

–1

–1

–1

–3.41 to –1.0 

–0.99 to 0

0.01 to 1.0

1.01 to 3.0

DRAWDOWN, IN FEET

EXPLANATION

CHANGE IN WATER LEVEL

 —Contour interval 1 foot

–1

37 ̊

38 ̊

77 ̊ 76 ̊
Simulated Conditions 45

Figure 22. Observed and simulated water levels in wells 57H 14, 61B 2, 62A 2, 62C 9, and 57B 8 in the Yorktown-Eastover 

aquifer of the Virginia Coastal Plain. (Well locations are shown in figure 19.)

St. Marys Aquifer

In southern Virginia, the St. Marys aquifer is present as 

a north-south trending sand layer approximately 50 ft thick 

between the underlying Calvert confining unit and the overlying St. Marys confining unit. Predevelopment water levels 

in the St. Marys aquifer simulated by the steady-state stress 

period are depicted in figure 23a. The simulated predevelopment water-level contours suggest that groundwater flowed 

generally from west to east through the St. Marys aquifer in 

Virginia and south into North Carolina. Water levels simulated 

in the St. Marys aquifer for 2003 (fig. 23b) are lower, but 

suggest flow directions similar to those simulated for predevelopment time. The simulated water-level drawdown from 

predevelopment through 2003 exceeds 40 ft near the northern 

terminus of the St. Mary’s aquifer in southern Virginia 

(fig. 23c). 

The St. Marys aquifer, (along with the Peedee aquifer), 

is one of the least-utilized aquifers in southern Virginia; the 

relatively small quantity of withdrawn groundwater (table 2) 

is used for private water supplies in the city of Suffolk. In 

comparison with other developed Coastal Plain aquifers, the 

St. Marys aquifer in southern Virginia has relatively limited 

area and low transmissivity (fig. 45b, p. 89), which may 

account for the disproportionately large simulated drawdown 

resulting from the relatively modest groundwater withdrawals 

from the St. Marys aquifer. A superposition analysis indicates 

that simulated domestic withdrawals account for up to 3 ft of 

the drawdown in the St. Marys aquifer in southern Virginia. 

The location of a water-level observation well (57C 21) 

that is screened in the St. Marys aquifer in southern Virginia is 

depicted in figure 19. Figure 24 depicts water levels measured 

in this well, along with the corresponding simulated values.

0

10

20

30

40

50

60

Simulated 57H 14

Observed 57H 14

Simulated 61B 2

Observed 61B 2

Simulated 62A 2

Observed 62A 2

Simulated 62C 9

Observed 62C 9

Simulated 57B 8

Observed 57B 8

1975 1980 1985 1990 1995 2000

YEAR

WATER-LEVEL ALTITUDE, IN FEET
46 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 23a. Simulated predevelopment water levels in the St. Marys aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

Figure 23a. Simulated predevelopment water level in the St. Marys aquifer of the Virginia Coastal Plain. 40

30

10

20

30

10

VIRGINIA

NORTH CAROLINA

MARYLAND

–5 to 0

1 to 10

11 to 20

21 to 30

31 to 40

41 to 52

WATER LEVEL, IN FEET

EXPLANATION

WATER LEVEL—Contour

interval 10 feet

10

CHESAPEAKE

BAY

ATLANTIC

OCEAN

37 ̊

38 ̊

77 ̊ 76 ̊
Simulated Conditions 47

Figure 23b. Simulated 2003 water levels in the St. Marys aquifer of the Virginia Coastal Plain.

10

0

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

Figure 23b. Simulated 2003 water levels in the St. Marys aquifer of the Virginia Coastal Plain. 10

20

10

VIRGINIA

NORTH CAROLINA

10

WATER LEVEL, IN FEET

EXPLANATION

WATER LEVEL—Contour

interval 10 feet

10

–27 to –10

–9 to 0

1 to 10

11 to 20

21 to 30

31 to 40

41 to 51

CHESAPEAKE

BAY

ATLANTIC

OCEAN

37 ̊

38 ̊

77 ̊ 76 ̊
48 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 23c. Simulated drawdown from predevelopment to 2003 in the St. Marys aquifer of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

Figure 23c. Simulated 2003 drawdown in the St. Marys aquifer of the Virginia Coastal Plain.

40

VIRGINIA

NORTH CAROLINA

40

30

20

10

10

0

–2 to 0

1 to 10

11 to 20

21 to 30

31 to 40

41 to 51

DRAWDOWN, IN FEET

EXPLANATION

CHANGE IN WATER LEVEL

 —Contour interval 10 feet

10

CHESAPEAKE

BAY

ATLANTIC

OCEAN

37 ̊

38 ̊

77 ̊ 76 ̊
Simulated Conditions 49

Figure 24. Observed and simulated water levels in well 57C 21 in the St. Marys aquifer of the Virginia Coastal Plain. 

(Well location is shown in figure 19.)

Piney Point Aquifer

Although the Piney Point aquifer is regionally extensive 

beneath the Coastal Plain, it is relatively thin south of the 

James River, with correspondingly low transmissivity. The 

simulated potentiometric surface of the Piney Point aquifer 

during the steady-state stress period is shown in figure 25a. 

Predevelopment groundwater flowed from western up-dip 

areas toward the east and southeast, to discharge through 

the overlying Calvert confining unit near the Atlantic Coast. 

By the year 2003, concentrated groundwater withdrawals 

resulted in several cones of depression (fig. 25b) that altered 

flow directions in the Piney Point aquifer in the vicinity of the 

Middle and York-James Peninsulas. The maximum simulated 

drawdown in the Piney Point aquifer near Williamsburg 

is 127 ft (fig. 25c). A superposition analysis indicates that 

simulated domestic withdrawals may account for approximately 10 ft of drawdown in the Piney Point aquifer near 

Williamsburg, and from 7 to 9 ft across the eastern portions of 

the Middle Peninsula and Northern Neck.

A hydrograph for observation wells 54Q 73, 56G 38, 

and 56H 29 in the Piney Point aquifer is depicted in figure 26. 

Except for a few observations near Williamsburg, the fit of 

simulated to observed water levels in Piney Point observation 

wells is good. Simulated water levels in the Piney Point 

aquifer with less-than-optimal fit to observed water levels 

near Williamsburg were sensitive to magnitudes of vertical 

hydraulic conductivity of the adjacent Calvert and NanjemoyMarlboro confining units. Although further refinement of 

the simulated spatial variation of the Calvert confining unit 

vertical hydraulic conductivity could improve the model fit 

for some observation wells, such an adjustment was limited 

during model calibration because of uncertainty in the 

magnitudes of domestic withdrawals in the Williamsburg area. 

An informal polling of domestic well drillers (Scott Bruce, 

Virginia Department of Environmental Quality, oral commun., 

2006) suggested that actual domestic withdrawals from the 

Piney Point aquifer in the Williamsburg area could be substantially greater than the simulated fluxes, which were based on 

the analysis of Pope and others (2007). Superposition analysis 

of drawdown rates resulting from the specified domestic 

withdrawals suggested that under-representation of the 

domestic-withdrawal flux may also account for the magnitude 

of differences between water levels simulated and observed in 

some Piney Point observation wells.

1975 1980 1985 1990 1995 2000

0

10

20

22

24

26

28

12

14

16

18

2

4

6

8

30

Simulated 57C 21

Observed 57C 21

YEAR

WATER-LEVEL ALTITUDE, IN FEET

Figure 24. Observed and simulated water levels in well 57C 21 in the St. Marys aquifer of the Virginia Coastal Plain.

Well location is shown on figure18.
50 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 25a. Simulated predevelopment water levels in the Piney Point aquifer of the Virginia Coastal Plain.

MARYLAND

10

20

10

40

50

60

70

90

30

10

30

30

40

50

60

20

10

40

50

60

20

10

30

50 70

0

0

25 50 MILES

25 50 KILOMETERS

Figure 25a. Simulated predevelopment water levels in the Piney Point aquifer of the Virginia Coastal Plain.

VIRGINIA

NORTH CAROLINA

0 to 20

21 to 40

41 to 60

61 to 80

81 to 100

101 to 117

WATER LEVEL,

IN FEET

EXPLANATION

WATER LEVEL—

Contour interval

10 feet

10

CHESAPEAKE

BAY

ATLANTIC

OCEAN

37 ̊

38 ̊

77 ̊ 76 ̊
Simulated Conditions 51

Figure 25b. Simulated 2003 water levels in the Piney Point aquifer of the Virginia Coastal Plain.

-

-

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

Figure 25b. Simulated 2003 water levels in the Piney Point aquifer of the Virginia Coastal Plain.

20

10

VIRGINIA

NORTH CAROLINA

20

10

10

20

20

40

30

50

60

10

30

40

80

–30

–40

–50–60

–20

–10

–30

–10

0

10

20

0

0

0

0

10

WATER LEVEL,

IN FEET

EXPLANATION

WATER LEVEL—

Contour interval 

10 feet

10

–93 to –80

–79 to –60

–59 to –40

–39 to –20

–19 to 0

1 to 20

21 to 40

41 to 60

61 to 80

81 to 100

101 to 117

CHESAPEAKE

BAY

ATLANTIC

OCEAN

37 ̊

38 ̊

77 ̊ 76 ̊
52 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 25c. Simulated drawdown from predevelopment to 2003 in the Piney Point aquifer of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

Figure 25c. Simulated 2003 water levels drawdowns in the Piney Point aquifer of the Virginia Coastal Plain.

VIRGINIA

NORTH CAROLINA

0

DRAWDOWN, IN FEET

EXPLANATION

CHANGE IN WATER LEVEL—

 Contour interval 10 feet

10

–2 to 10

11 to 20

21 to 30

31 to 40

41 to 50

51 to 60

61 to 70

71 to 80

81 to 90

91 to 100

101 to 110

111 to 130

CHESAPEAKE

BAY

ATLANTIC

OCEAN

40

20

30

20

10

10

20

30

40

50

60

70

80

10

20

30

40

40

60

70

37 ̊

50

38 ̊

77 ̊ 76 ̊
Simulated Conditions 53

Figure 26. Observed and simulated water levels in wells 56H 29, 54Q 73, and 56G 38 in the Piney Point aquifer 

of the Virginia Coastal Plain. (Well locations are shown in figure 19.)

Aquia Aquifer

The predevelopment groundwater flow within the Aquia 

aquifer, as inferred from contours of simulated steady-state 

water levels (fig. 27a), originated from western up-dip areas 

and flowed toward the east, where it bifurcated around the 

low-permeability hydrogeologic units filling the Chesapeake 

Bay impact crater. The steady-state flow to the northeast 

and southeast around the Chesapeake Bay impact crater 

discharged through adjacent hydrogeologic units, primarily 

the Nanjemoy-Marlboro confining unit. By 2003, groundwater 

withdrawals had produced a broad area of lower water levels 

west of the Chesapeake Bay impact crater (fig. 27b). The 

maximum water-level drawdown simulated in the Aquia 

aquifer just west of the Chesapeake Bay impact crater is about 

135 ft (fig. 27c).

Simulated groundwater withdrawals from the Aquia 

aquifer in Maryland have produced a well-defined drawdown 

bowl centered near longitude –76.4 ° in southern Maryland 

(fig. 27c). A groundwater divide near latitude 38 ° separates 

groundwater in the Aquia aquifer that flows from northern 

Virginia into Maryland from groundwater in the remainder of 

Virginia that flows generally toward the broad area of lower 

water levels west of the Chesapeake Bay impact crater. 

Water-level measurements from Aquia observation wells 

51M 18, 57G 2, and 55B 67 (fig. 28) are well-simulated in 

trend, although the water levels simulated for some observation wells are consistently greater than observed water-level 

altitudes by up to approximately 15 ft.

–80

–70

–60

–50

–40

–30

–20

–10

0

10

20

30

Simulated 56H 29

Observed 56H 29

Simulated 54Q 73

Observed 54Q 73

Simulated 56G 38

Observed 56G 38

1975 1980 1985 1990 1995 2000

YEAR

WATER-LEVEL ALTITUDE, IN FEET

Figure 26. Observed and simulated water levels in wells 54Q 73, 56G 38, and 56H 29 in the Piney Point aquifer

of the Virginia Coastal Plain. Well location is shown on figure18.
54 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 27a. Simulated predevelopment water levels in the Aquia aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

Figure 27a. Simulated predevelopment water levels in the Aquia aquifer of the Virginia Coastal Plain.

10

10

20

30

40

50

60

70

70

80

VIRGINIA

NORTH CAROLINA

40

50

60

70

80

30

40

50

20

30

20

0 to 20

21 to 40

41 to 60

61 to 80

81 to 100

101 to 120

121 to 145

WATER LEVEL, IN FEET

EXPLANATION

WATER LEVEL—Contour

interval 10 feet

10

CHESAPEAKE

BAY

ATLANTIC

OCEAN

37 ̊

38 ̊

77 ̊ 76 ̊
Simulated Conditions 55

Figure 27b. Simulated 2003 water levels in the Aquia aquifer of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

Figure 27b. Simulated 2003 water levels in the Aquia aquifer of the Virginia Coastal Plain.

–40

VIRGINIA

NORTH CAROLINA

30

40

50

20

60 70

80

10

100

0

–10

–20

–30

–40

–50

–60

–70

–80

90

–30

–10

0 50

40

30

10

10

37 ̊

38 ̊

77 ̊ 76 ̊

-

WATER LEVEL, IN FEET

EXPLANATION

WATER LEVEL—Contour

interval 10 feet

10

–90 to –60

–59 to –40

–39 to –20

–19 to 0

1 to 20

21 to 40

41 to 60

61 to 80

81 to 100

101 to 120

121 to 145

CHESAPEAKE

BAY

ATLANTIC

OCEAN
56 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 27c. Simulated drawdown from predevelopment to 2003 in the Aquia aquifer of the Virginia Coastal Plain.

DRAWDOWN, IN FEET

EXPLANATION

CHANGE IN WATER LEVEL—

 Contour interval 10 feet

10

–4 to 0

1 to 20

21 to 40

41 to 60

61 to 80

81 to 100

101 to 135

CHESAPEAKE

BAY

ATLANTIC

OCEAN

MARYLAND

Figure 27c. Simulated 2003 water-level drawdowns in the Aquia aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

60

VIRGINIA

NORTH CAROLINA

80

90

100

110

120

130

70

30

20

10

40

50

60

20

30

50

40

60

70 80

100 90

10

20

30

40

20

37 ̊

38 ̊

77 ̊ 76 ̊
Simulated Conditions 57

Figure 28. Observed and simulated water levels in wells 51M 18, 57G 2, and 55B 67 in the Aquia aquifer of the 

Virginia Coastal Plain. (Well locations are shown in figure 19.)

Peedee Aquifer

The Peedee aquifer is present in a limited area in 

southeast Virginia and extends to the south into North Carolina 

(figs. 29a, b). Contours of simulated water levels in the Peedee 

aquifer for the steady-state stress period (fig. 29a) indicate that 

predevelopment groundwater flowed from west to east through 

the Peedee aquifer. Lower water-level altitudes simulated 

for 2003 (fig. 29b) suggest that water levels have declined 

throughout the aquifer. The relation between simulated 

water levels and water levels measured in the Peedee aquifer 

observation well 61B 6 indicates that the model is accurately 

simulating actual water-level conditions in this aquifer 

(fig. 30). Because groundwater withdrawals from the Peedee 

aquifer have not been reported and unreported domestic 

withdrawals are considered not to occur (Pope and others, 

2007), no groundwater withdrawals from the Peedee aquifer 

were specified in the groundwater model. The simulated 

water-level decline in the Peedee aquifer, therefore, likely 

results from flow to hydraulically connected aquifers in 

North Carolina.

Simulated 51M 18

Observed 51M 18

Simulated 57G 2

Observed 57G 2

Simulated 55B 67

Observed 55B 67

1965 1970 1975 1980 1985 1990 1995 2000

YEAR

WATER-LEVEL ALTITUDE, IN FEET

–90

–80

–70

–60

–50

–40

–30

–20

–10

0

10

20

30

40

50

60

70

80

90

100

110

120

Figure 28. Observed and simulated water levels in wells 51M 18, 57G 2, and 55B 67 in the Aquia aquifer of the

Virginia Coastal Plain. Well location is shown on figure19.
58 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 29a. Simulated predevelopment water levels in the Peedee aquifer of the Virginia Coastal Plain.

MARYLAND

Figure 29a. Simulated predevelopment water levels in the Peedee aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

5

25

20

15

10

VIRGINIA

NORTH CAROLINA

WATER LEVEL, IN FEET

EXPLANATION

WATER LEVEL—Contour

interval 5 feet

5

1 to 5

6 to 10

11 to 15

16 to 20

21 to 25

26 to 29

CHESAPEAKE

BAY

ATLANTIC

OCEAN

37 ̊

38 ̊

77 ̊ 76 ̊
Simulated Conditions 59

Figure 29b. Simulated 2003 water levels in the Peedee aquifer of the Virginia Coastal Plain.

MARYLAND

Figure 29b. Simulated 2003 water levels in the Peedee aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

5

10

0

5

VIRGINIA

NORTH CAROLINA

0

–4 to 0

1 to 5

6 to 10

11 to 12

WATER LEVEL, IN FEET

EXPLANATION

WATER LEVEL—Contour

interval 5 feet

5

CHESAPEAKE

BAY

ATLANTIC

OCEAN

37 ̊

38 ̊

77 ̊ 76 ̊
60 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 30. Observed and simulated water levels in well 61B 6 in the Peedee aquifer of the Virginia Coastal Plain. 

(Well location is shown in figure 19.)

Virginia Beach Aquifer

The Virginia Beach aquifer is present in southeast 

Virginia, where it extends offshore beneath the continental 

shelf and into North Carolina (figs. 31a, b). Contours of the 

simulated water levels in the Virginia Beach aquifer for the 

steady-state stress period (fig. 31a) indicate that predevelopment groundwater flowed from western up-dip areas through 

the Virginia Beach aquifer toward the east. Reported and 

self-supplied withdrawals from the Virginia Beach aquifer in 

southeast Virginia caused simulated water levels to decline 

by the end of 2003, resulting in a relatively flat simulated 

potentiometric surface in the down-dip area of the aquifer 

(fig. 31b). Because water-level measurements from observation wells screened in the Virginia Beach aquifer were not 

available for model calibration or evaluation, the accuracy 

of these simulated water levels is unknown. A superposition 

analysis indicates that simulated domestic withdrawals may 

account for approximately 1 ft of drawdown in the Virginia 

Beach aquifer.

1975 1980 1985 1990 1995 2000

YEAR

WATER-LEVEL ALTITUDE, IN FEET

0

10

20

30

28

26

24

22

18

16

14

12

8

6

4

2

Simulated 61B 6

Observed 61B 6

Figure 30. Observed and simulated water levels in well 61B 6 in the Peedee aquifer of the Virginia Coastal Plain.

Well location is shown on figure18.
Simulated Conditions 61

Figure 31a. Simulated predevelopment water levels in the Virginia Beach aquifer of the Virginia Coastal Plain.

MARYLAND

Figure 31a. Simulated predevelopment water levels in the Virginia Beach aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

25

20

10

5

25

15

30

30

35

VIRGINIA

NORTH CAROLINA

WATER LEVEL, IN FEET

EXPLANATION

WATER LEVEL—Contour

interval 5 feet

5

1 to 10

11 to 20

21 to 30

31 to 40

41 to 50

51 to 60

CHESAPEAKE

BAY

ATLANTIC

OCEAN 37 ̊

38 ̊

77 ̊ 76 ̊
62 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 31b. Simulated 2003 water levels in the Virginia Beach aquifer of the Virginia Coastal Plain.

-

0

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

Figure 31b. Simulated 2003 water levels in the Virginia Beach aquifer of the Virginia Coastal Plain.

–10

–20

0

10

–10

0

VIRGINIA

NORTH CAROLINA

20

10

30

0

37 ̊

38 ̊

77 ̊ 76 ̊

WATER LEVEL, IN FEET

EXPLANATION

WATER LEVEL—Contour

interval 10 feet

10

–71 to –20

–19 to –10

–9 to 0

1 to 10

11 to 20

21 to 57

CHESAPEAKE

BAY

ATLANTIC

OCEAN
Simulated Conditions 63

Potomac Aquifer

The top of the Potomac aquifer is at a shallow depth near 

the Fall Line, where it has good hydraulic connection through 

the surficial aquifer to the water table and probably receives 

groundwater recharge (fig. 4a–d). The contours of simulated 

predevelopment water levels at the top of the Potomac aquifer 

(fig. 32a) suggest that groundwater-flow directions were 

influenced by the proximity of sea level along the James River 

near the Fall Line. Water-level contours east of the Fall Line 

in the northern and southern Coastal Plain indicate generally 

easterly-directed groundwater flow, which bifurcated to flow 

northeasterly and southeasterly around the Chesapeake Bay 

impact crater. Between these areas, east of Richmond near 

latitude 37.3°N, predevelopment groundwater flow from the 

north and south converged and discharged to the James River. 

Immediately northeast, within a broad area of relatively small 

water-level gradient, a groundwater divide near 37.4°N, 

–76.6°W separated this area of convergent flow from the area 

that flowed toward the north side of the Chesapeake Bay 

impact crater. Concentric contours near 37°N, –76.85°W 

(fig. 32a) delineate where the Potomac aquifer is near the 

land surface and proximal to the Nottoway River, which 

consequently largely controls water levels simulated at the top 

of the Potomac aquifer in that area.

Industrial and municipal groundwater withdrawals, 

which began to increase during the 1940s, have resulted in 

substantial decline of water levels in the Potomac aquifer and 

the formation of two major “cones of depression” in Virginia 

centered near Franklin and West Point (fig. 32b). At the end 

of the simulation in December 2003, the maximum drawdown 

simulated at the top of the Potomac aquifer was 235 ft near 

Franklin and 232 ft near West Point (fig. 32c). Regional 

groundwater-flow directions within the Potomac aquifer have 

changed substantially from their predevelopment pattern as 

a result of these declines. Groundwater within the Potomac 

aquifer generally flows toward one of the major withdrawal 

centers at either Franklin or West Point, although it may be 

intercepted by one of many Potomac wells elsewhere in the 

Coastal Plain. The water-level gradients near the Atlantic coast 

have reversed from their predevelopment direction, which 

has increased the potential for landward transport of saltwater 

toward production wells. As of 2006, consistent trends of 

increased chloride concentration have not been observed in 

wells, however.

Substantial groundwater withdrawals from the Potapsco 

aquifer in Maryland, which correlates by age, depth, and 

lithology with the Potomac aquifer in Virginia, have generated 

an area of substantial drawdown north of the model domain 

near longitude –77° in Maryland. Because these withdrawals 

cannot be explicitly simulated within the model domain, 

their effect on groundwater flow out of northern Virginia was 

simulated with a specified-flux boundary, which was discussed 

in the preceding section entitled “Transient Underflow to 

Maryland.” This flux has caused simulated groundwater 

flow in northern Virginia to be redirected under the Potomac 

River and north toward the Potapsco-withdrawal area in 

southern Maryland. 

Observed water levels from 37 Potomac aquifer observation wells (fig. 33) are shown with their simulated equivalents 

(figs. 34–39) at various locations throughout the Virginia 

Coastal Plain. In general, simulated water levels agree within 

several feet of observed water levels in the Potomac aquifer, 

with few differences in excess of 10 ft. One exception is 

observation well 58F 1 (fig. 37), located in James City County 

southeast of Williamsburg (fig. 33). Four observations beginning in 1998 are more than 10 ft lower than their simulated 

equivalents, and the observed water level in year 2002 is 41 ft 

lower than its simulated equivalent. These discrepancies may 

be because of unreported groundwater withdrawals that were 

not simulated in the model. Simulated equivalents to water 

levels in well 52J 34 (fig. 36), located in eastern Henrico 

County, are increasingly higher than the rising water levels 

observed after 1997. The observed rising water levels are 

likely because of decreased groundwater withdrawals in this 

area. Although a relatively low horizontal hydraulic conductivity (on the order of 10 ft/d) in this area yields a good overall 

fit to observed water levels, lithologies (such as a channel-fill 

facies) with higher horizontal hydraulic conductivity may be 

present locally and influence the observation-well response.
64 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 32a. Simulated predevelopment water levels in the 

Potomac aquifer of the Virginia Coastal Plain.

Franklin

MARYLAND

Figure 32a. Simulated predevelopment water levels in the Potomac aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

10

0

20

30

40

10

20

30

40

10

80

70

60

50

90 100

130

120

110

140

140

130

120

110

40

30

10

20

80

60

50

70 70

80

VIRGINIA

NORTH CAROLINA

80

70

37 ̊

38 ̊

77 ̊ 76 ̊

0

-

90

100

WATER LEVEL, IN FEET

EXPLANATION

WATER LEVEL—Contour

 interval 10 feet

10

–12 to 0

1 to 20

21 to 40

41 to 60

61 to 80

81 to 100

101 to 120

121 to 140

141 to 160

161 to 167

CHESAPEAKE

BAY

ATLANTIC

OCEAN

West Point
Simulated Conditions 65

Figure 32b. Simulated 2003 water levels in the Potomac aquifer 

of the Virginia Coastal Plain.

MARYLAND

NORTH CAROLINA

-

-

-

-

-

-

0

Figure 32b. Simulated 2003 water levels in the Potomac aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

–50

–60

–80

–90

–100

–110

–100

–120

0

10

30

50 70

90 120 100

–40

–30

–10

–20

0

–20 –20

–10

–20

30

20

10

40 50

60

VIRGINIA

WATER LEVEL, IN FEET

EXPLANATION

WATER LEVEL—Contour

 interval 10 feet

10

–185 to –160

–159 to –140

–139 to –120

–119 to –100

–99 to –80

–79 to –60

–59 to –40

–39 to –20

–19 to 0

1 to 20

21 to 40

41 to 60

61 to 80

81 to 100

101 to 120

121 to 140

141 to 160

161 to 180 

CHESAPEAKE

BAY

ATLANTIC

OCEAN

West Point

Franklin

–30

–10

–20

–30

–40

–50

–60

–70

–80

–90

–50

–60

–70

–40

–30 –20

–10

–20

37 ̊

38 ̊

77 ̊ 76 ̊
66 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 32c. Simulated drawdown from predevelopment to 2003 

in the Potomac aquifer of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

Figure 32c. Simulated 2003 water-level drawdowns in the Potomac aquifer of the Virginia Coastal Plain.

VIRGINIA

NORTH CAROLINA

80

70

60

50

40

30

20

10

70

80

90

100

110

120

130

140

60

50

30

40

10

20

90

100

110

120

80

70

60

50

40

30

130

140

10

20

30

40

50

160

150

170

180

DRAWDOWN, IN FEET

CHANGE IN WATER LEVEL—

 Contour interval 10 feet

10

–2 to 20

21 to 40

41 to 60

61 to 80

81 to 100

101 to 120

121 to 140

141 to 160

161 to 180

181 to 200

201 to 240

EXPLANATION

CHESAPEAKE

BAY ATLANTIC

OCEAN

West Point

Franklin

160

150

170

180

37 ̊

38 ̊

77 ̊ 76 ̊
Simulated Conditions 67

Figure 33. Locations of observation wells screened in the Potomac aquifer for hydrographs depicted in figures 34–39.

0

0

25 50 MILES

25 50 KILOMETERS

Delmarva Peninsula

VIRGINIA

NORTH CAROLINA

CHESAPEAKE

BAY

ATLANTIC

OCEAN

54R 2

58H 6

56J 11

54L 5

55L 1

59K 1

56M 6

54G 11 56G 8

52G 29

52J 34

55H 1 55H 20

52K 3

52F 1

61C 1

59C 30

60C 41

60B 3

58A 77

56A 1

55B 16 54B 12

55B 51

52A 1

54C 1 55C 3

53D 6

53E 5

57D 1

57E 10

58E 3

58F 1

56F 1

57G 1

55B 45

Figure 33. Locations of selected observation wells screened in the Potomac aquifer.

MARYLAND

37 ̊

38 ̊

77 ̊ 76 ̊
68 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 34.

Observed and 

simulated water 

levels of selected 

wells screened 

in the Potomac 

aquifer in the 

southwest model 

area. (Locations of 

wells are shown in 

figure 33.)

Figure 35.

Observed and 

simulated water 

levels of selected 

wells screened 

in the Potomac 

aquifer in the 

southeast model 

area. (Locations of 

wells are shown in 

figure 33.)

–190

–180

–170

–160

–150

–140

–130

–120

–110

–100

–90

–80

–70

–60

–50

–40

–30

–20

–10

0

10

20

Simulated 52A 1

Observed 52A 1

Simulated 54B 12

Observed 54B 12

Simulated 55B 16

Observed 55B 16

Simulated 55B 51

Observed 55B 51

Simulated 55B 45

Observed 55B 45

Simulated 55C 3

Observed 55C 3

1960 1965 1970 1975 1980 1985 1990 1995 2000

YEAR

WATER-LEVEL ALTITUDE, IN FEET

Figure 34. Observed and simulated water levels of selected wells screened in the Potomac aquifer in the

southwest model area. (Locations of wells are shown on figure 32.)

–110

–100

–90

–80

–70

–60

–50

–40

–30

–20

–10

0

Simulated 56A 1

Observed 56A 1

Simulated 58A 77

Observed 58A 77

Simulated 60B 3

Observed 60B 3

Simulated 59C 30

Observed 59C 30

Simulated 60C 41

Observed 60C 41

Simulated 61C 1

Observed 61C 1

1965 1970 1975 1980 1985 1990 1995 2000

YEAR

WATER-LEVEL ALTITUDE, IN FEET

Figure35. Observed and simulated water levels of selected wells screened in the Potomac aquifer in the

southeast model area. (Locations of wells are shown on figure 32.)
Simulated Conditions 69 Figure 36. Observed and simulated water levels of selected wells screened in the Potomac aquifer in the western model area. (Locations of wells are shown in figure 33.)

–10

–20–30–40–50

0

70

605040302010

Simulated 54C 1Observed 54C 1

Simulated 53D 6

Observed 53D 6

Simulated 53E 5

Observed 53E 5Simulated 52F 1Observed 52F 1

Simulated 52G 29

Observed 52G 29

Simulated 52J 34

Observed 52J 34

1905 1910 1915 1920 1925 1930 1935 1940 1945 1950 1955 1975 1980 1985 1960 1965 1970 1995 2000 1990

YEAR

WATER-LEVEL ALTITUDE, IN FEET

Figure 36. Observed and simulated water levels of selected wells screened in the Potomac aquifer in the western model area.(Locations of wells are shown on figure 32.)
70 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 37. Observed and simulated water levels of selected wells screened in the Potomac aquifer in the eastern model area. (Locations of wells are shown in figure 33.)

–100–110–120–130

–10–20–30–40–50–60–70–80–90

0

2010

Simulated 57D 1

Observed 57D 1Simulated 57E 10Observed 57E 10Simulated 58E 3Observed 58E 3

Simulated 58F 1

Observed 58F 1

Simulated 56F 1

Observed 56F 1Simulated 57G 1

Observed 57G 1

1945 1955 1960 1950 1965 1970 1990 1985 1980 1975 2000 1995

YEAR

WATER-LEVEL ALTITUDE, IN FEET

Figure 37. Observed and simulated water levels of selected wells screened in the Potomac aquifer in the eastern model area.(Locations of wells are shown on figure 32.)
Simulated Conditions 71 Figure 38. Observed and simulated water levels of selected wells screened in the Potomac aquifer in the central model area. (Locations of wells are shown in figure 33.)

–100–110–120

–10–20–30

–40

–50–60–70–80–90

0

302010

Simulated 56G 8

Observed 56G 8

Simulated 55H 1

Observed 55H 1

Simulated 55H 20

Observed 55H 20Simulated 54G 11Observed 54G 11

Simulated 56J 11

Observed 56J 11Simulated 58H 6Observed 58H 6

1950 1960 1965 1955 1970 1990 1995 1985 1980 1975 2000

YEAR

WATER-LEVEL ALTITUDE, IN FEET

Figure 38. Observed and simulated water levels of selected wells screened in the Potomac aquifer in the central model area.(Locations of wells are shown on figure 32.)
72 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 39. Observed and simulated water levels of selected wells screened in the Potomac aquifer in the northern model area. (Locations ofwells are shown in figure 33.)

1965 1975 1970 1980 1995 1990 1985 2000

YEAR

WATER-LEVEL ALTITUDE, IN FEET

–10–20–30–40–50–60–70–80

0

302010

Simulated 52K 3

Observed 52K 3Simulated 54L 5Observed 54L 5Simulated 55L 1

Observed 55L 1

Simulated 59K 1

Observed 59K 1Simulated 56M 6

Observed 56M 6

Simulated 54R 2

Observed 54R 2

Figure 39. Observed and simulated water levels of selected wells screened in the Potomac aquifer in the northern model area.(Locations of wells are shown on figure 32.)
Simulated Conditions 73

Water Budget

Recharge and evapotranspiration are large compared to 

other steady-state budget components and the pumpage and 

storage components that appear in the transient stress periods. 

Predevelopment water-budget components for the model 

domain (table 7) are discussed in the following section on the 

steady-state water budget. Water-budget components for the 

portion of the model in Virginia (table 8) were calculated 

with the program ZONEBUDGET (Harbaugh, 1990) for the 

steady-state and two transient stress periods. 

Table 7. Steady-state water-budget components for model domain.

Component

Volume 

(cubic feet 

per day)

Volume 

(million

gallons 

per day)

Percent 

of total

Inflows

Recharge to water table 1.47 x109 11,013 97.17

Simulated rivers above sea level 3.85 x107 288 2.54

Tidal rivers 2.30 x106 17 .15

Underflow from Piedmont 2.16 x106 16 .14

TOTALS 1.51 x109 11,334 100.00

Outflows

Evapotranspiration and non-simulated streams 1.34 x109 10,048 88.66

Simulated rivers above sea level 1.27 x108 953 8.40

Tidal rivers, Chesapeake Bay, ocean 4.45 x107 333 2.94

TOTALS 1.51 x109 11,334 100.00

Table 8. Steady-state and transient water-budget components for Virginia sub-domain.

Component

Predevelopment

volume 

(million gallons 

per day) 

2001

volume 

(million

gallons 

per day)

2003

volume 

(million

gallons 

per day)

Inflows

Recharge to water table 7,636 5,617 10,770

Net from storage 0 817 –1,157

Underflow from Piedmont 16 16 16

Underflow from North Carolina 2 23 23

TOTALS 7,654 6,473 9,652

Outflows

Evapotranspiration and non-simulated streams 6,891 5,725 8,650

Net leakage to rivers 492 389 617

Net leakage to tidewaters and bay 269 236 266

Underflow to Maryland 2 3 3

Reported pumpage 0 94 89

Domestic pumpage 0 26 27

TOTALS 7,654 6,473 9,652

Steady-State Water Budget

Components of the groundwater budget for simulated 

predevelopment (steady-state) conditions are summarized in 

table 7. For the initial (steady-state) model stress period, the 

specified annual recharge flux to the water table (22.95 in.) 

is equivalent to 50 percent of the average Coastal Plain 

precipitation. The 50 percent of precipitation that does not 

enter the groundwater model is considered overland runoff 

and evapotranspiration from the vadose zone. Relatively 

small quantities of recharge are simulated as underflow from 

the Piedmont Province to the west 

and seepage from rivers. Of the 

simulated recharged groundwater, 

a substantial fraction is lost as 

evapotranspiration from the water 

table and seepage to small streams, 

both of which are accounted for by 

the outflow flux simulated with the 

EVT package. As discussed in the 

section on rivers, only 4.2 percent 

of the total length of mapped rivers 

in the Virginia Coastal Plain is 

explicitly represented with the RIV 

package. Because seepage to the 

remaining rivers, which exist in 

every finite-difference cell, is not 

explicitly simulated, some of this 

seepage flux is effectively compensated for by increased simulated 

flux through the EVT package. 

Although the relative magnitude of 

these two components simulated 

with the EVT package is uncertain, 

approximately 35 percent of the 

simulated evapotranspiration flux 

(about 7.1 in.) probably occurs as 

seepage to the small streams shown 

in figure 7. Under that assumption, 

75 percent of the total base flow 

(9.8 in.) discharges to non-simulated 

small streams. In reality, these 

streams drain into the larger rivers, 

the seepage to which is explicitly 

simulated in the model. Although 

their combined length is 24 times 

that of the major rivers explicitly 

represented in the model, this flux 

is only 3.7 times greater than the 

explicitly simulated river flux. The 

flux per unit length to non-simulated 

streams is 15 percent of the flux per 

unit length to explicitly represented 

rivers, which is consistent with their 

smaller stream widths and associated seepage areas. 
74 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Transient Water Budget

As the magnitude of simulated recharge increases or 

decreases between stress periods during the transient simulation, the simulated evapotranspiration changes in response 

to the different head in the uppermost active model layer. 

Variations in recharge and ET result in smaller changes in 

storage and flow through head-dependent boundaries, which 

also respond to changes in specified pumpage.

To illustrate budget differences between the steady-state 

and transient stress periods, budget components from the 

Virginia model sub-domain have been summarized in table 8. 

Although both inflows and outflows occur from storage and 

head-dependent boundaries, the net volume for these components was computed in order to simplify table 8. The simulated recharge rate during the steady-state stress period was 

specified equal to 50 percent of the average historical Coastal 

Plain precipitation (fig. 11). In contrast, recharge specified 

for transient stress periods representing 2001 and 2003 were 

below and above the annual average, respectively. Extensive 

groundwater development resulted in lower water levels 

in southern Virginia by 2001, which also caused increased 

groundwater underflow from North Carolina toward the areas 

of drawdown in southern Virginia. Similarly, groundwater 

development in Maryland had increased the underflow of 

groundwater from Virginia into Maryland. (Flow leaving the 

model domain was not included as simulated underflow from 

Virginia to Maryland in table 8.) In comparison to the predevelopment conditions, by 2001 groundwater pumpage had 

caused withdrawals from storage and decreased net outflows 

to rivers and the Chesapeake Bay.

Following several years of drought conditions, 2003 was 

an unusually “wet” year in Virginia. Specified recharge to the 

surficial aquifer increased water levels in the surficial aquifer, 

and water-level declines in underlying aquifers generally 

decreased as well. The resulting addition to storage in the 

groundwater system is represented by the negative value in the 

“from storage” component of the 2003 Virginia water budget 

(table 8). In comparison to the predevelopment and 2001 

rates, the relative abundance of simulated water in the shallow 

aquifer system in 2003 resulted in more evapotranspiration 

and outflows to rivers, and outflow rates to tidewater bodies 

(including the Chesapeake Bay) close to predevelopment rates.

Calibration Assessment

If the weights of the water-level observations accurately 

represent observation error and the model correctly represents 

the relevant hydrologic processes influencing the simulated 

equivalents to the observations, the distribution of weighted 

residuals should be random in space and time. Parameter 

values should not be strongly correlated if they are to be 

independently estimated. The final estimated parameter values 

should be within the range of expected reasonable values, or 

their confidence intervals should overlap with the range of 

expected reasonable values. 

During the course of model calibration, many configurations of the model were tested with various parameterization 

schemes. For each of these model realizations, the optimal 

parameter values were estimated by minimization of the objective function with either PEST or UCODE–2005. Following 

each parameter-estimation run, the magnitude and distribution 

of residuals (observed minus simulated water levels) and 

estimated parameter values were evaluated to identify possible 

model error or bias. Model bias is a form of model error that 

may appear in the spatial or temporal pattern of residuals. 

Areas with patterns or magnitudes of residuals suggestive of 

model bias were identified for alternative parameter schemes 

or other possible approaches of improving model fit. The 

criteria used for these evaluations are reviewed in the following sections as they apply to the calibrated model.

Residual Plots

The model fit of simulated to observed water levels 

over the entire range of groundwater levels measured in the 

Virginia Coastal Plain is illustrated in figure 40a. The outlying 

group with observed water-level altitudes between 80 and 

110 ft, which are greater than their simulated equivalents, are 

from three Potomac aquifer observation wells. Two of these 

observation wells (52P 10 and 52N 1), which are 

close to mapped faults in Caroline County, are anomalous 

compared with other nearby Potomac observation wells, 

which have simulated and observed water levels in close 

agreement. Water levels observed in the third well (52H 17), 

located about 1.5 mi northeast of the James River in eastern 

Henrico County, are anomalously high and may represent 

perched conditions. These observations were assigned a large 

variance (low weight), so as not to influence aquifer-system 

parameter estimates.

The relation between weighted simulated and weighted 

observed water levels (fig. 40b) is a preferred indicator of 

model error because the expected error variance has been 

accounted for by the observation weights. Weighted residuals 

are calculated as the difference between weighted observed 

and weighted simulated water levels. Trends resulting from 

model bias are more easily detected when weighted residuals 

versus weighted simulated water levels are graphed as in 

figure 41. Ideally, the residual distribution should appear randomly distributed about the horizontal line with ordinate = 0. 

The weighted residuals in figure 41 do not appear entirely 

random. A pattern that appears non-random occurs for positive 

weighted residuals with weighted simulated equivalents near 

zero. These residuals correspond to observations in cells 

where the water level is controlled by head-dependent flux 

boundaries representing tidal river areas or the Chesapeake 

Bay. Several similar vertical clusters of weighted residuals 

along abscissa values greater than zero are likely because of 

the control of river boundaries on the simulated water levels 
Calibration Assessment 75 Figure 40a. Relation between observed and simulated water levels in the Virginia Coastal Plain. SIMULATED WATER-LEVEL ALTITUDE, IN FEET 120 100 80 60 40 20 0 –20 –40 –60 –80 –100 –120 –140 –160 –180

OBSERVED WATER-LEVEL ALTITUDE, IN FEET

130140120110100

908070605040302010

0

–10–20–30–40–50–60–70–80–90

–100–110–120–130–140–150–160–170–180–190

Figure 40a. Relation between observed and simulated water levels in the Virginia Coastal Plain.
76 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 40b. Relation between weighted observed water levels and weighted simulated water levels in the Virginia Coastal Plain.

WEIGHTED SIMULATED WATER-LEVEL ALTITUDE, IN FEET

100 90 80 70 60 50 40 30 20 10 0 –10 –20 –30 –40 –50 –60 –70 –80 –90

WEIGHTED OBSERVED WATER-LEVEL ALTITUDE, IN FEET

50

100105

959085807570656055504540353025201510

–5

–10–15–20–25–30–35–40–45–50–55–60–65–70–75–80–85–90–95

–100

Figure 40b. Relation between weighted observed water levels and weighted simulated water levels in the Virginia Coastal Plain.
Calibration Assessment 77

Figure 41. Relation between weighted residuals and simulated water levels in the Virginia Coastal Plain.

in cells containing observation wells. In these cases, residuals 

result from topographically influenced observed water levels 

not being accurately simulated at the regional model scale.

Weighted Residual Maps

If the weights of the water-level observations adequately 

represent observation error and the model sufficiently 

represents the relevant hydrologic processes influencing 

the simulated equivalents to the observations, the residual 

distribution should be random in space. Ideally, positive and 

negative weighted residuals should be approximately equal in 

number and appear randomly distributed. In a transient model, 

many water-level observations from each observation well 

may be used for model calibration. The observation with the 

maximum weighted residual from each observation well was 

selected to generate maps of weighted residuals. To depict 

the spatial distribution of model error within each aquifer 

(fig. 42a–f), the weighted residuals were sorted according to 

the aquifer screened by each corresponding observation well.

Because the residual (observed minus simulated water 

level) is weighted by the inverse of the error variance 

associated with each measurement, the magnitudes of 

weighted residuals differ from the raw residual magnitudes. 

For an observation with a large measurement-error variance 

(representing large uncertainty in well altitude or water-level 

measurement), the corresponding small regression weight 

may result in a weighted residual less than the difference 

between the measured water level and its simulated equivalent. Conversely, an accurate water-level measurement in a 

surveyed well can have a weighted residual greater than the 

difference between the measured water level and its simulated 

equivalent. For example, well 52G 29 on the map of Potomac 

weighted residuals (fig. 42f) appears to be among the largest 

positive weighted residuals, yet the actual differences between 

measured and simulated water levels in this well are only 

about 4 ft (fig. 36). The average (root mean square error) of 

the 5,087 weighted residuals for Virginia water-level observations is 3.6 ft.

Positive maximum weighted residuals outnumber 

negative maximum weighted residuals in the surficial aquifer 

(fig. 42a). This suggests that either (1) model error affects 

simulated water-table altitudes, or (2) a bias in the sampling 

of water-table altitudes by surficial-aquifer observation wells. 

Simulated water-table altitudes are calculated for the center 

of each finite-difference cell and are interpolated between 

the centers of the finite-difference cells to calculate the 

simulated equivalent to the water-table observation. Because 

of the strong influence of simulated evapotranspiration on the 

simulated water-table altitude within each finite-difference 

cell, the specified evapotranspiration-surface altitude for 

each cell largely controls the simulated water-table altitude. 

The evapotranspiration-surface altitude for each model cell 

was specified to be equal to the median altitude within each 

finite-difference cell determined from digital elevation model 

Figure 41. Relation between weighted residuals and simulated water levels in the

Virginia Coastal Plain.

SIMULATED WATER-LEVEL ALTITUDE, IN FEET

–175 –150 –125 –100 –75 –50 –25 0 25 50 75 100 125

WEIGHTED RESIDUAL

25

20

15

10

5

0

–5

–10

–15

–20

–25

–30
78 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 42a. Maximum residual magnitudes for observation wells screened in the surficial aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

Delmarva Peninsula

VIRGINIA

NORTH CAROLINA

CHESAPEAKE

BAY

ATLANTIC

OCEAN "

!!

"

!

!

!

!

!

!

!

!

! !

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

! ! !

!

"

"

"

"

""

"

"

"

"

"

"

"

"

"

!

!

!

!

Figure 42a. Maximum residual magnitudes for observation wells screened in the surficial aquifer of the

Virginia Coastal Plain.

–26.1 to –11.0

–10.9 to –7.0

–6.9 to –3.0

–2.9 to 0.0

0.1 to 3.0

3.1 to 7.0

7.1 to 11.0

11.1 to 15.8

WEIGHTED

RESIDUAL

EXPLANATION

MARYLAND

37 ̊

38 ̊

77 ̊ 76 ̊
Calibration Assessment 79

Figure 42b. Maximum residual magnitudes for observation wells screened in the Yorktown-Eastover aquifer of the 

Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

Delmarva Peninsula

VIRGINIA

NORTH CAROLINA

CHESAPEAKE

BAY

ATLANTIC

OCEAN "!"!

!

!

"

"

"

"

"

"

"

"

!!

!

!

!

!

!

!

!

! !

!

!

!

!

!

!

!

"

""

"

"

"

" "

"

"

"

"

"

"

"

"

"

!

!

!

!

Figure 42b. Maximum residual magnitudes for observation wells screened in the Yorktown-Eastover aquifer

of the Virginia Coastal Plain.

WEIGHTED

RESIDUAL

EXPLANATION

–20.6 to –11.0

–10.9 to –7.0

–6.9 to –3.0

–2.9 to 0.0

0.1 to 3.0

3.1 to 7.0

7.1 to 11.0

11.1 to 12.7

MARYLAND

37 ̊

38 ̊

77 ̊ 76 ̊
80 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 42c. Maximum residual magnitudes for observation wells screened in the Piney Point aquifer of the Virginia 

Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

Delmarva Peninsula

VIRGINIA

NORTH CAROLINA

CHESAPEAKE

BAY

ATLANTIC

OCEAN

37 ̊

38 ̊

77 ̊ 76 ̊

"

"

" "

"

"

"

"

"""

!

!

!

"

"

""""

"

""

!

!

!

!

"

"

"

"

"

!

!

!

!

"

Figure 42c. Maximum residual magnitudes for observation wells screened in the Piney Point aquifer of the

Virginia Coastal Plain.

WEIGHTED

RESIDUAL

EXPLANATION

–14.3 to –11.0

–10.9 to –7.0

–6.9 to –3.0

–2.9 to 0.0

1.3 to 3.0

3.1 to 7.0

7.1 to 9.2
Calibration Assessment 81

Figure 42d. Maximum residual magnitudes for observation wells screened in the Aquia aquifer of the Virginia 

Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

Delmarva Peninsula

VIRGINIA

NORTH CAROLINA

CHESAPEAKE

BAY

ATLANTIC

OCEAN

"

"

"

"

!

!

!

!

!

!

!

!

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

Figure 42d. Maximum residual magnitudes for observation wells screened in the Aquia aquifer of the

Virginia Coastal Plain.

WEIGHTED

RESIDUAL

EXPLANATION

–13.5 to –11.0

–10.9 to –7.0

–6.9 to –3.0

–2.9 to 0.0

0.1 to 3.0

3.1 to 7.0

7.1 to 10.2

MARYLAND

37 ̊

38 ̊

77 ̊ 76 ̊
82 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 42e. Maximum residual magnitudes for observation wells screened in the Peedee and St. Marys aquifers of 

the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

Delmarva Peninsula

VIRGINIA

NORTH CAROLINA

CHESAPEAKE

BAY

ATLANTIC

OCEAN

!

!

Figure 42e. Maximum residual magnitudes for observation wells screened in the Peedee and St. Marys

aquifers of the Virginia Coastal Plain.

6

9

!

!

H

WEIGHTED RESIDUAL—

ST. MARYS

2

WEIGHTED RESIDUAL—

PEEDEE

!

H

EXPLANATION

MARYLAND

37 ̊

38 ̊

77 ̊ 76 ̊
Calibration Assessment 83

Figure 42f. Maximum residual magnitudes for observation wells screened in the Potomac aquifer of the Virginia 

Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

Delmarva Peninsula

VIRGINIA

NORTH CAROLINA

CHESAPEAKE

BAY

ATLANTIC

OCEAN

Figure 42f. Maximum residual magnitudes for observation wells screened in the Potomac aquifer of the

Virginia Coastal Plain.

"

"

"

"

!

!

!

!

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

" " "

"

! !

!

!

!

! !

!

! ! ! !

! !

!

! !

!

!

!!

!

!

!!

! !

!

!

!

!

!

!

!

!

!

!

!

"

""

"

"

" "

"

"

"

"

"

"

"

"

!

!

!

!

"

"

"

" "

"

" "

"

" "

"

"

"

"

!

!

"

"

"

"

"

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

! !

!

!

!

!

!

! ! !

!

!

"

"

"

"

"

"

"

"

"

"

"

"

"

" "

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

"

" "

"

"

"

WEIGHTED

RESIDUAL

EXPLANATION

–20.7 to –11.0

–10.9 to –7.0

–6.9 to –3.0

–2.9 to 0.0

0.1 to 3.0

3.1 to 7.0

7.1 to 11.0

11.1 to 20.2

37 ̊

38 ̊

77 ̊ 76 ̊
84 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

(DEM) data. Analysis of DEM data for the Virginia Coastal 

Plain indicates that the range of altitudes within each model 

cell is of a similar order to the range of residual magnitudes. 

Comparison of observation well-head altitudes with their 

corresponding median cell altitudes indicates that 67 percent 

of observation wells used for model calibration have wellhead altitudes greater than the median altitude in the cell, 29 

percent have well-head altitudes less than the median altitude 

in the cell, and 4 percent have well-head altitudes equal to 

the median altitude in the cell. In general, observation-well 

altitudes are greater than the median land-surface altitudes 

calculated for the 1-mi2 model grid cells, indicating that 

wells tend to be located on higher-altitude areas. Because 

water-table altitudes correlate with topographic altitude, there 

is probably a positive bias in the sampling of water-table 

altitudes by surficial-aquifer observation wells. If such a bias 

exists, a larger number of positive residuals than negative 

residuals would be expected for the surficial aquifer. Although 

unequal in number, the positive and negative residuals appear 

to have random spatial distribution, suggesting the residualsign inequality more likely results from observation-sampling 

bias than model error.

The spatial distribution of positive and negative maximum weighted residuals within the Yorktown-Eastover aquifer 

(fig. 42b) appears to be random, indicating it is not affected by 

substantial systematic model error.

A clustering of negative maximum weighted residuals 

within the Piney Point aquifer near Williamsburg (fig. 42c) 

indicates that the spatial distribution of positive and negative 

maximum weighted residuals is not random and that simulated 

water levels are affected by some model error. As discussed in 

the preceding section on simulated conditions within the Piney 

Point aquifer, the cause of this model error is uncertain but 

may be either insufficient simulated domestic-withdrawal flux 

or inaccurate simulated confining-unit leakance in the area of 

the clustered negative residuals.

A prevalence of negative weighted residuals within 

the Aquia aquifer indicates that the spatial distribution of 

positive and negative maximum weighted residuals is not 

random (fig. 42d). Hydrographs of simulated and observed 

Aquia water levels show that simulated water levels in some 

observation wells are greater than observed water levels by 

up to approximately 15 ft. A complex distribution of Aquia 

horizontal hydraulic conductivity with 5 parameters associated 

with 20 pilot points did not substantially improve the fit of 

simulated to observed water levels during model calibration 

using regression. It is possible that localized inaccuracy 

of spatial leakance through either the Nanjemoy-Marlboro 

confining unit or Potomac confining zone may contribute 

to the over-prediction of water levels in some of the Aquia 

observation wells, but that hypothesis was not tested during 

model calibration.

The number of observation wells was insufficient to 

evaluate the spatial distribution of residuals for the Peedee, 

Virginia Beach, or St. Marys aquifers. Residuals from observation wells in the Peedee and St. Marys aquifers are depicted in 

figure 42e to show observation-well locations.

Frequent analysis of the spatial distribution of Potomac 

aquifer weighted residuals during model calibration was 

necessary to evaluate different configurations of pilot points or 

zones, each with their associated parameters, which were used 

to specify horizontal variations of hydraulic conductivity. The 

spatial clustering of positive and negative maximum weighted 

residuals within the Potomac aquifer (fig. 42f) indicates 

that this residual distribution for the calibrated model is not 

entirely random. 

Parameter Correlation

The parameter correlation coefficients computed by 

UCODE–2005 are useful for designing parameterization 

schemes, as well as determining which parameters may be 

effectively estimated by the parameter-estimation regression 

process. Examination of the parameter correlation coefficients following each regression during model calibration 

occasionally suggested modifications to the parameterization 

scheme. Several parameter pairs in the final model version 

have large correlation coefficients (table 9). The “recharge” 

parameter, which represents the fraction of precipitation that 

infiltrates to the water table, is completely correlated with 

the “et_max_rate” parameter, which represents the maximum 

evapotranspiration flux from the water table. This correlation 

indicates that it is not possible to estimate the magnitude of 

either parameter by regression. Recharge was also somewhat 

correlated with the parameter representing the evapotranspiration extinction depth. The regression was sensitive to the ratio 

of “recharge” and “et_max_rate”; final values were obtained 

for “et_max_rate” by regression after fixing the “recharge” 

parameter and the evapotranspiration extinction depth at 

reasonable values.

The parameter “kdepth1,” which is used to simulate the 

change with depth of Potomac aquifer horizontal hydraulic 

conductivity, is correlated with the parameter “k_islewight,” 

Table 9. Parameter pairs with absolute correlation 

coefficients greater than 0.8.

[et_max_rate, maximum evapotranspiration flux; recharge, specified 

recharge flux; k_islewight, horizontal hydraulic conductivity at pilot 

point in Isle of Wight County; kdepth1,  in exponent in equation 8; 

k_potom_cz, horizontal hydraulic conductivity of Potomac confining 

zone; k_megabloc, horizontal hydraulic conductivity of Potomac aquifer 

within Chesapeake Bay impact crater; k_kinggeorg, horizontal hydraulic 

conductivity at pilot point in King George County; mdflow, specified flux 

out of model row 1] 

Parameter pair Correlation coefficient

 et_max_rate recharge 1.00

 k_islewight kdepth1 0.95

 k_potom_cz k_megabloc –0.94

 k_kinggeorg mdflow 0.86
Calibration Assessment 85

representing Potomac aquifer horizontal hydraulic conductivity at a pilot point in Isle of Wight County. 

The parameter “k_kinggeorg,” representing Potomac 

aquifer horizontal hydraulic conductivity at a pilot point in 

model row 1 and in King George County, is correlated with 

the parameter “mdflow,” which represents the magnitude 

of groundwater flux out of the model domain in this area 

toward Maryland.

The parameter “k_potom_cz,” which represents horizontal hydraulic conductivity within the Potomac confining zone, 

is inversely correlated with the parameter “k_megabloc,” 

which represents horizontal hydraulic conductivity of the 

Potomac aquifer within the Chesapeake Bay impact crater. 

(The composite scaled sensitivities of “k_megabloc” indicates 

it is not a sensitive parameter; this apparent inverse correlation may be a numerical artifact.) Regression runs with 

UCODE–2005 and PEST were made with sub-sets of the total 

model-parameter set. Possible difficulties estimating unique 

values for parameters in pairs with correlation coefficients 

between 0.8 and 0.95 (considered “somewhat correlated” 

parameters) were avoided by including members of the 

correlated pairs in separate parameter-estimation runs.

Parameter Sensitivities

The composite scaled sensitivity reflects the total amount 

of information provided by the water-level observations for 

the estimation of one model parameter (Hill and Tiedeman, 

2007). The composite scaled sensitivities for estimated parameters were computed by UCODE–2005 and are displayed 

graphically in figure 43. In general, parameters with composite 

scaled sensitivities greater than about 1 can be well-estimated 

and have correspondingly smaller 95-percent confidence 

intervals. Parameters with composite scaled sensitivities 

smaller than about 1 may be more difficult to estimate and(or) 

have larger 95-percent confidence intervals. 

Estimated Parameters

The model-parameter values estimated by non-linear 

regression are shown with their respective reasonable-value 

ranges and individual linear 95-percent confidence intervals 

in figure 44. The confidence intervals were computed by 

assuming that the model is linear, true observation errors are 

Figure 43. Composite scaled sensitivities for model parameters. (Parameter definitions are on p. x.) Figure 43. Composite scaled sensitivities for model parameters. k_islewight kdepth1 ss_potomac vk_poto_cz ss_upper_ce vk_nanjemo k_franklin k_sohampt_n vk_smary_c ss_other

piedflow

k_newportne

k_aquia

recharge

k_jamescity

vk_upper_c

vk_surface

k_chesapeak

vk_potomac

mdflo

w

evapotransp

vk_calve_c

k_aquitard

k_princegeo

k_chesterfi

k_pineyp_jc

river_cond

K_SUSSEX

k_york_e_a

k_sohampton

k_potom_cz

k_henrico

k_westmore

k_charlesci

k_kinggeorg

k_columbia

k_virgin_a

vk_chickah

k_kingwilli

ss_aquitard

k_lancaster

k_westwestm

spec_yield

k_peedee_a

k_maryland

k_s_mary_a

vk_james_c

vk_virgb_c

vk_york_cz

rw

k_megabloc

vk_peede_c

k_aquiclud

vk_aquifer

COMPOSITE SCALED SENSITIVITY

4.6

4.8

4.4

4.2

4.0

3.8

3.6

3.4

3.2

3.0

2.8

2.6

2.4

2.2

2.0

1.8

1.6

1.4

1.2

1.0

0.8

0.6

0.4

0.2

0
86 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 44. Estimated values, confidence intervals, and reasonable ranges for parameters representing horizontal hydraulic conductivity, 

vertical hydraulic conductivity, specific storage, and other fluxes. (Parameter definitions are on p. x.)

Figure 44. Estimated values, confidence intervals, and reasonable ranges for parameters representing horizontal hydraulic

conductivity, vertical hydraulic conductivity, specific storage, and other fluxes.

PARAMETER

k_columbia

k_york_e_a

k_s_mary_a

k_pineyp_jc

k_aquia

k_peedee_a

k_virgin_a

k_kinggeorg

k_maryland

k_westwestm

k_chesterfi

k_lancaster

k_westmore

k_jamescity

k_newportne

k_sussex

k_islewight

k_sohamptonk_sohampt_n

k_princegeo

k_kingwilli

k_franklin

k_chesapeak

k_potom_cz

k_aquitard

vk_surface

vk_york_cz

vk_smary_c

vk_calve_c

vk_james_c

vk_chickah

vk_nanjemo

vk_upper_c

vk_poto_cz

vk_potomac

ss_other

ss_aquitard

ss_upper_ce

ss_potomac

spec_yield

mdflow

piedflow

kdepth1

et_max_rate

recharge

PARAMETER VALUE

1x10–7

1x10–6

1x10–5

1x10–4

1x10–3

1x10–2

1x10–1

1x100

1x10+1

1x10+2

1x10+3

1x10+4

Range of 95-percent linear

confidence interval 

Estimated value

Expected range of reasonable

values
Model Limitations 87

random, and weighted true errors are independent. Because the 

groundwater model is non-linear, the depicted intervals are not 

entirely accurate. Calculation of the more accurate non-linear 

confidence intervals (Vecchia and Cooley, 1987) requires 

months of computer execution time, which was prohibitive 

to this project. Despite their inherent inaccuracy, the depicted 

confidence intervals are somewhat useful indicators of the 

relative confidence of the estimated parameter values.

All estimated parameter values fall within their respective 

ranges of reasonable values for hydraulic properties or model 

parameters of their type. Horizontal hydraulic conductivity 

(Kh, designated by the prefix ‘k’ in fig. 44) values are larger 

than values for vertical hydraulic conductivity (prefix: ‘vk’), 

specific storage (prefix: ‘ss’), and other parameters; parameters 

representing Kh in hydrogeologic units and at different pilotpoint locations are shown on the left side of figure 44. With 

the exception of parameters representing horizontal hydraulic 

conductivity in the St. Marys and Peedee aquifers, where few 

observations exist, and at a pilot-point location in Westmoreland county (representing Kh

 in the Potomac aquifer), the 

computed 95-percent confidence intervals for these parameters 

are also within a range of reasonable values for horizontal 

hydraulic conductivity of these hydrogeologic units.

The magnitudes of parameter values for specific yield, 

Piedmont and Maryland specified fluxes, evapotranspiration 

and recharge rates, and the decrease of horizontal hydraulic 

conductivity with depth in the Potomac aquifer, while not 

meaningful compared to each other, are also shown with 

respect to their confidence intervals and reasonable ranges on 

the right side of figure 44. With the exception of the parameter 

representing simulated well radii, which is insensitive, limits 

of the 95-percent confidence intervals for each parameter are 

within an order of magnitude of their corresponding estimated 

parameter values.

Parameters representing the vertical hydraulic conductivities of aquifers and the Peedee and Virginia Beach confining 

units have both very small composite scaled sensitivities and 

large linear confidence intervals and were not well-estimated 

by the model-calibration regression. The parameters representing horizontal hydraulic conductivity in the Chickahominy and 

Exmore confining units and the horizontal hydraulic conductivity in the structurally disrupted “mega block” zone of the 

Potomac aquifer similarly had very small composite scaled 

sensitivities and large linear confidence intervals and also were 

not well-estimated by the model-calibration regression. 

Aquifer and Confining-Unit Spatial Variability

Because the thickness (b) of every aquifer is spatially 

variable, the transmissivity (T = Khb) and storage coefficient 

(Ssb) of every aquifer varies spatially. Similarly, the thickness 

of each confining unit and confining zone is spatially variable; 

consequently, the leakance (Kv/b) of these hydrogeologic 

units also varies spatially. Maps of model-calculated transmissivity for each aquifer are presented in figures 45a–f. Maps 

of model-calculated leakance for each confining unit are 

presented in figures 46a–k. Maps of aquifer transmissivity 

may be useful for visualizing differences in the relative 

potential productivity of the thinner confined aquifers above 

the Potomac aquifer system. These maps are also useful for 

visualizing the extent and variations in thickness of each 

hydrogeologic unit. 

Model Limitations

The model was designed to simulate regional 

groundwater flow and water levels in the confined aquifers 

of the Virginia Coastal Plain. For most of the study area, 

water-table altitudes and horizontal flow directions within the 

surficial aquifer are affected by topographic variations and 

streams that cannot be adequately simulated with the 1-mi2

finite-difference-cell size used for the model. Therefore, the 

model may not be suitable, for example, to study site-specific 

groundwater contamination problems in the surficial aquifer. 

The groundwater model discretization, which employs 

1-mi2 finite-difference cells, is most appropriate for evaluating 

the drawdown of withdrawal wells at distances greater than 

approximately 1 mi from those wells. Because specified withdrawals are simulated from the center of each finite-difference 

cell, the simulated withdrawal may be horizontally displaced 

up to 3,733 ft from the actual location of the withdrawal 

well. A Cooper-Jacob (1946) approximation of water-level 

drawdown resulting from a withdrawal of 0.1 Mgal/d 

(3,000,000 gal per month) from an aquifer with values of 

transmissivity T = 10,000 ft2/d and storage S = 0.0001, which 

are typical of those estimated from Potomac aquifer tests 

(table 1), is 0.7 ft at a distance of 1 mi after 1 year. The similar 

drawdown approximation for the withdrawal of 0.01 Mgal/d 

(300,000 gal per month) is less than 0.1 ft at 1 mi after 1 year. 

Drawdown resulting from well withdrawals between 0.01 and 

0.1 Mgal/d may be of concern within a mile or two from a 

pumping well. Because the finite-difference approximation is 

less accurate at this scale, it may be more effective to employ 

a local-scale model or analytical solution when evaluating 

drawdown at this range. 

Because the minimum stress-period length of the model 

is 1 year, seasonal variations of groundwater levels are not 

simulated. However, it would be possible to simulate future 

seasonal variations by modifying the temporal discretization 

and specifying appropriate fluxes in the model stress packages.

Although the change of groundwater density resulting 

from spatial variation of solute concentration is simulated, 

the model is not currently capable of accurately simulating or 

predicting temporal changes of solute concentration resulting 

from saltwater intrusion. A useful transport model of areas of 

potential seawater intrusion will likely require finer horizontal 

discretization and a computationally intensive solution 

algorithm, such as Total-Variation-Diminishing (Zheng and 

Wang, 1999) that minimizes the effects of numerical dispersion and oscillation on simulated solute concentrations.
88 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 45a. Estimated transmissivity of the Yorktown-Eastover aquifer of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

TRANSMISSIVITY, IN

FEET SQUARED PER DAY

EXPLANATION

Figure 45a. Estimated transmissivity of the Yorktown-Eastover aquifer of the Virginia Coastal Plain.

37 ̊

38 ̊

77 ̊ 76 ̊

20 to 1,000

1,001 to 2,000

2,001 to 3,000

3,001 to 5,000

5,001 to 10,360

VIRGINIA

NORTH CAROLINA

CHESAPEAKE

BAY

ATLANTIC

OCEAN
Model Limitations 89

Figure 45b. Estimated transmissivity of the St. Marys aquifer of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

29 to 1,000

1,001 to 2,000

2,001 to 3,000

3,001 to 5,000

5,001 to 7,200

Figure 45b. Estimated transmissivity of the St. Marys aquifer of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

TRANSMISSIVITY, IN

FEET SQUARED PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
90 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 45c. Estimated transmissivity of the Piney Point aquifer of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

23 to 1,000

1,001 to 2,000

2,001 to 3,000

3,001 to 5,000

5,001 to 8,262

Figure 45c. Estimated transmissivity of the Piney Point aquifer of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

TRANSMISSIVITY, IN

FEET SQUARED PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
Model Limitations 91

Figure 45d. Estimated transmissivity of the Aquia aquifer of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

78 to 2,000

2,001 to 4,000

4,001 to 8,000

8,001 to 12,000

12,001 to 16,852

Figure 45d. Estimated transmissivity of the Aquia aquifer of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

TRANSMISSIVITY, IN

FEET SQUARED PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
92 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 45e. Estimated transmissivity of the Peedee aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

489 to 1,000

1,001 to 2,000

2,001 to 3,000

3,001 to 5,000

5,001 to 10,275

Figure 45e. Estimated transmissivity of the Peedee aquifer of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

TRANSMISSIVITY, IN

FEET SQUARED PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

MARYLAND

37 ̊

38 ̊

77 ̊ 76 ̊
Model Limitations 93

Figure 45f. Estimated transmissivity of the Virginia Beach aquifer of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

10 to 150

151 to 300

301 to 600

601 to 1,000

1,001 to 1,640

Figure 45f. Estimated transmissivity of the Virginia Beach aquifer of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

TRANSMISSIVITY, IN

FEET SQUARED PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
94 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 46a. Estimated leakance of the Yorktown confining zone of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

4.8x10–4 to 1.0x10–3

1.1x10–3 to 2.0x10–3

2.1x10–3 to 3.0x10–3

3.1x10–3 to 1.0x10–2

1.1x10–2 to 4.4x10–2

Figure 46a. Estimated leakance of the Yorktown confining zone of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

LEAKANCE PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
Model Limitations 95

Figure 46b. Estimated leakance of the St. Marys confining unit of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

1.5x10–7 to 5.0x10–7

5.1x10–7 to 1.0x10–6

1.1x10–6 to 2.0x10–6

2.1x10–6 to 4.0x10–6

4.1x10–6 to 7.3x10–5

Figure 46b. Estimated leakance of the St. Marys confining unit of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

LEAKANCE PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

MARYLAND

37 ̊

38 ̊

77 ̊ 76 ̊
96 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 46c. Estimated leakance of the Calvert confining unit of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

2.4x10–8 to 5.0x10–7

5.1x10–7 to 1.0x10–6

1.1x10–6 to 2.0x10–6

2.1x10–6 to 4.0x10–6

4.1x10–6 to 3.4x10–5

Figure 46c. Estimated leakance of the Calvert confining unit of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

LEAKANCE PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
Model Limitations 97

Figure 46d. Estimated leakance of the Chickahominy confining unit of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

7.6x10–8 to 5.0x10–7

5.1x10–7 to 1.0x10–6

1.1x10–6 to 2.0x10–6

2.1x10–6 to 4.0x10–6

4.1x10–6 to 5.7x10–6

Figure 46d. Estimated leakance of the Chickahominy confining unit of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

LEAKANCE PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

MARYLAND

37 ̊

38 ̊

77 ̊ 76 ̊
98 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 46e. Estimated leakance of the Exmore matrix confining unit of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

9.9x10–8 to 5.0x10–7

5.1x10–7 to 1.0x10–6

1.1x10–6 to 2.0x10–6

2.1x10–6 to 4.0x10–6

4.1x10–6 to 4.0x10–5

Figure 46e. Estimated leakance of the Exmore matrix confining unit of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

LEAKANCE PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
Model Limitations 99

Figure 46f. Estimated leakance of the Exmore clast confining unit of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

8.5x10–9 to 5.0x10–7

5.1x10–7 to 1.0x10–6

1.1x10–6 to 2.0x10–6

2.1x10–6 to 4.0x10–6

4.1x10–6 to 4.0x10–5

Figure 46f. Estimated leakance of the Exmore clast confining unit of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

LEAKANCE PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
100 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 46g. Estimated leakance of the Nanjemoy-Marlboro confining unit of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

2.8x10–8 to 5.0x10–7

5.1x10–7 to 1.0x10–6

1.1x10–6 to 2.0x10–6

2.1x10–6 to 4.0x10–6

4.1x10–6 to 3.0x10–5

Figure 46g. Estimated leakance of the Nanjemoy-Marlboro confining unit of the Virginia Coastal Plain.

CHESAPEAKE

BAY ATLANTIC

OCEAN

LEAKANCE PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

MARYLAND

37 ̊

38 ̊

77 ̊ 76 ̊
Model Limitations 101

Figure 46h. Estimated leakance of the Peedee confining zone of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

5.1x10–6

5.2x10–6 to 7.0x10–6

7.1x10–6 to 9.0x10–6

9.1x10–6 to 2.0x10–6

2.1x10–6 to 6.3x10–6

Figure 46h. Estimated leakance of the Peedee confining zone of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

LEAKANCE PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
102 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 46i. Estimated leakance of the Virginia Beach confining zone of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

8.5x10–9 to 4.0x10–8

4.1x10–8 to 6.0x10–8

6.1x10–8 to 7.0x10–8

7.1x10–8 to 1.0x10–7

1.1x10–7 to 1.0x10–6

Figure 46i. Estimated leakance of the Virginia Beach confining zone of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

LEAKANCE PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
Model Limitations 103

Figure 46j. Estimated leakance of the Upper Cenomanian confining unit of the Virginia Coastal Plain.

MARYLAND

0

0

25 50 MILES

25 50 KILOMETERS

7.5x10–8 to 5.0x10–7

5.1x10–7 to 1.0x10–6

1.1x10–6 to 2.0x10–6

2.1x10–6 to 4.0x10–6

4.1x10–6 to 1.5x10–5

Figure 46j. Estimated leakance of the Upper Coenomanian confining unit of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

LEAKANCE PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

37 ̊

38 ̊

77 ̊ 76 ̊
104 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure 46k. Estimated leakance of the Potomac confining zone of the Virginia Coastal Plain.

0

0

25 50 MILES

25 50 KILOMETERS

6.0x10–8 to 5.0x10–7

5.1x10–7 to 1.0x10–6

1.1x10–6 to 2.0x10–6

2.1x10–6 to 4.0x10–6

4.1x10–6 to 3.5x10–5

Figure 46k. Estimated leakance of the Potomac confining zone of the Virginia Coastal Plain.

CHESAPEAKE

BAY

ATLANTIC

OCEAN

LEAKANCE PER DAY

EXPLANATION

VIRGINIA

NORTH CAROLINA

MARYLAND

37 ̊

38 ̊

77 ̊ 76 ̊
Summary 105

Because the transient release of water stored in lowpermeability confining units is simulated, the drawdown 

resulting from simulated pumping stresses may change 

substantially through time before reaching steady state. Model 

testing indicates that several decades may be required for the 

aquifer system to equilibrate, or approximate steady-state 

conditions, after an appreciable change in pumping conditions. 

One procedure sometimes used for assessing future groundwater levels for management purposes is to include all pumping 

stresses in a steady-state simulation. Such an evaluation of the 

“eventual water levels” does not consider the timing of influence of hydraulic boundaries and may, therefore, over- and 

under-predict water levels at different locations through time. 

A more accurate evaluation of probable future aquifer-system 

responses to proposed pumping should be attainable from 

transient simulations of water levels at different future times.

Summary

Groundwater withdrawals have lowered water levels 

in Virginia Coastal Plain aquifers and have resulted in 

drawdown in the Potomac aquifer exceeding 200 ft in some 

areas. The discovery of the buried impact crater beneath the 

Chesapeake Bay is the most dramatic of several changes in the 

conceptualization of the Virginia Coastal Plain hydrogeologic 

framework that affects the simulation of groundwater flow. In 

order to issue groundwater withdrawal permits, State watermanagement authorities need an improved groundwater model 

capable of predicting aquifer-system responses to proposed 

groundwater withdrawals. In response to this need, a transient, 

three-dimensional variable-density groundwater flow model 

of the Virginia Coastal Plain aquifer system was developed 

and calibrated. 

The model simulates historical water levels and groundwater drawdown resulting from development since 1890. The 

113-year historical transient simulation was time-discretized 

into 34 stress periods of various lengths. The spatial scale 

of the model is appropriate for assessing the regional 

water-level responses of the confined aquifers of the Coastal 

Plain. Because groundwater flow in the unconfined surficial 

aquifer is influenced by topography and streams that cannot 

be represented at the regional model scale, local horizontal 

groundwater flow through the surficial aquifer is not intended 

to be accurately simulated. However, vertical flow through 

the surficial aquifer from and to the water table in areas of 

groundwater recharge and discharge is simulated and allows 

transmission of groundwater to and from underlying confined 

aquifers. Representation of recharge, evapotranspiration, and 

interaction with surface-water features, such as major rivers, 

lakes, the Chesapeake Bay, and the Atlantic Ocean, enable 

simulation of shallow flow-system details that influence 

centers of recharge to and discharge from the deeper confined 

flow system. The increased groundwater density associated 

with the transition from fresh to salty groundwater near the 

Atlantic Ocean affects regional groundwater flow and was 

simulated with the VDF Process of SEAWAT. The specified 

groundwater density distribution was generated by a separate 

108,000-year simulation of Pleistocene freshwater flushing 

around the Chesapeake Bay impact crater during transient 

sea-level changes, which is described in the appendix. 

Specified-flux boundaries simulate increasing groundwater underflow out of the model domain into Maryland and 

a small amount of recharge underflow from the Piedmont 

Physiographic Province over the Norfolk Arch. Historically 

reported municipal, industrial, commercial, and public-supply 

groundwater withdrawals accounted for approximately 

75 percent of total Coastal Plain groundwater withdrawals 

in the year 2000, the locations of which are known and were 

specified in the transient model. The remaining fraction of the 

total Coastal Plain groundwater withdrawals are unreported 

self-supplied domestic withdrawals, the locations of which 

are uncertain but were spatially estimated and allocated to 

individual aquifers. 

The groundwater flow model was calibrated to 7,183 

historic water-level observations from 497 observation 

wells with the parameter-estimation codes UCODE–2005 

and PEST. Most water-level observations were from the 

Potomac aquifer system, which consequently enabled a more 

detailed estimation of the spatial distribution of hydraulic 

parameters than was possible for other hydrogeologic units. 

The zone-boundary approach was used to distribute hydraulic 

parameters according to known hydrogeologic boundaries and 

geologic features, such as the extent of the Chesapeake Bay 

impact crater. Spatial interpolation of horizontal hydraulic 

conductivity within the Potomac aquifer and vertical hydraulic 

conductivity within the Calvert confining unit was facilitated 

with the use of pilot points. The average (RMSE) of 5,087 

weighted water-level residuals in Virginia is 3.6 ft, indicating 

a good fit of simulated to observed water levels and that the 

model is a good representation of the physical groundwater 

flow system. Magnitudes and distributions of water-level 

residuals indicate that, while not entirely absent, model bias is 

not problematic. 

The simulated water-table altitude in the Virginia 

Coastal Plain varies with a higher spatial frequency than the 

potentiometric surface of any underlying confined aquifer. The 

simulated evapotranspiration flux largely controls the simulated water-table altitude, resulting in a water-table surface 

that is a smoothed and subdued replica of the land-surface 

altitude represented in the model EVT package. It is unlikely 

that the surficial aquifer conducts shallow groundwater flow 

horizontally for distances substantially greater than 1 mi 

because of the numerous rivers or non-simulated streams 

available for discharge.

Simulated contours of water levels in the YorktownEastover aquifer indicate that predevelopment and current 

directions of groundwater flow are away from the higheraltitude areas in southern Virginia and near the midlines of the 

Northern Neck, Middle Peninsula, and York-James Peninsula 
106 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

toward discharge areas in lower altitude areas, including major 

rivers and the Chesapeake Bay.

Simulated steady-state water-level contours suggest that 

predevelopment groundwater flowed from west to east through 

the St. Marys aquifer. Water levels simulated in the St. Marys 

aquifer for 2003 are up to 40 ft lower than those simulated for 

predevelopment time, but suggest similar flow directions.

Simulated predevelopment groundwater flow in the 

Piney Point aquifer is from western up-dip areas toward 

the east and southeast, where it discharges through the 

overlying Calvert confining unit near the Atlantic coast. By 

2003, groundwater withdrawals had formed several cones of 

depression, with a simulated maximum drawdown of 127 ft. 

Groundwater flow within the Piney Point aquifer is directed 

toward these drawdown cones in the vicinity of the Middle 

and York-James Peninsulas.

The simulated steady-state groundwater flow within the 

Aquia aquifer originated from western up-dip areas and flowed 

toward the east, where flow paths split around low-permeability hydrogeologic units within the Chesapeake Bay impact crater. By 2003, groundwater withdrawals had produced a broad 

area of lower water levels west of the Chesapeake Bay impact 

crater, where the maximum groundwater level drawdown 

was 135 ft. Substantial simulated groundwater withdrawals 

from the Aquia aquifer in southern Maryland have produced a 

well-defined cone of depression there.

Contours of the simulated steady-state water levels in 

the Peedee aquifer indicate that predevelopment groundwater 

flowed from west to east through the aquifer. Observed water 

levels and the simulated potentiometric surface for 2003 

indicate that water levels have declined in the western area of 

the Peedee aquifer, possibly because of flow to hydraulically 

connected aquifers in North Carolina.

Contours of the simulated steady-state water levels in the 

Virginia Beach aquifer indicate that predevelopment groundwater flowed from western up-dip areas through the Virginia 

Beach aquifer toward the east. Reported and self-supplied 

withdrawals from the Virginia Beach aquifer in southeast 

Virginia have caused simulated water levels to decline, 

resulting in a relatively flat simulated potentiometric surface 

at the end of 2003. Simulated conditions in the Virginia Beach 

aquifer are poorly constrained because of the lack of observation wells.

Simulated steady-state water levels at the top of the 

Potomac aquifer indicate a complex pattern of predevelopment 

groundwater flow that was influenced by the proximity of 

sea level along the James River near the Fall Line. Simulated 

water-level contours east of the Fall Line in the northern and 

southern Coastal Plain indicate west to east paths of groundwater flow, which became northeasterly and southeasterly 

where flow bifurcated on either side of the Chesapeake Bay 

impact crater. Between these areas and east of Richmond, 

groundwater flow from the north and south converged to flow 

west and discharge to the James River.

Industrial and municipal groundwater withdrawals, 

which began to increase during the 1940s, have resulted in a 

substantial decline of water levels in the Potomac aquifer and 

in the formation of two major cones of depression in Virginia 

centered near Franklin and West Point. The maximum simulated drawdown near the top of the Potomac aquifer was 235 ft 

near Franklin and 232 ft near West Point by the year 2003. As 

a result of these declines, regional groundwater flow velocities 

within the Potomac aquifer have changed substantially from 

their predevelopment state and generally are directed toward 

either of the major withdrawal centers at Franklin or West 

Point. The water-level gradients near the Atlantic coast have 

reversed from their predevelopment direction, which has 

increased the potential for landward transport of salty water 

toward production wells.

Acknowledgments

Discussions with Robin Patton, Beverly Quinlan, 

T. Scott Bruce, and Terry D. Wagner of the Virginia DEQ, 

and Scott Emry and Whitney Katchmark of HRPDC provided 

useful guidance for this study. Additional discussions with 

and assistance from Leonard Orzol, Mary Hill, Eileen Poeter, 

Steffen Mehl, Richard Yager, Christian Langevin, Randall 

McFarland, Ward Sanford, and other USGS colleagues are 

greatly appreciated. A substantial portion of this study was 

supported by the USGS Office of Groundwater (OGW) 

through the Groundwater Resources Program. The authors 

thank the OGW for its support of the study.

References Cited

Anderman, E.R., and Hill, M.C., 2000, MODFLOW–2000, 

The U.S. Geological Survey modular ground-water 

model—Documentation of the Hydrogeologic-Unit Flow 

(HUF) package: U.S. Geological Survey Open-File Report 

00–342, 89 p.

Anderman, E.R., and Hill, M.C., 2003, MODFLOW–2000, 

the U.S. Geological Survey modular ground-water model—

Three additions to the Hydrogeologic-Unit Flow (HUF) 

Package—Alternative storage for the uppermost active 

cells, flows in hydrogeologic units, and the hydraulicconductivity depth-dependence (KDEP) capability: U.S. 

Geological Survey Open-File Report 03–347, 36 p.

Canadell, J., Jackson, R.B., Ehleringer, J.R., Mooney, H.A., 

Sala, O.E., and Schulze, E.D., 1996, Maximum rooting 

depth of vegetation types at the global scale: Oecologia, 

v. 108, p. 583–595.

Cederstrom, D.J., 1943, Chloride in groundwater in the coastal 

plain in southeastern Virginia: Virginia Geological Survey 

Bulletin 58, 36 p. 4 pl.
References Cited 107

Cederstrom, D.J., 1945, Geology and ground-water resources 

of the Coastal Plain in southeastern Virginia: Virginia Geological Survey Bulletin 63, 384 p.

Cederstrom, D.J., 1957, Geology and ground-water resources 

of the York-James Peninsula, Virginia: U.S. Geological 

Survey Water-Supply Paper 1361, 237 p.

Cederstrom, D.J., 1969, Geology and ground-water resources 

of the Middle Peninsula, Virginia: U.S. Geological Survey 

Open-File Report 69–37, 231 p.

Cooper, H.H., and Jacob, C.E., 1946, A generalized graphical 

method for evaluating formation constants and summarizing well field history: American Geophysical Union Trans., 

v. 27, p. 526–534.

Daniel, C.C., III, 1989, Statistical analysis relating well yield 

to construction practices and siting of wells in the Piedmont 

and Blue Ridge Provinces of North Carolina: U.S. Geological Survey Water-Supply Paper 2341–A, 27 p.

Doherty, John, 2003, Ground water model calibration using 

pilot points and regularization: Ground Water, v. 41, no. 2, 

p. 170–177.

Doherty, John, 2004, PEST model-independent parameter 

estimation (5th ed.): Corinda, Australia, Watermark Numerical Computing, 336 p.

Domenico, P.A., and Mifflin, M.D., 1965, Water from low-permeability sediments and land subsidence: Water Resources 

Research, v. 1, no. 4, p. 563–576.

Domenico, P.A., and Schwartz, F.W., 1990, Physical and 

chemical hydrogeology: New York, John Wiley & Sons, 

824 p.

Epstein, V.J., 1987, Hydrologic and geologic factors affecting 

land subsidence near Eloy, Arizona: U.S. Geological Survey 

Water-Resources Investigations Report 87–4143, 28 p.

Gardner, W.R., and Fireman, M., 1958, Laboratory studies of 

evaporation from soil columns in the presence of a water 

table: Soil Science, v. 85, p. 244–249.

Guo, W., and Langevin, C.D., 2002, User’s guide to SEAWAT; 

a computer program for simulation of three-dimensional 

variable-density ground-water flow: U.S. Geological Survey 

Techniques of Water-Resources Investigations 6–A7, 77 p.

Halford, K.J., and Hanson, R.T., 2002, User guide to the 

drawdown-limited, multi-node well (MNW) package for 

the U.S. Geological Survey’s modular three-dimensional 

finite-difference ground-water flow model, versions MODFLOW–96 and MODFLOW–2000: U.S. Geological Survey 

Open-File Report 02–293, 33 p.

Harbaugh, A.W., 1990, A computer program for calculating 

subregional water budgets using results from the U.S. Geological Survey modular three-dimensional finite-difference 

ground-water flow model: U.S. Geological Survey OpenFile Report 90–392, 24 p.

Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, 

M.G., 2000, MODFLOW–2000, the U.S. Geological Survey modular ground-water model—User guide to modularization concepts and the ground-water flow process: U.S. 

Geological Survey Open-File Report 00–92, 121 p.

Harsh, J.F., and Laczniak, R.J., 1990, Conceptualization and 

analysis of ground-water flow system in the Coastal Plain of 

Virginia and adjacent parts of Maryland and North Carolina: 

U.S. Geological Survey Professional Paper 1404–F, 100 p.

Helm, D.C., 1974, Evaluation of stress-dependent aquitard 

parameters by simulation of observed compaction from 

known stress history: Berkeley, University of California, 

Ph.D. dissertation, 175 p.

Heywood, C.E., 1995, Investigation of aquifer-system 

compaction in the Hueco Basin, El Paso, Texas, USA, in

Proceedings of the Fifth International Symposium on Land 

Subsidence: The Hague, IAHS Publication no. 234, p. 

35–45.

Heywood, C.E., 1998, Piezometric—Extensometric estimations of specific storage in the Albuquerque Basin, New 

Mexico, in Borchers, J.W., ed., Land subsidence case studies and current research; Proceedings of the Dr. Joseph F. 

Poland symposium on land subsidence, Association of 

Engineering Geologists Special Publication No. 8, Belmont, 

Calif., Star Publishing Company, p. 435–440.

Hill, M.C., 1998, Methods and guidelines for effective model 

calibration: U.S. Geological Survey Water-Resources Investigations Report 98–4005, 90 p.

Hill, M.C., Banta, E.R., Harbaugh, A.W., and Anderman, E.R., 

2000, MODFLOW–2000, The U.S. Geological Survey 

modular ground-water model—User guide to the observation, sensitivity, and parameter-estimation processes and 

three post-processing programs: U.S. Geological Survey 

Open-File Report 00–0184, 209 p.

Hill, M.C., and Tiedeman, C.R., 2007, Effective groundwater model calibration—with analysis of data, sensitivities, 

predictions, and uncertainty: Hoboken, NJ, John Wiley & 

Sons, 455 p.

Ireland, R.L., Poland, J.F., and Riley, F.S., 1984, Land subsidence in the San Joaquin Valley, California, as of 1980: U.S. 

Geological Survey Professional Paper 437–I, 105 p., 1 pl.

Johnson, A.I., 1967, Specific yield—Compilation of specific 

yields for various materials: U.S. Geological Survey WaterSupply Paper 1662–D, 74 p.
108 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Kernodle, J.M., and Philip, R.D., 1988, Using a geographic 

information system to develop a ground-water flow model: 

American Water Resources Association Monograph, series 

no. 14, p. 191–202.

Langevin, C.D., Shoemaker, W.B., and Guo, W., 2003, MODFLOW–2000, The U.S. Geological Survey modular groundwater model—Documentation of the SEAWAT–2000 

version with the variable-density flow process (VDF), and 

the integrated MT3DMS transport process (IMT): U.S. 

Geological Survey Open-File Report 03–426, 43 p.

Larsen, C.E., and Clark, I., 2003, Scale in sea level studies 

in Book of Abstracts, Coastal Sediments ‘03, Crossing 

Disciplinary Boundaries: Fifth International Symposium 

on Coastal Engineering and Science of Coastal Sediment 

Processes, May 18–23, 2003, Clearwater, FL, p. 432–433.

Leake, S.A., and Claar, D.V., 1999, Procedures and computer 

programs for telescopic mesh refinement using MODFLOW: U.S. Geological Survey Open-File Report 99–238, 

53 p.

Leake, S.A., and Lilly, M.R., 1997, Documentation of a computer program (FHB1) for assignment of transient specifiedflow and specified-head boundaries in applications of the 

Modular Finite-Difference Ground-Water Flow Model 

(MODFLOW): U.S. Geological Survey Open-File Report 

97–571, 56 p.

Leahy, P.P., and Martin, M., 1993, Geohydrology and simulation of ground-water flow in the northern Atlantic Coastal 

Plain aquifer system: U.S. Geological Survey Professional 

Paper 1404–K, 81 p.

McDonald, M.G., and Harbaugh, A.W., 1988, A modular 

three-dimensional finite-difference ground-water flow 

model: U.S. Geological Survey Techniques of WaterResources Investigations, book 6, chap. A1, 586 p. 

McFarland, E.R., 1998, Design, revisions, and considerations 

for continued use of a ground-water-flow model of the 

Coastal Plain aquifer system in Virginia: U.S. Geological 

Survey Water-Resources Investigations Report 98–4085, 

49 p.

McFarland, E.R., and Bruce, T.S., 2005, Distribution, origin, 

and resource-management implications of ground-water 

salinity along the western margin of the Chesapeake Bay 

impact structure in eastern Virginia, in Horton, J.W., Jr., 

Powars, D.S., and Gohn, G.S., eds., Studies of the Chesapeake Bay impact structure—The USGS–NASA Langely 

corehole, Hampton, Virginia, and related coreholes and 

geophysical surveys: U.S. Geological Survey Professional 

Paper 1688, chap. K, 32 p.

McFarland, E.R., and Bruce, T.S., 2006, The Virginia Coastal 

Plain hydrogeologic framework: U.S. Geological Survey 

Professional Paper 1731, 118 p. 25 pls.

Meisler, H., 1989, The occurrence and geochemistry of salty 

ground water in the northern Atlantic Coastal Plain: U.S. 

Geological Survey Professional Paper 1404–D, 51 p.

Meisler, H., Leahy, P.P., and Knobel, L.L., 1984, Effect of 

eustatic sea-level changes on saltwater-freshwater in the 

Northern Atlantic Coastal Plain: U.S. Geological Survey 

Water-Supply Paper 2255, 28 p.

Meng, A.A., III, and Harsh, J.F., 1988, Hydrogeologic framework of the Virginia Coastal Plain: U.S. Geological Survey 

Professional Paper 1404–C, 82 p.

National Climate Data Center, 2002, accessed March 20, 

2009, at http://www.ncdc.noaa.gov/oa/ncdc.html

Nelms, D.L., Harlow, G.E., Plummer, N., and Busenberg, E., 

2003, Aquifer susceptibility in Virginia, 1998–2000: U.S. 

Geological Survey Water-Resources Investigations Report 

03–4278, 58 p.

Orzol, L.L., 1997, User’s guide for MODTOOLS—Computer programs for translating data of MODFLOW and 

MODPATH into geographic information system files: U.S. 

Geological Survey Open-File Report 97–240, 86 p.

Peltier, R.W., 1994, Ice age paleotopography: Science, v. 265, 

p. 195–201.

Poag, C.W., 1999, Chesapeake invader—Discovering America’s giant meteorite crater: Princeton, NJ, Princeton University Press, 183 p.

Poag, W.C., Powars, D.S., Poppe, L.J., and Mixon, R.B., 1994, 

Meteoroid mayhem in Ole Virginny—Source of the North 

American tektite strewn field: Geology, v. 22, p. 691–694.

Poeter, E.P., Hill, M.C., and Banta, E.R., 2005, UCODE_2005 

and six other computer codes for universal sensitivity analysis, calibration, and uncertainty evaluation: U.S. Geological 

Survey Techniques and Methods, book 6, chap. 11, 283 p.

Pollock, D.W., 1994, Users guide for MODPATH/MODPATHPLOT, version 3: A particle tracking post-processing package for MODFLOW, the U.S. Geological Survey finite-difference ground-water flow model: U.S. Geological Survey 

Open-File Report 94–464, 249 p.

Pope, J.P., and Burbey, T., 2004, Multiple-aquifer characterization from single borehole extensometer records: Groundwater, v. 42, no. 1, p. 45–58.

Pope, J.P., McFarland, E.R., and Banks, R.B., 2007, Private domestic-well characteristics and the distribution 

of domestic withdrawals among aquifers in the Virginia 

Coastal Plain: U.S. Geological Survey Scientific Investigations Report 2007–5250, 47 p., http://pubs.usgs.gov/

sir/2007/5250/
References Cited 109

Powars, D.S., 2000, The effects of the Chesapeake Bay impact 

crater on the geologic framework and the correlation of 

hydrogeologic units of Southeastern Virginia, south of the 

James River: U.S. Geological Survey Professional Paper 

1622, 53 p., 1 pl.

Powars, D.S., and Bruce, T.S., 1999, The effects of the Chesapeake Bay impact crater on the geological framework and 

correlation of hydrogeologic units of the lower York-James 

Peninsula, Virginia: U.S. Geological Survey Professional 

Paper 1612, 82 p., 9 pls.

Powell, J.D., and Abe, J.M., 1985, Availability and quality of ground water in the Piedmont province of Virginia: 

U.S. Geological Survey Water-Resources Investigations 

Report 85–4235, 32 p. 

Quinlan, B.J., 2004, Virginia Coastal Plain model 2002 withdrawals simulation, reported use: Virginia Department of 

Environmental Quality, Water Resources Division, Office of 

Ground-Water Withdrawal Permitting, 37 p.

Quinlan, B.J., 2005, Virginia Coastal Plain model 2004 withdrawals simulation: Virginia Department of Environmental 

Quality, Water Resources Division, Office of Ground-Water 

Withdrawal Permitting, 43 p.

Richardson, D.L., 1994a, Hydrogeology and analysis of the 

ground-water flow system of the Eastern Shore, Virginia: 

U.S. Geological Survey Water-Supply Paper 2401, 82 p., 

1 pl.

Richardson, D.L., 1994b, Ground-water discharge from the 

Coastal Plain of Virginia: U.S. Geological Survey WaterResources Investigation Report 93–4191, 15 p.

Riley, F.S., 1998, Mechanics of aquifer systems—The scientific legacy of Joseph F. Poland, in Borchers, J.W., ed., Land 

subsidence case studies and current research, Proceedings 

of the Dr. Joseph F. Poland symposium on land subsidence: 

Belmont, CA, Star Publishing Company, Association of 

Engineering Geologists Special Publication No. 8, p. 13–27.

Rutledge, A.T., 1993, Computer programs for describing the 

recession of ground-water discharge and for estimating 

mean ground-water recharge and discharge from streamflow 

records: U.S. Geological Survey Water-Resources Investigation Report 93–4121, 45 p.

Sanford, S., 1913, The underground water resources of the 

Coastal Plain province of Virginia: Virginia Geological 

Survey Bulletin No. 5, 361 p. 

Sanford, W.E., Pope, J.P., and Nelms, D.L., 2009, Simulation 

of ground-water level and salinity changes in the Eastern 

Shore, Virginia: U.S. Geological Survey Scientific Investigations Report 2009–5006, 125 p.

Thornthwaite, C.W., 1948, An approach toward a rational classification of climate: Geographic Review, v. 38, p. 55–94.

University of Virginia Climatology Office, 2007; accessed on 

March 6, 2008, at http://climate.virginia.edu/

U.S. Geological Survey, 2002, National Water-Use Information Program 2000 data for counties; accessed in May 2002 

at http://water.usgs.gov/watuse/

Vecchia A.V., and Cooley, R.L., 1987, Simultaneous confidence and prediction intervals for nonlinear regression models with application to a groundwater flow model: Water 

Resources Research, v. 22, no. 2, p. 95–108.

Wright, J., 2006, Virginia Coastal Plain model 2005 withdrawals simulation: Virginia Department of Environmental 

Quality, Water Resources Division, Office of Ground-Water 

Withdrawal Permitting, 50 p.

Zheng, C., and Wang, P., 1999, MT3DMS—A modular threedimensional multispecies transport model for simulation of 

advection, dispersion, and chemical reactions of contaminants in groundwater systems—Documentation and user’s 

guide: U.S. Army Corps of Engineers, U.S. Army Engineer 

Research and Development Center, Strategic Environmental 

Research and Development Program, 221 p.
110 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Appendix

Transient Saltwater-Transport Simulation

A three-dimensional digital representation of the Virginia 

Coastal Plain groundwater density distribution is needed for 

SEAWAT to realistically simulate regional groundwater flow 

near the saltwater transition zone. A preliminary simulation 

of Coastal Plain saltwater transport generated a saltwater 

distribution in reasonable agreement with observed concentration data and supported the “differential flushing” hypothesis 

for the persistence of saltwater within the Chesapeake Bay 

impact crater (Heywood, 2003). A possible alternative 

approach for generating a representative density distribution is 

extrapolation of observed concentration data, which is sparse 

in Virginia and preferentially samples low-concentration 

areas in relatively high-permeability hydrogeologic units. 

The modeling approach described in this appendix was 

considered superior because it facilitated generation of a 

three-dimensional density distribution consistent with the 

geometry of aquifers and confining units in the hydrogeologic 

framework. The model described in this section employs 

different boundary conditions than are employed for the model 

described in the main report. To avoid confusion, the 108,000-

year saltwater-transport model described in this appendix 

is termed the “glacial model,” whereas the term “historical 

model” refers to the 113-year transient model described in the 

main report body.

The objective of the glacial model was to generate a 

reasonable approximation of the present saltwater transition 

zone, across which groundwater density increases by about 

2.5 percent. The late Pleistocene through Holocene epochs 

are simulated in the glacial model, during which multiple ice 

ages caused sea-level transgressions and regressions across 

the continental shelf and current Coastal Plain of Virginia 

(fig. A1). This history of late Quaternary sea-level change 

has influenced the configuration of the resulting (present 

day) saltwater transition zone. During Pleistocene glacial 

periods, when sea level was as much as 300 ft below its 

current level, the average hydraulic gradient through Coastal 

Plain hydrogeologic units was three to four times greater than 

present, and local gradients near the present day Atlantic coast 

were substantially greater. These gradients enabled enhanced 

flushing of meteoric groundwater through relatively hydraulically conductive hydrogeologic units. Saltwater was retained 

in areas of low regional horizontal hydraulic conductivity, 

such as the crater-filling clast- and matrix-supported Exmore 

confining units (Exmore Tsunami Breccia). In the glacial 

simulation described in this appendix, recharged freshwater 

displaced and mixed with saltwater to form a freshwatersaltwater transition zone that migrates toward a position of 

hydrodynamic equilibrium.

A subsidiary use of the glacial model was for simulation 

of the fraction of Carbon-14 (14C) in the dissolved carbon in 

Coastal Plain groundwater. Groundwater age dates can be 

used as a surrogate for flow data, such as streamflow gains and 

losses, in constraining values of groundwater flux (and hence 

recharge and hydraulic conductivities) through the aquifer 

system. Simulations of 14C ages using particle tracking with 

the steady-state stress period of the historical Coastal Plain 

model were unsatisfactory, in part, because the simulated 

hydraulic gradients were not representative of the probable 

Coastal Plain hydraulic gradients during the Pleistocene ice 

ages, during which older-observed waters were transported 

from their recharge areas. Apparent groundwater ages calculated from Coastal Plain 14C observations generally increase 

from west to east and with depth and range from about 

9,000 to 40,000 radiocarbon years (Nelms and others, 2003). 

The time scale of the glacial model was sufficiently long to 

simulate the effect of sea-level changes on hydraulic gradients 

and consequent transport velocities. Dissolved 14C concentrations were simulated by specifying the 14C radioactive-decay 

rate in the MT3DMS (Zheng and Wang, 1999) Chemical 

Reaction Package in SEAWAT–2000. Because the constituent 

lithologies of Virginia Coastal Plain hydrogeologic units are 

predominantly siliciclastic, with negligible carbonate rocks, 

the effects of carbon-exchange reactions in the aquifers and 

confining units were not simulated. The simulated 14C concentration distribution (expressed as Percent-Modern-Carbon) is 

in general agreement with 14C observations (considering the 

probable range of effective porosity and other uncertainties) 

and served to validate the conceptual flow model and the 

simulated recharge (within an order of magnitude).

Figure A1. Sea levels during the Pleistocene and Holocene epochs (modified from Zellmer, 1979). Figure A1. Sea levels during the Pleistocene and Holocene (modified from Zellmer, 1979).

100

0

–400

–300

–200

–100

SEA LEVEL, IN FEET, 

REFERENCED TO NGVD 29

100

–400

–300

–200

–100

0 100 200 300 400 500 600 700 800 900

TIME, IN THOUSANDS OF YEARS BEFORE PRESENT

0
Appendix 111

Boundary Condition Modifications

The model domain, hydrogeologic framework, and 

associated hydraulic properties of the glacial model described 

in this appendix are identical to those of the historical 

model. The purpose of the glacial model required changing 

the top-surface boundary conditions of the standard model 

to simulate changes in sea level. The head-dependent flux 

boundaries at the top surface of the historical model used to 

simulate recharge, evapotranspiration, and seepage to and 

from rivers, the Chesapeake Bay, and the Atlantic Ocean floor 

were replaced by time-varying specified-head boundaries, 

which were simulated with the FHB package (Leake and 

Lilly, 1997). Specified-head boundaries on the continental 

shelf and slope represented lowered sea levels at various times 

during the late-Pleistocene simulation. The water level in 

the top-most-active finite-difference cell at each row-column 

location was specified as the greater of either the cell-surface 

altitude or the Pleistocene sea level. Cell-surface altitudes 

were referenced to NGVD 29 and are negative for cells under 

the Chesapeake Bay or Atlantic Ocean. The FHB package 

internally interpolated sea levels between the specified levels 

(table A1) at times selected to represent the Pleistocene 

sea-level history depicted in figure A1. 

The solute concentrations at specified-head boundaries 

were specified to simulate appropriate recharge concentrations. 

Recharge through areas with top-surface altitudes above the 

simulated sea level (table A1) was specified at the freshwater 

concentration; boundary concentrations at areas below 

the simulated sea level were specified at the seawater 

concentration. 

Initial Conditions

Solute initial conditions were specified homogeneously 

at seawater concentration. This simple initial condition 

represents the salinity distribution that might have existed following a marine transgression similar to that which occurred 

during the Pliocene epoch. The extent of saltwater inundation 

into Coastal Plain sediments resulting from the Pleistocene 

marine transgression during the Eemian-Sangamon interglacial 

era (130,000–115,000 years before present) is uncertain and 

may not have been as extensive as during the Pliocene epoch. 

Because simulated saltwater is transported from shallow and 

western Coastal Plain sediments after several thousand years, 

the specified homogeneous initial condition does not substantially alter the final simulated configuration of the saltwater 

transition zone.

Simulation of Mass Transport

Values of aquifer-system properties required for variabledensity flow and mass transport simulation are summarized 

in table A2. To facilitate model calibration to different 

salinity indicators, the groundwater-solute concentration C

was normalized to the concentration for a seawater density 

ρs. The groundwater density was calculated from the equation of state where ∂ ∂ =  C 1.56 lb/ft3 for the normalized 

solute concentration C. Numerical dispersion of simulated 

concentrations may generate unrealistic densities; such effects 

were controlled by limiting calculated groundwater densities 

to range between the values corresponding to freshwater and 

seawater concentrations (table A2). Neither chemical reactions 

nor sorption were simulated because the major solute species 

comprising seawater are not considered to react or adsorb on 

the porous media.

The governing equations for variable-density groundwater flow and mass transport (see Numerical Method section 

of main report) were explicitly coupled and iteratively solved 

(with a one time step lag) by the Variable-Density Flow 

(VDF) and Integrated Mass Transport (IMT) processes of 

SEAWAT–2000 (Langevin and others, 2003). 

The mass-transport equation was solved with implicit 

finite differences using upstream weighting and a Courant 

number of 1. Upstream weighting was computationally 

Table A1. Specified PleistoceneEpoch sea levels.

Simulation 

time

(years) 

Sea level

(feet 

NGVD 29) 

Time 

before 

present

(years)

 0 –30 108,000

 28,000 –30 80,000

 38,000 –270 70,000

 48,000 –130 60,000

 63,000 –300 45,000

 94,000 –300 14,000

101,000 –30 7,000

108,000 0 0

Table A2. Values used for transport simulation.

 [ft, feet; ft2/d, feet squared per day; lb/ft3, pound per cubic foot]

Symbol Parameter name Value

αL longitudinal dispersivity 100 ft

αTH transverse horizontal dispersivity 10 ft

αTV transverse vertical dispersivity 1 ft

Dm molecular diffusion coefficient 10–3 ft2/d

ρf freshwater density 62.43 lb/ft3

ρs seawater density 63.99 lb/ft3

 effective porosity .2

C0 initial solute concentration 1

Cr recharge concentration 0
112 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

efficient (simulations ran in 1 to 2 days) but introduced 

numerical dispersion. Other solution techniques less 

susceptible to numerical dispersion were tested but required 

unacceptably long computer-execution times.

Figure A2. Simulated water density near the saltwater transition zone of the Virginia Coastal Plain.

The simulated groundwater densities in the vicinity of the 

saltwater transition zone are depicted in plan and cross-section 

views in figures A2 and A3, respectively. Although they 

were not used in the formal estimation of model parameters 

MARYLAND

A A’

A A’

Figure A2. Simulated water density near the present-day saltwater transition zone of the Virginia Coastal Plain.

63.99

63.60

63.21

62.82

62.43

DENSITY, IN POUNDS 

PER CUBIC FEET

EXPLANATION

0

0

25 50 MILES

25 50 KILOMETERS

Line of section

Location of chlorideobservation well

37 ̊

76 ̊ 75 ̊

38 ̊

CHESAPEAKE

BAY

ATLANTIC OCEAN

VIRGINIA

NORTH CAROLINA
Appendix 113

Figure A3. Simulated water density near the saltwater transition zone of the Virginia Coastal Plain. (Location of cross section shown in figure A2.)

63.99

63.6063.2162.8262.43

A A’

DENSITY, IN POUNDS 

EXPLANATION

PER CUBIC FEET
114 Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia

Figure A4. Relation between simulated and observed normalized chloride 

concentrations.

described in the main report body, 117 chloride-concentration 

observations helped to constrain the modern salinity distribution calculated for the end of the 108,000-year simulation 

described in this appendix. The relation between observed 

normalized chloride concentrations and simulated equivalents 

is shown in figure A4. In general, the simulated concentrations 

agree reasonably well with observed concentrations near the 

transition zone in deeper parts of the aquifer system. However, 

shallow wells with observed chloride concentrations up to 

about 13 percent that of seawater are simulated as essentially 

fresh and appear along the x-axis in figure A4. Because 

configuration of the simulated saltwater transition zone within 

the deeper aquifers, particularly the Potomac, was considered 

the more important influence on regional groundwater flow, 

the less accurate fit to observed salinities in other areas is 

acceptable for the purpose of this simulation.

Figure A4. Relation between simulated and observed normalized chloride concentrations.

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

OBSERVED NORMALIZED CHLORIDE CONCENTRATION

SIMULATED NORMALIZED CHLORIDE CONCENTRATION

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3
Appendix 115

Appendix References Cited

Heywood, C.E., 2003, Influence of the Chesapeake Bay 

impact structure on ground-water flow and salinity in the 

Atlantic Coastal Plain aquifer system of Virginia, in MODFLOW and More 2003—Understanding through modeling: Proceedings of the 5th International Conference of the 

International Ground Water Modeling Center, Golden, CO, 

Colorado School of Mines, p. 871–875.

Langevin, C.D., Shoemaker, W.B., and Guo, W., 2003, 

MODFLOW–2000, The U.S. Geological Survey 

modular ground-water model—Documentation of the 

SEAWAT–2000 version with the variable-density flow 

process (VDF), and the integrated MT3DMS transport 

process (IMT): U.S. Geological Survey Open-File 

Report 03–426, 43 p.

Leake, S.A., and Lilly, M.R., 1997, Documentation of a 

computer program (FHB1) for assignment of transient 

specified-flow and specified-head boundaries in applications of the Modular Finite-Difference Ground-Water Flow 

Model (MODFLOW): U.S. Geological Survey Open-File 

Report 97–571, 56 p.

Nelms, D.L., Harlow, G.E., Plummer, N., and Busenberg, E., 

2003, Aquifer susceptibility in Virginia, 1998–2000: U.S. 

Geological Survey Water-Resources Investigations Report 

03–4278, 58 p.

Zellmer, L.R., 1979, Development and application of a Pleistocene sea level curve to the Coastal Plain of southeastern 

Virginia: Williamsburg, Virginia, College of William and 

Mary, School of Marine Science, M.S. thesis, 85 p.

Zheng, C., and Wang, P., 1999, MT3DMS—A modular threedimensional multispecies transport model for simulation of 

advection, dispersion, and chemical reactions of contaminants in groundwater systems—Documentation and user’s 

guide: U.S. Army Corps of Engineers, U.S. Army Engineer 

Research and Development Center, Strategic Environmental 

Research and Development Program, 221 p.

Prepared by:

USGS Publishing Network

Raleigh Publishing Service Center

3916 Sunset Ridge Road

Raleigh, NC 27607

For additional information regarding this publication, contact:

Director

USGS Virginia Water Science Center

1730 East Parham Road

Richmond, VA 23228

phone: 804–261–2600

email: dc_va@usgs.gov

Or visit the Virginia Water Science Center Web site at:

http://va.water.usgs.gov

This publication is available online at http://pubs.usgs.gov/sir/2009/5039/.
Heywood and Pope—Simulation of Groundwater Flow in the Coastal Plain Aquifer System of Virginia—Scientific Investigations Report 2009–5039

Printed on recycled paper