#

308

Royal Valley Creek

Crouch Creek

Bjorgen

Creek

BROWNSVILLE HIGHWAY

Liberty 

Bay

Port

Orchard

Dogfish Bay

 Tide

Flats

Shallow

Lagoon

303

KEYPORT

VIRGINIA POINT

Naval Base

Kitsap

Liberty Bay

surface-water 

site

U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigation Report 2020-5066

Prepared in cooperation with the Department of the Navy, Naval Facilities Engineering 

Command, Northwest

Variable-Density Groundwater Flow and Contaminant 

Transport, Operable Unit 1, Naval Base Kitsap, Keyport, 

Washington
Cover: Map showing location of study area topography and surface drainage, transport domain and 

Operable Unit 1, Naval Base Kitsap, Keyport, Washington.
Variable-Density Groundwater Flow and 

Contaminant Transport, Operable Unit 1, 

Naval Base Kitsap, Keyport, Washington

By Richard M. Yager, Wendy B. Welch, Alexander Headman, and Richard S. 

Dinicola

Prepared in cooperation with the Department of the Navy, Naval Facilities 

Engineering Command, Northwest

Scientific Investigations Report 2020–5066

U.S. Department of the Interior

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

DAVID BERNHARDT, Secretary

U.S. Geological Survey

James F. Reilly II, Director

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

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

natural hazards, and the environment—visit https://www.usgs.gov or call 1–888–ASK–USGS.

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

https://store.usgs.gov/.

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

Government.

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

as noted in the text. Permission to reproduce copyrighted items must be secured from the copyright owner.

Suggested citation:

Yager, R.M., Welch, W.B., Headman, A., and Dinicola, R.S., 2020, Variable-density groundwater flow and contaminant 

transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington: U.S. Geological Survey Scientific Investigations 

Report 2020–5066, 58 p., https://doi.org/10.3133/sir20205066.

ISSN 2328-0328 (online)
iii

Acknowledgments

The authors acknowledge Carlotta Cellucci (Remedial Project Manager for Naval Facilities 

Engineering Command, Northwest) and Michael Meyer (Battelle Memorial Institute) for their 

assistance in providing historical reports and data to U.S. Geological Survey authors. We thank 

Chuck Heywood and Eric Swain (U.S. Geological Survey) for comments and discussions that 

sharpened our thinking and presentation.
This page left intentionally blank
v

Contents

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

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

Introduction.....................................................................................................................................................3

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

Previous Studies ...................................................................................................................................5

Contaminant History .............................................................................................................................5

Description of Study Area ............................................................................................................................8

Geologic Setting ....................................................................................................................................8

Hydrogeologic Units .............................................................................................................................8

Hydrologic Setting ................................................................................................................................9

Water-Level Response to Tidal Fluctuations .........................................................................12

Recharge .....................................................................................................................................13

Simulation of Constant-Density Groundwater Flow ...............................................................................15

Model Design.......................................................................................................................................15

Hydrogeologic-Unit Flow (HUF2) Package ............................................................................16

Surface Water ............................................................................................................................18

Domain Boundaries ...................................................................................................................21

Surface Drainage .......................................................................................................................21

Groundwater Withdrawals .......................................................................................................21

Recharge .....................................................................................................................................21

Model Calibration................................................................................................................................22

Iterative Calibration Approach ................................................................................................22

Representation of Storage Values ..........................................................................................26

Simulation Results ..............................................................................................................................26

Parameter Values ......................................................................................................................26

Model A ..............................................................................................................................26

Model B ..............................................................................................................................26

Model Fit ......................................................................................................................................27

Simulated Groundwater Flow ..................................................................................................27

Simulation of Variable-Density Flow and Transport of Chlorinated Ethenes .....................................35

Model Design.......................................................................................................................................35

Seawater Simulation .................................................................................................................35

CVOC Simulation ........................................................................................................................37

Simulation Results ..............................................................................................................................37

Seawater Simulation .................................................................................................................37

Chlorinated Volatile Organic Compounds Simulation ..........................................................37

Discussion of Simulation Results ..............................................................................................................46

Containment of Chlorinated Volatile Organic Compounds Contaminants .....................<?>

Model Uncertainty .....................................................................................................................46

Effect of Biodegradation ...........................................................................................................46

Alternative Models ....................................................................................................................49

Effect of Tidal Fluctuations .......................................................................................................49

Implications for Remediation of Contaminant Sources .......................................................50
vi

Summary........................................................................................................................................................52

Constant-Density Groundwater-Flow Simulation ..........................................................................52

Variable-Density Groundwater Flow and Contaminant-Transport Simulations .......................54

Discussion of Simulation Results .....................................................................................................55

References Cited..........................................................................................................................................56

Appendix 1. Soil-Water Balance (SWB) Model ..................................................................................59

Figures

1. Map showing study area showing topography and surface drainage, transport 

domain and Operable Unit 1, Naval Base Kitsap, Keyport, Washington .............................4

2. Map showing locations of former landfill and data-collection sites used to 

monitor tidal fluctuations at Operable Unit 1, Naval Base Kitsap, Keyport, 

Washington ....................................................................................................................................6

3. Map showing generalized extent of total concentrations of chlorinated volatile 

organic compounds in groundwater at Operable Unit 1, Naval Base Kitsap, 

Keyport, Washington, 2004, and locations of persistent hotspots in 2017 ..........................7

4. Map showing surficial hydrogeology, cross section traces, and framework well 

locations near Operable Unit 1, Naval Base Kitsap, Keyport, Washington ......................10

5. Generalized geologic sections near Operable Unit 1, Naval Base Kitsap, 

Keyport, Washington ..................................................................................................................11

6. Generalized conceptual section showing relationship between Units 7, 8, and 9 

beneath Operable Unit 1, Naval Base Kitsap, Keyport, Washington .................................12

7. Graph showing predicted tides for Liberty Bay at Poulsbo, Washington (station 

9445719), in July and August 2017 ............................................................................................13

8. Graph showing surface water levels in Liberty Bay, the tide flats, the marsh 

pond and creek, and the shallow lagoon measured from July 12 to 22, 2018, at 

Operable Unit 1, Naval Base Kitsap, Keyport, Washington .................................................14

9. Graph showing water-level response to tidal fluctuations at selected 

monitoring wells from July 12 to 22, 2018 at Operable Unit 1, Naval Base Kitsap, 

Keyport, Washington ..................................................................................................................16

10. Map showing distribution of average annual recharge from precipitation based 

on the Soil-Water Balance (SWB) model near Keyport, Washington ...............................17

11. Graph showing tidal fluctuation in Dogfish Bay from July 19 to July 22, 2018, and 

tidal fluctuation in Keyport groundwater-flow model ...........................................................18

12. Map showing boundaries specified within Keyport, Washington, 

groundwater-flow model domain .............................................................................................19

13. Map showing recharge areas based on surficial geology delineated for model 

calibration, Naval Base Kitsap, Keyport, Washington .........................................................23

14. Map showing observations used in model calibration, Naval Base Kitsap, 

Keyport, Washington ..................................................................................................................24

15. Residual plots showing relations between observed and simulated values in 

models A and B ...........................................................................................................................28

16. Residual plots showing relations between simulated values and residuals in 

models A and B ...........................................................................................................................29
vii

17. Maps showing Keyport groundwater-flow model results ...................................................30

18. Vertical flow from (positive) to (negative) Qv sediments simulated by Keyport 

groundwater-flow Model A, Operable Unit 1, Naval Base Kitsap, Keyport, 

Washington ..................................................................................................................................33

19. Map showing domain and boundary conditions in transport simulations, 

Operable Unit 1, Naval Base Kitsap, Keyport, Washington .................................................36

20. Graph showing computed mass of seawater during 5,000-year, single-solute 

simulations with models A and B, Operable Unit 1, Naval Base Kitsap, Keyport, 

Washington ..................................................................................................................................38

21. Three-dimensional perspective view showing computed seawater 

concentrations in the Qv sediments at the end of 5,000-year, single-solute 

simulation with Model A, Operable Unit 1, Naval Base Kitsap, Keyport, Washington ...38

22. Vertical section C–C’ through tide flats showing computed seawater 

concentrations and transition zone between fresh water and seawater at 

the end of 5,000-year, single-solute simulation with Model A, Operable Unit 1, 

Naval Base Kitsap, Keyport, Washington ...............................................................................39

23. Graph showing computed mass of CVOCs during 50-yr, dual-solute simulations 

using total-variation-diminishing and implicit finite-difference solvers with 

Model A, and total-variation-diminishing solver with Model B, Operable Unit 1, 

Naval Base Kitsap, Keyport, Washington ...............................................................................39

24. Map and three-dimensional perspective view showing computed CVOC 

concentrations in the Qv sediments at the end of 50-year, dual-solute 

simulation with Model A ............................................................................................................41

25. Map and graphic showing computed CVOC concentrations in the Qv sediments 

at the end of 50-yr, dual-solute simulation with Model B ....................................................43

26. Vertical section from marsh pond through tide flats showing computed CVOC 

and seawater concentrations at the end of 50-yr, dual-solute simulation 

with Model A ...............................................................................................................................44

27. Map showing comparison of the extents of chlorinated volatile organic 

compounds migration with total-variation-diminishing and implicit 

finite-difference solvers with Model A showing effect of numerical dispersion ............45

28. Graph showing computed mass of chlorinated volatile organic compounds 

during original 50-yr, dual-solute simulation with Model A, and alternate 

simulation with biodegradation ................................................................................................47

29. Map showing computed chlorinated volatile organic compound (CVOC) 

concentrations in the Qv sediments at the end of 50-yr, dual-solute simulation 

with Model A including biodegradation: A, plan view, B, three-dimensional 

perspective view .........................................................................................................................48

30. Graph showing computed mass of chlorinated volatile organic 

compounds during original 50-year, dual-solute simulations using 

total-variation-diminishing and implicit finite-difference solvers with Model A, 

and alternate steady-state simulation C with lower tidal stages and alternate 

transient simulation D with transient representation of tidal fluctuations .......................51

31. Graph showing reduction in Chlorinated volatile organic compound mass with 

time resulting from three contaminant source-control scenarios with Model A 

that reduce CVOC mass-loading rates ....................................................................................51
viii

Tables

1. Hydrogeologic units defined in this study and correlation between 

hydrogeologic units and model layers of the Keyport and Kitsap groundwater 

flow models and previous investigations .................................................................................9

2. Groundwater level response to tidal fluctuations in monitoring wells at 

Operable Unit 1, Naval Base Kitsap, Keyport, Washington .................................................15

3. Water levels specified for surface water bodies at high, low, and mean tides ...............21

4. Parameter values in Keyport groundwater-flow models A and B specified or 

estimated through model calibration ......................................................................................25

5. Water budgets for the groundwater-flow and transport model domains 

simulated for the steady-state period with Model A by the confined 

MODFLOW-2005 and unconfined MODFLOW-NWT simulations ........................................34

6. Reduction in chlorinated volatile organic compounds mass after 50 years of 

contaminant source-control measures with Model A, Operable Unit 1, Naval 

Base Kitsap, Keyport, Washington ..........................................................................................52

Conversion Factors

U.S. customary units to International System of Units

Multiply By To obtain

Length

inch (in.) 2.54 centimeter (cm)

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Area

acre 0.4047 hectare (ha)

square foot (ft2) 0.09290 square meter (m2)

square mile (mi2) 2.590 square kilometer (km2)

Volume

gallon (gal) 3.785 liter (L)

gallon (gal) 0.003785 cubic meter (m3)

cubic foot (ft3) 0.02832 cubic meter (m3)

Flow rate

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

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

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

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

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

gallon per day (gal/d) 0.003785 cubic meter per day (m3/d)
ix

Mass

ounce, avoirdupois (oz) 28.35 gram (g)

pound, avoirdupois (lb) 0.4536 kilogram (kg)

ton, short (2,000 lb) 0.9072 metric ton (t)

ton, long (2,240 lb) 1.016 metric ton (t)

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)

Hydraulic gradient

foot per mile (ft/mi) 0.1894 meter per kilometer (m/km)

Transmissivity

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

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

 °F = (1.8 × °C) + 32.

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

 °C = (°F – 32) / 1.8.

Datums

Vertical coordinate information is referenced to the North American Vertical Datum of 1988 

(NAVD 88).

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

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

Supplemental Information

Specific conductance is given in microsiemens per centimeter at 25 degrees Celsius (μS/cm 

at 25 °C).

Concentrations of chemical constituents in water are given in either milligrams per liter (mg/L) 

or micrograms per liter (μg/L).
x

Abbreviations

1,1-DCA 1,1-dichloroethane

1,1-DCE 1,1-dichloroethene

1,1,1-TCA 1,1,1-trichloroethane

CA chloroethane

CVOC chlorinated volatile organic compounds

DOD U.S. Department of Defense

EPA U.S. Environmental Protection Agency

FD finite difference

PCE tetrachloroethene

ppm parts per million

SWB soil water balance

TCE chloroethene compounds trichloroethene

cis-DCE cis-1,2-dichloroethene

trans-DCE trans-1,2-dichloroethene

TVD total-variation-diminishing

USGS U.S. Geological Survey

VC vinyl chloride
Variable-Density Groundwater Flow and Contaminant 

Transport, Operable Unit 1, Naval Base Kitsap, Keyport, 

Washington

By Richard M. Yager, Wendy B. Welch, Alexander Headman, and Richard S. Dinicola

Abstract

Chlorinated volatile organic compounds (CVOCs) have 

migrated to groundwater beneath a former 9-acre landfill 

at Operable Unit 1 (OU-1) on Naval Base Kitsap, which 

was active from the 1930s through 1973 on the Keyport 

Peninsula, in Kitsap County, Washington. Biodegradation of 

CVOCs at OU-1 limits the mass of dissolved-phase CVOCs 

in groundwater that discharges to surface water, but contaminant concentrations up to 630 milligrams per liter persist 

in localized areas, likely from the dissolution of residual, 

non-aqueous phase liquids. Variable-density groundwaterflow and contaminant-transport models were developed using 

the SEAWAT-Version 4 computer program to simulate the 

direction and rate of groundwater flow in a 5.9 square-mile 

(mi2) - area surrounding the Keyport Peninsula, to estimate the 

CVOC mass in groundwater and the rate of mass loading, and 

to assess possible remedial activities at OU-1.

The study area is underlain by Quaternary deposits consisting of alternating glacial and interglacial sediments ranging 

from 500 to 1,500 feet (ft) thick. A hydrogeologic model 

delineated a sequence of 10 units including a relatively thin 

package (less than 100 ft) of recent sediments (Vashon Stade 

and younger) beneath the Keyport Peninsula that are underlain 

by the much thicker (more than 300 ft) Clover Park Aquitard, 

which overlies a confined, sea-level aquifer.

A shallow lagoon, tide flats, marsh (including a pond 

and outlet creek), and Dogfish Bay (a segment of Liberty Bay 

connected to Puget Sound) adjacent to OU-1 are the principal 

areas of groundwater discharge. Water levels in these surfacewater bodies fluctuate in response to semi-diurnal tidal cycles. 

Daily ranges in water-level fluctuations range from 10.5 ft 

in Dogfish Bay to 0.5 ft in the marsh. Groundwater-level 

response to tidal fluctuations in and near OU-1 is as much 

as 7.3 ft.

Three-dimensional (3D), constant-density groundwater 

flow was simulated using MODFLOW-2005 and MODFLOWNWT models that were run alternately in an iterative 

sequence to converge on a common numerical solution. 

MODFLOW-2005 employed the HUF2 package to represent 

thin, discontinuous hydrogeologic units beneath OU-1, and the 

aquifer system was treated as confined to obtain the numerical 

solution for hydraulic head. MODFLOW-NWT was employed 

to obtain a hydraulic head distribution while treating the aquifer system as unconfined, but it does not support the HUF2 

package. The water table computed by the MODFLOW-NWT 

model was used to specify as inactive the model cells in the 

MODFLOW-2005 model that were above the simulated water 

table, thereby reducing the saturated thickness in the confined 

simulation.

The groundwater-flow model utilized a variably spaced 

grid with 15 model layers, 120 rows, and 88 columns; cell 

widths ranged from 25 ft near OU-1 to 500 ft near model 

boundaries. Ten model layers of equal thickness represent the 

Quaternary Vashon Stade or younger deposits, while the bottom four model layers represents the Clover Park Aquitard and 

deeper units. The top model layer represents surface water. 

The groundwater-flow simulation included a steady-state 

period that provided initial conditions to a 31-day transient 

period that represented approximately 30 tidal cycles. Tidal 

fluctuations in seawater were approximated by a diurnal function that ranged from a single low-tide stage to a high-tide 

stage, corrected for the specific gravity of seawater. Water 

levels in the tide flats and the marsh pond and its outlet were 

specified as the average high, low, and mean levels measured 

during July and August 2018.

Model boundaries coincide with the grid of a previously 

developed regional Kitsap Peninsula groundwater-flow model, 

which was used to set head and flow boundaries in model 

layers representing the permeable interbeds within the Clover 

Park Aquitard and an underlying aquifer. Hydraulic heads at 

model boundaries within the Quaternary Vashon sediments 

were assigned to maintain the simulated water table below 

land surface. Groundwater discharge to surface water and 

groundwater withdrawals from deep wells were represented 

as internal boundaries, and net outflow from shallow groundwater withdrawals and infiltration from septic systems was 

included in the estimated recharge values. Recharge values 

computed by a Soil-Water Balance (SWB) model for the 

model domain were adjusted during model calibration.

Two alternative models (A and B) with differing storage properties were constructed to assess the ramifications 

of representing changes in storage from the fluctuating water 
2 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

with the confined numerical model. In Model A, a specific 

storage (Ss) value of 2.0×10-6 ft-1 was assigned to all hydrogeologic units, while in Model B, the Ss value of the uppermost unit (Qvr sediments) was increased to 0.1 to approximate 

an unconfined storage value associated with specific yield. A 

multiplicative scaling factor based on the Qvr unit thickness 

was applied in Model B to reduce the Ss values where the 

unit spanned multiple model layers. Unconfined storage near 

the water table at OU-1 was underestimated in Model A and 

overestimated in Model B.

Values of hydraulic conductivities and recharge were 

adjusted in both models A and B during calibration to match 

99 observations of 4 types of data: groundwater levels (69), 

water-level fluctuations in wells in response to tides (26), 

vertical hydraulic gradients (3) and streamflow (1). Pseudo 

water-table altitudes computed by MODFLOW-NWT were 

also specified at 23 upland locations. Observation weights 

accounted for differences in units or ranges of measurements 

such that the observations were weighted equally. Parameter 

values were initially estimated in the confined transient 

MODFLOW-2005 Model And then updated in the unconfined 

steady-state MODFLOW-NWT model. The MODFLOWNWT water table was then used to update inactive model cells 

in the MODFLOW-2005 simulation, and parameter estimation 

was reapplied. This procedure was continued until the volumetric water budgets of both simulations were nearly equal.

A total of 24 parameter values were specified, of which 

15 were estimated through nonlinear regression using UCODE 

in Model A and 16 were estimated in Model B. The remaining 

parameters were estimated through trial and error or assigned 

from literature values. The average recharge rate computed 

by the SWB model (16.7 inches per year [in/yr]) was reduced 

to 4.2 in/yr in Model A (comparable to the rate of 4 in/yr 

estimated for the northern part of the Kitsap groundwater-flow 

model) and 6.7 in/yr in Model B. This difference indicates 

that the SWB model likely overestimates recharge in uplands 

where relatively permeable soils at land surface are underlain by lower permeability sediments at depth. In these areas, 

infiltration through the soil profile likely discharges as interflow to local drainage, rather than entering the groundwaterflow system.

The standard error in model fit in both models A and B 

was similar for groundwater levels and pseudo water-table 

altitudes (7.8 and 11.1 ft, respectively), but better for waterlevel fluctuations in Model A (1.0 ft) than in Model B (1.3 ft). 

The directions of the three vertical hydraulic-gradient observations were simulated correctly, with errors ranging from 25 

to 66 percent, but simulated base flow to a perennial stream 

70 percent less than the measured value in Model A and 40 

percent less in Model B. The simulated hydraulic-head distribution in both models indicates that OU-1 is just north of a 

regional discharge area that limits groundwater flow from the 

landfill toward uplands located south and west of the tide flats 

and the marsh. Simulated vertical groundwater flow is downward beneath the uplands and upward beneath the lowlands 

and surface-water bodies.

Two SEAWAT simulations of steady-state, variabledensity flow and transient solute transport were conducted 

for both models A and B in a 1.3 mi2 domain surrounding the 

Keyport Peninsula. A 5,000-yr simulation with a single solute 

(seawater) established the initial distribution of seawater for a 

50-year simulation with two solutes (seawater and CVOCs). 

Both simulations used the design and parameter distributions 

from the corresponding groundwater-flow model. The transport model domain extends from surface water to the upper 

part of the Clover Park Aquitard.

Seawater concentrations in both transport simulations 

ranged from 0 percent for fresh water (specific gravity 1.0) to 

100 percent for seawater (specific gravity 1.022). Constantconcentration boundaries were specified in model layer 1 to 

represent seawater. An implicit finite-difference (FD) method 

with upstream weighting was used in the single-solute simulation which has 100-day time steps in the 5,000-year simulation. The predicted distributions of seawater equilibrated and 

the simulated transition zone between fresh water and seawater was about 100 ft wide near the tide flats west of OU-1 in 

Model A, where saline water overlies fresh water at shallow 

depths, and 125 ft in Model B.

Inflow of CVOCs at OU-1 was represented in dualsolute simulations through mass-loading rates (MT-1) at 

seven persistent contaminant sources, adjusted to reproduce 

the maximum observed CVOC concentrations in 2017. The 

expected biodegradation of CVOCs, initially neglected in the 

transport simulation, was addressed in an alternative transport 

simulation. The total-variation-diminishing (TVD) method 

was used to minimize numerical dispersion using 0.5-day time 

steps in the 50-yr simulation. The predicted CVOC plume 

stabilized at a contaminant mass of 256 kilograms (kg) and 

a mass-loading rate of 150 grams per day (g/d) in Model A. 

Most of the contaminant mass (134 g/d) discharged to small 

creeks draining the southern edge of OU-1, while the remainder discharged to the marsh and tide flats. A northern plume 

with 30 percent of the CVOC mass migrated toward the marsh 

pond and creek, and a southern plume with 70 percent of the 

mass that migrated toward small creeks south of OU-1. The 

predicted CVOC mass was 310 kg with a mass-loading rate of 

73 g/d in Model B, in which groundwater flow through deeper 

sediments beneath OU-1 was restricted by higher vertical 

anisotropy. The two contaminant plumes intermingled and 

were diverted westward to Dogfish Bay by higher hydraulic 

conductivity values in Model B, and less CVOC mass (48 

g/d) discharged to the small creeks and more mass discharged 

to the marsh pond and creek (7.8 g/d) and the tide flats 

(160 mg/d).

Results of groundwater-flow and contaminant-transport 

simulations in both models A and B indicate that geologic 

and hydraulic barriers limit the migration of CVOCs from 

the former landfill at OU-1, including the low-permeability 

Clover Park Aquitard, upward hydraulic gradients beneath 

the Keyport Peninsula, underflow from upland areas south 

of OU-1, and the density difference between fresh and saline 

water beneath the tide flats and in Dogfish Bay. Uncertainty 
Introduction 3

in the rates of recharge and hydraulic properties produces 

uncertainty in the estimated mass-loading rate of CVOCs at 

OU-1. Higher rates of recharge and groundwater flow require 

larger mass-loading rates to maintain a stable configuration of 

the contaminant plume.

Reductive dechlorination of CVOCs at OU-1 was represented in an alternative 50-yr transport simulation using a 

first-order, removal rate-constant of 10-3 d-1 that corresponds 

to a half-life of 850 days, as computed from the CVOC 

contaminant mass (256 kg) predicted by the transport simulation with Model B and an estimated biodegradation rate (150 

g/d). A stable 77-kg CVOC plume formed within 10 years of 

migration, and one half of the CVOC mass was lost through 

biodegradation (75 g/d). Nearly all the remainder (71 g/d) 

was discharged by the southern plume to a small creek south 

of OU-1. The areas of the northern and southern plumes were 

67 percent and 32 percent less, respectively, than the areas 

computed by the original transport simulation.

Both models A and B provide an adequate fit to hydraulic and contaminant data, but additional data are required 

to determine which model more accurately predicts actual 

contaminant migration at OU-1. Both models predict similar 

contaminant masses and biodegradation rate constants. Model 

predictions differ in terms of some of the hydraulic parameters 

and in the CVOC mass-loading rates and contaminant plume 

configurations. While Model B better represents the effects 

of tidal fluctuations on unconfined storage and uses recharge 

rates closer to the SWB estimated rates, Model A better predicts the groundwater-level response to tidal fluctuations.

Two alternative 50-yr, contaminant-transport simulations 

with Model A assessed the sensitivity of predicted CVOC 

migration to the depiction of tidal fluctuations. Mean tides 

were lowered in seawater, the tide flats, and Marsh Creek 

in a 50-yr steady-state simulation with the total-variationdiminishing (TVD) solver to match average water levels 

measured during July and August 2018. Diurnal tidal fluctuations were also represented in a 50-yr transient simulation with 

the finite difference solver that computed 600 repetitions of a 

31.05-day transient period. Both alternative representations of 

tidal fluctuations had a small effect on the simulated CVOC 

plumes, because the period of the fluctuations (12 hr) is rapid 

compared with the groundwater velocity (20 ft/yr) beneath 

OU-1, so tidal fluctuations only slightly alter groundwaterflow paths.

Three 50-yr, steady-state contaminant-transport simulations with Model A predicted the effects of potential sourcecontrol scenarios on the contaminant plumes by reducing 

CVOC mass-loading rates by factors of 50, 90, and 99 

percent, respectively. Under each scenario the mass of CVOCs 

decreased rapidly and then slowly approached a new dynamic 

equilibrium. The reductions of CVOC mass in groundwater 

were 46, 86, and 95 percent, respectively, but sizes of the 

contaminant plumes were largely unchanged.

Introduction

Chlorinated volatile organic compounds (CVOCs) have 

migrated to groundwater beneath a former 9-acre landfill 

at Operable Unit 1 (OU-1) of Naval Base Kitsap in Kitsap 

County, Washington. Naval Base Kitsap is located on the 

Keyport Peninsula, which protrudes into an extension of Puget 

Sound called Liberty Bay (fig. 1). The three predominant 

groundwater contaminants are the chloroethene compounds 

trichloroethene (TCE), cis-1,2-dichloroethene (cis-DCE), 

and vinyl chloride (VC); six other contaminants of concern include the chloroethene and chloroethane compounds 

tetrachloroethene (PCE), trans-1,2-dichloroethene (transDCE), 1,1,1-trichloroethane (1,1,1-TCA), 1,1-dichloroethane (1,1-DCA), chloroethane (CA) and 1,1-dichloroethene 

(1,1-DCE). Recently another solvent 1-4 dioxane was also 

detected at the site (Department of the Navy, 2015). A need for 

remedial action was identified because some of the contaminants present a potential risk to humans, primarily through 

drinking contaminated groundwater or through ingesting 

seafood harvested from contaminated surface water (URS 

Consultants, Inc., 1998).

The U.S. Navy began a cooperative effort with the U.S. 

Geological Survey (USGS) in 1995 to evaluate the effectiveness of natural attenuation processes for removing and 

controlling the migration of CVOCs in groundwater at OU-1 

through field and laboratory studies conducted from 1996 

through 2000 (Dinicola and others, 2002). Additional collection of geochemical and contaminant concentration data was 

conducted by the USGS from 2001 through 2004 (Dinicola, 

2006). The USGS studies demonstrated that biodegradation 

of CVOCs in shallow groundwater at OU-1 is substantial 

and prevents most of the mass of dissolved-phase CVOCs in 

groundwater beneath the landfill from discharging to surface 

water (Dinicola, 2006). However, dissolved-phase contaminant concentrations in the hundreds of milligrams per liter 

continue to persist in localized areas of OU-1. These data 

indicate that residual sources of chloroethenes in the form of 

non-aqueous phase liquid remain at the site, and that biodegradation is only partly effective at reducing the dissolved-phase 

contaminants that are generated from these sources.
4 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

#

308

Royal Valley

Creek

Crouch Creek

Bjorgen

Creek

BROWNSVILLE HIGHWAY

Liberty 

Bay

Port

Orchard

Dogfish Bay

 Tide

Flats

Shallow

Lagoon

T

26

N

T

25

N

R 01 E

303

KEYPORT

VIRGINIA POINT

Naval Base

Kitsap

Liberty Bay

surface-water 

site

0 0.5 1 MILES

0 0.5 1 KILOMETERS

122°39' 122°38' 122°37'

47°

42'

47°

41'

Base map modified from U.S. Geological Survey and other digital

data at 1:26,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

WASHINGTON

Bremerton

Seattle

Poulsbo

Figure location

Kitsap Peninsula 

regional 

groundwater 

model domain

EXPLANATION

Naval Base Kitsap (NBK)

Flow model domain

Transport model domain

Former landfill - Operable Unit 1 (OU-1)

Figure 1. Location of study area showing topography and surface drainage, transport domain and Operable Unit 1, Naval Base Kitsap, 

Keyport, Washington.
Introduction 5

In 2018, an additional USGS effort was begun to simulate 

variable-density groundwater flow and contaminant transport 

in the vicinity of OU-1 using a revised hydrogeologic model 

of the site and a refined delineation of persistent contaminant 

sources, based on data collected by the U.S. Navy in 2016 

and 2017 (Battelle Memorial Institute, 2018). The groundwater flow and contaminant transport models were developed 

using SEAWAT-Version 4, a computer program based on 

MODFLOW and MT3DMS to simulate three-dimensional, 

variable-density groundwater flow coupled with multi-species 

solute transport (Langevin and others, 2007). These models 

were used to simulate the direction and rate of groundwater 

flow near OU-1, estimate the CVOC mass in groundwater and 

the rate of mass loading, and assess possible remedial activities at OU-1. Additional tasks were undertaken in support of 

the model development, including measurements of the timing 

and magnitude of water-level changes in response to tides and 

estimates of groundwater recharge on and around the Keyport 

Peninsula.

Purpose and Scope

This report describes the hydrogeology in the vicinity of 

the Keyport Peninsula, including a discussion of the surficial 

geology and the hydrologic setting. The report includes generalized geologic sections that show the principal hydrogeologic 

units that underlie the study area. The report also presents 

measured water-level changes in response to tides and the estimation of groundwater recharge using the Soil Water-Balance 

(SWB) model (Westenbroek and others, 2010). The report 

documents the design and calibration of a three-dimensional, 

groundwater-flow model that was used to simulate the pattern and rate of groundwater flow and to estimate values of 

hydraulic properties of the hydrogeologic units that control the 

flow of groundwater. The report also documents the design of 

variable-density, groundwater flow and contaminant transport 

models that were used to simulate the migration of CVOCs 

from localized sources at OU-1. Applications of additional 

model simulations are presented to show the effects of tidal 

fluctuations and biodegradation on the migration of CVOCs. 

Finally, applications of model simulations are presented that 

predict the fate of contaminants at OU-1 after three hypothesized reductions of localized sources through potential 

remedial activities.

Previous Studies

Geologic, water-level, and contaminant-concentration 

data are available from site characterization studies and 

long-term monitoring conducted for the U.S. Navy using 

direct-push and auger drilling and groundwater sampling from 

monitoring wells (Department of the Navy, 2015; Battelle 

Memorial Institute, 2018). CVOC mass-degradation rates and 

fluxes to surface water were calculated by Dinicola and others 

(2002) and Dinicola (2006) from contaminant concentrations 

and water-level contour maps. A three-dimensional characterization of hydraulic head and groundwater flow rates in the 

575 square mile (mi2) area surrounding the Keyport Peninsula 

is available from a steady-state, numerical groundwaterflow model of the Kitsap Peninsula developed by Frans and 

Olsen (2016).

Contaminant History

The landfill at OU-1 is on the narrow strip of land connecting the Keyport Peninsula to the mainland and is adjacent 

to tide flats that are an extension of Dogfish Bay, which is a 

segment of Liberty Bay (figs. 1 and 2). The OU-1 landfill is 

unlined at the bottom and was constructed in a marsh. The 

landfill was the primary disposal area for domestic and industrial wastes generated by Naval Base Kitsap from the 1930s 

through 1973. Paints, thinners, solvents, acids, dried sludge 

from a wastewater-treatment plant, and other industrial wastes 

were disposed in the landfill (Department of the Navy, 1993).

The most concentrated disposal area for waste paints 

and solvents was at the southern end of the landfill where the 

maximum total CVOC concentrations in groundwater were 

over 100 milligrams per liter (mg/L) in 2017 (fig. 3). The 

term “total CVOCs,” as used in this report, is the sum of the 

positively detected concentrations (in either mass or moles 

per unit volume) of six chloroethene and three chloroethane 

compounds of concern at OU-1. Total CVOC concentrations 

were less, but still over 10 mg/L in the central part of the 

landfill in 2017. Results from long-term monitoring indicate 

that the extent of the CVOC plume in groundwater at OU-1 

has been relatively stable, and that a phytoremediation remedy 

employed at both the northern and southern ends of the landfill 

has not been as effective as originally anticipated (Department 

of the Navy, 2015). Non-aqueous phase liquids in the subsurface are assumed responsible for both the persistence and 

magnitude of CVOC concentrations in these “hotspots,” 

locations where CVOC concentrations are substantially higher 

than in surrounding areas. Additional characterization of the 

hotspots in 2017 through vertical profiling using direct-push 

sampling revealed maximum CVOC concentrations of 630 

mg/L as much as 28 ft below land surface in the southern part 

of OU-1, and 120 mg/L more than 15 ft below land surface in 

the central part of OU-1.
6 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

Shapely

Creek

Marsh Creek

Seasonal Creek

308

Shallow

lagoon

Dogfish Bay

Tide

flats

Marsh 

Pond

Keyport Saltwater

Park

Tide flats

Marsh Creek

Marsh

Pond

Shallow

Lagoon

A ST

KEYS RD

NE STATE HWY 308

PRICHARD ST

BRADLEY RD

SHARPLEY RD

BROWNSVILLE HWY NE

CANAGA RD

TORPEDO RD

JENSON RD

MCKITTRRICK RD

UNNAMED

KESTREL PL NE

NE TAGHOLM RD

WRIGHT RD

HUNNICUTT RD

2ND ST

GADBERRY ST

NE FIR ST

GARNETT WAY

NE PLEASANT ST

MW1-9

MW1-2

1MW-4

1MW-1

MW1-60

MW1-51

MW1-50

MW1-47

MW1-45

MW1-43

MW1-39

MW1-38

MW1-29

MW1-25

MW1-20

0 250 500 FEET

0 50 100 METERS

122°37'45" 122°37'30"

47°

42'

00"

47°

41'

45"

Base map modified from U.S. Geological Survey and other digital

data at 1:4,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

Tidal well

Surface-water site

Former landfill - Operable Unit 1 (OU-1)

Phytoremediation plantation

Model grid

Naval Base Kitsap (NBK)

Figure 2. Locations of former landfill and data-collection sites used to monitor tidal fluctuations at Operable Unit 1, Naval Base 

Kitsap, Keyport, Washington.
Introduction 7

Shapely

Creek

Marsh Creek

Seasonal Creek

Tide

Flats

Marsh 

Pond

P1-9

P1-8

P1-6 P1-7

P1-5

P1-4

P1-3

P1-1

P1-10

MW1-6

MW1-5

MW1-4

MW1-2

1MW-1

MW1-3

MW1-60

MW1-58

MW1-56 MW1-57

MW1-55

MW1-54

MW1-53

MW1-52

MW1-51

MW1-50 MW1-49

MW1-48

MW1-47

MW1-46

MW1-45

MW1-44

MW1-43

MW1-42

MW1-41

MW1-29

MW1-28

MW1-25

MW1-15

MW1-41

MW1-18

MW1-17

MW1-20

KEYS RD

BRADLEY RD

BROWNSVILLE HWY NE

SHARPLEY RD

NE TAGHOLM RD

NE STATE HWY 308

GADBERRY ST

TORPEDO RD

0 125 250 FEET

0 25 50 METERS

122°37'45" 122°37'40" 122°37'35" 122°37'30"

47°

41'

55"

47°

41'

50"

47°

41'

45"

Base map modified from U.S. Geological Survey and other digital

data at 1:2,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

Well with concentration data

Former landfill - Operable Unit 1 (OU-1)

Naval Base Kitsap (NBK)

Chlorinated volatile organic compound 

concentration in milligrams per liter

(downstream extents are approximate)

0.01 to 0.5

0.51 to 1.0

1.01 to 10

10.01 to 100

More than 100

Figure 3. Generalized extent of total concentrations of chlorinated volatile organic compounds (CVOCs) in groundwater at 

Operable Unit 1, Naval Base Kitsap, Keyport, Washington, 2004, and locations of persistent hotspots in 2017. [Concentrations 

reflect the likely maximums expected at any depth based on 2017 data (Battelle Memorial Institute, 2018).]
8 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

Description of Study Area

Groundwater flow was simulated in a 5.9 mi2 area surrounding the Keyport Peninsula extending from Liberty Bay 

to the uplands 2 mi south of the peninsula. The area includes 

Virginia Point, a residential area with domestic wells located 

0.6 mi northwest of OU-1 across Dogfish Bay. Study area 

boundaries were delineated to coincide with the model grid of 

the Kitsap Peninsula used by Frans and Olsen (2016), so simulated groundwater levels from their model could be used to set 

boundary conditions for the Keyport models. Variable-density 

groundwater flow and contaminant transport was simulated 

within a smaller 1.3 mi2 transport domain within the larger 

5.9 mi2 study area that is centered on OU-1 and surrounds the 

Keyport Peninsula (fig. 1).

Geologic Setting

A summary of the major Quaternary events on the 

Kitsap Peninsula surrounding Naval Base Kitsap is given 

below and is based on the work of Garling and others (1965), 

Vaccaro and others (1998), Jones (1999) and Welch and others (2014). Plate tectonics and Pleistocene glaciations are the 

two major processes that have shaped the present-day Puget 

Lowland. The advance and retreat of continental glaciers in 

the Pleistocene Epoch of the Quaternary Period beginning 

over 100,000 years ago left behind as much as 3,000 ft of 

unconsolidated deposits in the Puget Lowland with 500–1,500 

ft of unconsolidated deposits commonly found beneath the 

Kitsap Peninsula. Underlying these unconsolidated deposits 

are Tertiary bedrock units consisting primarily of sedimentary claystone, siltstone, sandstone, beds of coal, and volcanic rocks.

The Puget Lobe of the Cordilleran ice sheet advanced 

and retreated several times into the Puget Lowland from the 

mountains of British Columbia since the beginning of the 

Quaternary Period and left behind a complex sequence of 

alternating glacial and interglacial sediments. The Vashon 

Stade of the Fraser Glaciation (12,000–18,000 years ago) was 

the most recent and extensive of the major advances that covered the Kitsap Peninsula. The sequence of glacial deposits in 

each glacial interval typically contains advance outwash sand 

and gravels, overlain by glacial till (poorly sorted mixture of 

indurated clay, silt, sand and gravel) and covered by recessional outwash sand and gravels. Each major glacial interval 

was followed by an extended interglacial period where fluvial, 

lacustrine, bog and marsh deposition dominated. Interglacial 

deposits typically consist of clay, silt, or discontinuous lenses 

of sand and gravel or peat. Since the end of the last glaciation, 

erosion has been the dominant geomorphic process affecting 

the Kitsap Peninsula. Alluvial sediments, typically sands and 

gravels, have been deposited by streams in valleys and marsh 

deposits formed in low-lying poorly drained areas.

Hydrogeologic Units

For this study, unconsolidated deposits were described 

by a hydrogeologic model that grouped the deposits into 

ten hydrogeologic units consisting of aquifers and confining units (table 1; figs. 4 and 5) that were based on lithologic 

(depositional facies, grain size, and sorting) and hydrologic 

characteristics, and relative stratigraphic position. Generally, 

in the Puget lowland, saturated glacial outwash or coarsegrained interglacial deposits form aquifers, whereas deposits 

such as till and fine-grained interglacial deposits form confining units. The hydrogeologic units defined in this study 

are similar to units defined by previous investigations of the 

Kitsap Peninsula (Welch and others, 2014) and the Keyport 

Peninsula (URS Consultants Inc., 1997); however, they have 

been reinterpreted for this study, and the extents or thicknesses 

of individual units are not directly comparable. Drillers’ logs 

for wells containing descriptions of lithology were the primary 

source of information for the definition of hydrogeologic units 

in the vicinity of the Keyport Peninsula. The hydrogeologic 

analysis was conducted using a geographic information system 

(GIS) that included locations and lithologic information for 

more than 300 wells, surficial hydrogeology from Welch and 

others (2014), and lidar-derived digital land-surface altitudes.

Unit 1. Unit 1 is present at land surface and is primarily 

composed of Vashon recessional outwash (Qvr), artificial fill 

and topsoil found near OU-1, and surficial alluvial deposits. 

This unit consists of silt, sand, and gravel deposits that closely 

follow river valleys and the shoreline, as well as modified 

land/artificial fill composed of heterogeneous material ranging from clay to gravel and trash. The average thickness of 

Unit 1 throughout the study area is about 8 ft and is 10 ft 

beneath OU-1.

Unit 2. Unit 2 consists of recent marsh and estuary 

deposits (Qm) which are composed of uncompacted silt or 

silty sand that generally underlie Unit 1 and are distinguished 

in drillers’ logs by dark organic material. Unit 2 extends from 

Dogfish Bay, underlies the tide flats and OU-1, and extends 

eastward beneath the shallow lagoon. The average thickness 

of Unit 2 throughout the study area is 14 ft and is 6 ft beneath 

OU-1. Unit 2 was divided into marsh (Qmarsh) and estuary 

(Qestuary) units for the purpose of estimating hydraulic conductivity in model calibration.

Unit 3. Unit 3 is present at land surface throughout much 

of the study area and consists primarily of Vashon till (Qvt). 

The till is a dense mix of sand and gravel in a clay matrix. 

Unit 3 mantles the hilltops south of the Keyport Peninsula 

and is not present at OU-1. The average thickness of Unit 3 

throughout the study area is 32 ft but varies spatially from a 

thin veneer on hilltops to more than 150 ft near Port Orchard.

Unit 4. Unit 4 primarily consists of Vashon advance 

outwash (Qva). Unit 4 is found beneath Unit 3 throughout 

the study area and is exposed at land surface along coastal 

bluffs and hillsides. It consists of well-sorted sand or sand 
Description of Study Area 9

and gravel with lenses of silt and clay. The average thickness 

of Unit 4 throughout the study area is 40 ft but is only 6 ft 

beneath OU-1.

Unit 5. Unit 5 is a laterally discontinuous sandy silt with 

some clay (Qvsemi) that only occurs in the area surrounding 

OU-1. Unit 5 forms a semi-confining layer between Units 4 

and 6 that can only be distinguished from those units by its 

finer grain size. The average thickness of Unit 4 throughout 

the study area is 8 ft and is 6 ft beneath OU-1.

Unit 6. Unit 6 consists of sand, gravel, and silt deposits 

(Qvinter) in the area surrounding OU-1 that overlies a continuous clay layer known locally as the Clover Park Aquitard, 

Lawton Clay (Battelle Memorial Institute, 2017) or QC1 

(Welch and others, 2014). The average thickness of Unit 6 

throughout the study area is 12 ft and is 11 ft beneath OU-1.

Units 7, 8 and 9. The Clover Park Aquitard is divided 

into three parts: the upper clay (Unit 7), the permeable 

interbeds (Unit 8), and the lower clay (Unit 9). Units 7 and 

9 (QC1) are thick and laterally extensive consisting of early 

Vashon, low-permeability lacustrine silt and clay and underlying interglacial fine-grained deposits. Units 7 and 9 are present 

throughout the study area with average thicknesses of 232 ft 

and 121 ft beneath OU-1 for Unit 7, and average thicknesses 

of 92 ft and 60 ft within OU-1 for Unit 9. The permeable 

interbeds (QCpi) in Unit 8 are locally continuous deposits of 

sand and gravel that lie between Units 7 and 9 (fig. 6). The 

average thickness of Unit 8 where it exists in the study area 

is 65 ft.

Unit 10. Unit 10 is a confined aquifer commonly referred 

to as the sea-level aquifer (QA1) or Salmon Springs Drift 

(Welch and others, 2014) that is present beneath most of 

study area. Unit 10 consists primarily of sand and gravel with 

silt interbeds and is the lowermost unit represented in the 

groundwater-flow model.

Hydrologic Setting

The topography of the study area indicates that surface 

water in the northern two-thirds of the area drains northward 

to Liberty Bay, while the southern third drains southward to 

Royal Valley and Crouch Creek (fig. 1). Steep bluffs border 

the eastern edge of the study area along Port Orchard where 

groundwater discharges through springs at discrete locations. 

Groundwater level data are sparse in the uplands, but the 

direction of groundwater flow probably parallels the direction 

of surface-water drainage. A small perennial stream drains 

part of the uplands to a shallow lagoon on the east side of the 

Keyport Peninsula.

Table 1. Hydrogeologic units defined in this study and correlation between hydrogeologic units and model layers of the Keyport and 

Kitsap groundwater flow models and previous investigations.

[U.S. Geological Survey Kitsap Model: From Welch and others (2014). Stratigraphic units: From URS Consultants, Inc., 1997, table 3-2. Abbreviations: —, 

not applicable]

Keyport model (this report) U.S. Geological Survey 

Kitsap model Stratigraphic units

Transport model Flow model Hydrogeologic model Hydrogeologic units Unit identifier and 

hydrostatigraphy

Layer 1–Water Layer 1–Water — — —

Layer 2–11 Qv 

sediments

Layer 2–11 Qv 

sediments Unit 1 (Qvr, fill)

Qvr, Vashon recessional aquifer

A, fill and topsoil

Er, recessional aquifer

Unit 2 (Qm) B, marsh

C, estuary

Unit 3 (Qvt) Qvt, Vashon till confining unit Et, Vashon till

Unit 4 (Qva)

Qva, Vashon advance aquifer

Ea, Upper aquifer

H1, Upper aquifer

Unit 5 (Qvsemi) H2, middle aquitard

Unit 6 (Qvinter) H4, intermediate aquifer

Jo, intermediate aquifer

Layer 12 Layer 12 Unit 7 (QC1up) QC1, upper confining unit Km, Clover Park aquitard

— Layer 13 Unit 8 (QC1pi) QC1pi, permeable interbeds Ks, Clover Park coarse-grained 

zone

Layer 14 Unit 9 (QC1low) QC1, upper confining unit Km, Clover Park aquitard

Layer 15 Unit 10 (QA1) QA1, sea-level aquifer —
10 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

303

308

Liberty 

Bay

Port

Orchard

Dogfish Bay

Tide

Flats

Shallow

Lagoon

Brownsville Highway

Royal Valley

Creek

Crouch Creek

Bjorgen

Creek

Hoot Creek

Unit 4

Unit 4

Unit 1

Unit 1

Unit 3

Unit 3

Unit 3

Unit 3

Unit 7

Unit 7

Unit 2

Unit 1

Unit 7

Unit 1 A

A'

B

B'

0 0.5 1 MILES

0 0.5 1 KILOMETERS

122°39' 122°38' 122°37' 122°36'

47°

42'

47°

41'

Base map modified from U.S. Geological Survey and other digital

data at 1:26,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

Framework well

Cross section trace

U.S. Geological Survey model extent

Former landfill - Operable Unit 1 (OU-1)

Naval Base Kitsap (NBK)

Surficial hydrogeology

Unit 1 (Qvr, fill)

Unit 2 (marsh)

Unit 2 (estuary)

Unit 3 (Qvt)

Unit 4 (Qva)

Unit 7 (QC1up)

A A'

Figure 4. Surficial hydrogeology, cross section traces, and framework well locations near Operable Unit 1, Naval Base Kitsap, 

Keyport, Washington.
Description of Study Area 11

Unit 1

Unit 1

Unit 2

Unit 2

Unit 3

Unit 4

Unit 5 Unit 4

Unit 5 Unit 6

Unit 6

Unit 7

Unit 7

Unit 8

Unit 8

OU1 extent

OU1 extent

NORTHWEST SOUTHEAST

SOUTHWEST NORTHEAST

A A'

B B'

Feet

Feet

0 

50 

−50

0

50

−50

−100

−100

Unit 9

Dogfish 

Bay Tide Flats Shallow Lagoon

Marsh 

Pond

VERTICAL EXAGGERATION X 7

NORTH AMERICAN VERTICAL DATUM OF 1988

0 0.1 0.2 MILES

0 0.1 0.2 KILOMETERS

EXPLANATION

Unit 1 (Qvr,fill)

Unit 2 (Qmarsh)

Unit 2 (Qestuary)

Unit 3 (Qvt)

Unit 4 (Qva)

Unit 5 (Qvsemi)

Unit 6 (Qvinter) 

Unit 7 (QC1up) 

Unit 8 (QC1pi)

Unit 9 (QC1low)

Hydrogeologic unit

Figure 5. Generalized geologic sections near Operable Unit 1, Naval Base Kitsap, Keyport, Washington.
12 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

Two creeks flow into the marsh where the Keyport 

Peninsula connects to the mainland. A single perennial creek 

(referred to as Marsh Creek) drains a perennial pond in the 

marsh and discharges into the tide flats of Dogfish Bay (fig. 2). 

Marsh Creek flows through a culvert and tide gate and into the 

tide flats. Water exchange between the tide flats and Dogfish 

Bay is constrained by a weir in the State Route 308 causeway 

that has a sill elevation of 0.3 ft. During high tides, water 

levels in the tide flats are at the same water-level altitude as 

Dogfish Bay, but during low tides the water level in Dogfish 

Bay is lower than the sill elevation and the water level in the 

tide flats. Also, the tide gate at Marsh Creek outlet automatically closes during higher tides and prevents seawater from 

flowing upstream beyond the marsh pond, resulting in a higher 

water level in the tide flats than in the marsh. During periods 

of no surface-water runoff, there is perennial freshwater inflow 

into the marsh derived from groundwater discharge.

The shallow lagoon, the marsh (including the pond and 

associated creeks), the tide flats, and Dogfish Bay are the 

principal areas of groundwater discharge near OU-1. Shallow 

groundwater at OU-1 originates as recharge from precipitation that falls on OU-1, in addition to lateral flows from 

the Keyport Peninsula to the northeast and from the Kitsap 

Peninsula uplands to the southwest, and vertical flow from the 

underlying aquifer of Unit 10 (QA1). Groundwater level data 

collected by the Navy and USGS (Welch and others, 2014; 

URS Consultants Inc., 1998) indicate that groundwater flow 

northeast of OU-1 converges with flow from the southwest, 

and then flows beneath the landfill toward the marsh. Some 

groundwater discharges to the shallow lagoon southeast of 

OU-1, but the flow divide between discharge areas in the 

marsh and shallow lagoon is located to the southeast of OU-1. 

In general, hydraulic gradients indicate groundwater flow is 

vertically downward beneath most of the Kitsap Peninsula 

and the uplands to the southwest, and flow is upward beneath 

the Keyport Peninsula and marsh. Deep wells drilled below 

the Clover Park Aquitard (QC1) for Navy and public supplies 

were reported to be artesian.

Water-Level Response to Tidal Fluctuations

The water-level response to tidal fluctuations in the 

vicinity of OU-1 was measured by Opatz and Dinicola (2019)

to determine optimal times to sample groundwater flowing 

from the landfill. Groundwater and surface-water levels were 

measured at 19 monitoring wells and 5 adjacent surface-water 

sites to provide 15-min time-series data from July 12, 2018, to 

August 8, 2018, a period that included both neap and spring 

tides. Data from 15 of those wells and all 5 surface-water sites 

were used during model calibration (figs. 1 and 2). The five 

surface-water sites included Liberty Bay on the north side of 

the Keyport Peninsula, the tide flats, the pond in the marsh and 

its outlet creek, and the shallow lagoon.

Daily tidal fluctuations in Liberty Bay follow a semidiurnal pattern that includes two high tides and two low 

tides each day, with a higher and lower high tide and a 

higher and lower low tide (fig. 7). The sequence of the 

high-high, low-high, high-low, and low-low tides within a 

tidal cycle changes over the course of a year. The means of 

the high-high tides and the low-low tides were computed 

from the 2017 predicted tides for Liberty Bay (National 

Oceanic and Atmospheric Administration [NOAA] station 

9445719 at Poulsbo, Washington; fig. 1;National Oceanic and 

Atmospheric Administration, 2020), and were corrected to the 

North American Vertical Datum of 1988 (NAVD88) datum at 

the site (Opatz and Dinicola, 2019) for comparison to surface 

water and groundwater levels near the Keyport Peninsula. The 

mean high-high tide in Liberty Bay was 7.94 ft NAVD88 and 

the mean low-low tide was -2.54 ft, resulting in a mean range 

of 10.53 ft. The mean tide in Liberty Bay in 2017 was 5.24 

ft NAVD88.

Water-level data for surface-water bodies near the 

Keyport Peninsula (fig. 8) indicate that the daily tidal fluctuation in Liberty Bay ranged from 9 to 16 ft during the observed 

period. The daily tidal cycle had a mean length of 24.84 hr. 

The semi-diurnal tidal pattern was reflected in water-levels 

measured in the tide flats and Marsh Creek, with the minimum 

water level in the tide flats restricted to 0.3 ft by the sill that 

connects the tide flats to Dogfish Bay, yielding a daily waterlevel fluctuation range in the tide flats of only 7–8 ft. The 

daily range in water-level fluctuation in Marsh Creek was only 

about 1.5 ft, and the water level in the marsh pond showed 

less than 0.5 ft fluctuations in response to tides around a mean 

pond elevation of 6.0 ft. The water-level response to tides 

in the shallow lagoon was greatly attenuated by a sill at the 

lagoon outlet and was only discernible at the highest tides.

Groundwater-level response to tidal fluctuations at the 

19 monitoring wells ranged from 7.3 ft at MW1-39 to 0.1 ft 

at MW1-50 and MW1-51 (table 2). Four wells (MW1-10, 

MW1-20, MW1-41, and P1-10) exhibited a groundwaterlevel response to tidal fluctuations less than 0.1 ft. The lag in 

groundwater-level response ranged from 15 min at MW1-39 

to 255 min at MW1-50, with the lag time increasing with 

diminishing groundwater-level response (fig. 9). Groundwaterlevel data collected at high and low tides during a previous 

Unit 7

Unit 7 Unit 8

Unit 9

Figure 6. Relationship between Units 7, 8, and 9 beneath 

Operable Unit 1, Naval Base Kitsap, Keyport, Washington.
Description of Study Area 13

study (URS Consultants, Inc., 1997) indicated responses to 

tidal fluctuations ranging from 5.1 to 0.5 ft in an additional 14 

monitoring wells (table 2).

The horizontal propagation of changes in groundwater 

level from tidal fluctuations in Dogfish Bay to the monitoring 

wells is governed by the hydraulic diffusivity of the aquifer 

sediments, as determined using the equation:

(D = K/Ss) (1)

where

K is the hydraulic conductivity [LT-1], and ‘

Ss is the specific storage [T-1].

The hydraulic diffusivity is related to the period of the 

tidal fluctuations and the attenuation in amplitude of the 

fluctuations at an observed distance x from the shoreline. The 

hydraulic diffusivity can be calculated from the following 

equation (Keller and others, 1989; Yager and Kappel, 1998), 

given as a solution by Carslaw and Jaeger (1959) to the governing equation for one-dimensional flow induced through isotropic and homogeneous media by a fluctuating head boundary 

approximated as sinusoidal function:

 D = [

x

_2π

T0

 ][ln (

_a

A)

−2] (2)

where

T0 is the period of the fluctuations, and

a/A is the attenuation of the original amplitude A.

The hydraulic diffusivity was computed for the 12 wells 

where fluctuations in groundwater levels were greater than 

0.1 ft, assuming a 12-hr period in tidal fluctuation. Values of 

D ranged from 2.2×105 to 18×105 feet squared per day (ft2/d), 

with a geometric mean value of 5.4×105 ft2/d (table 2).

Recharge

The magnitude and distribution of groundwater recharge 

in the 5.9 mi2 study area was computed using national spatial 

datasets and the SWB model (Westenbroek and others, 2010), 

which calculates spatial and temporal variations in recharge 

based on climatological data, and soil and landscape properties. The SWB model is a deterministic model that uses gridded data and physically based parameters to apportion water 

derived from daily precipitation and snowmelt into surface 

runoff, evapotranspiration (ET), recharge, and water storage 

in the soil column. Recharge is computed as surplus water in 

excess of the soil field capacity (storage), a product of available water capacity and root depth in the soil column.

July August

2017

-4

-2

0

2

4

6

8

10

12

14

Water level, in feet

Figure 7. Predicted tides for Liberty Bay at Poulsbo, Washington (station 9445719), in July and August 2017.
14 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

SWB calculates water-balance components in the following equation using a modified Thornthwaite Mather soil-waterbalance approach at daily time steps in a grid of model cells:

 recharge = (precipitation + snowmelt + inflow)

− (interception + outflow + ET)−  Δst

where

outflow (Surface runoff) is determined by the 

National Resource Conservation (NRCS) 

curve numbers (Miller and White, 1998; 

Natural Resources Conservation Service, 

2018) that have been modified for the 

Pacific Northwest by Gendasezek and 

Welch (2018);

inflow (runoff from upgradient adjacent cells) 

is computed using runoff directions 

determined with a flow-direction 

grid derived from a digital elevation 

model (DEM);

Interception is the specified amount of precipitation that is 

trapped and transpired by vegetation; and

ET (evapotranspiration) was computed using the method of 

Hargreaves and Samani (1985).

The SWB model constructed for the study area utilizes 

a 25-ft grid spacing, contains 8,652 model cells and requires 

five spatially distributed datasets including: daily climate data 

(precipitation, maximum and minimum air temperature), land 

cover, flow direction, hydrologic soil group and available 

soil-water capacity (app. 1). Daily calculations of the waterbalance components were made for the period 1980–2015. The 

daily recharge rates were then summed and averaged to obtain 

spatially distributed, mean-annual recharge rates for the study 

area. The recharge rates reflect the spatial distribution of land 

cover and soil properties in the study area (fig. 10). The average recharge rate computed for the study area was 16.7 in/yr.

-8

-6

-4

-2

0

2

4

6

8

10

12

12 13 14 15 16 17 18 19 20 21 22

Water-level elevation, in feet above the North American Vertical Datum of 1988

July 2018

Tide flats

Marsh creek

Marsh pond

Shallow lagoon

Liberty Bay

EXPLANATION

Figure 8. Surface water levels in Liberty Bay, the tide flats, the marsh pond and creek, and the shallow lagoon measured from 

July 12 to 22, 2018, at Operable Unit 1, Naval Base Kitsap, Keyport, Washington.
Simulation of Constant-Density Groundwater Flow 15

Simulation of Constant-Density 

Groundwater Flow

Constant-density groundwater flow and contaminant 

transport within the Keyport study area were simulated using 

three MODFLOW-based models. Three-dimensional (3D) 

groundwater flow in the entire 5.9 mi2 study area was simulated using both MODFLOW-2005 (Harbaugh, 2005) and 

MODFLOW-NWT (Niswonger and others, 2011) during 

model calibration, as described below. Three-dimensional, 

variable-density groundwater flow and contaminant transport in the 1.3 mi2 transport domain surrounding OU-1 were 

simulated using SEAWAT-Version 4 (Langevin and others, 

2007), a computer program that combines a modified version of MODFLOW-2000 (Harbaugh and others, 2000) with 

MT3DMS (Zheng and Wang, 1999). The model input, output, 

and executable files are available from Yager (2020).

The 3D constant-density, transient groundwater-flow 

model was calibrated using the parameter-estimation program UCODE (Poeter and Hill, 1998) to groundwater-level 

responses in 26 wells during a generalized 31-day period that 

represented 30 tidal cycles. A steady-state simulation that 

provided initial conditions to the 31-day transient period was 

also included in model calibration to match groundwater levels 

observed in 69 wells. Steady-state, groundwater levels for 54 

of the 69 wells at or near OU-1 were determined by averaging 5 rounds of measurements (each including measurements 

at low tide and high tide) made from August 1995 through 

September 1996 (URS Consultants Inc., 1997). Steady-state, 

groundwater levels for the other 15 wells farther from OU-1 

were approximated from one or two measurements made 

during well installation, or at various times over the past few 

decades (URS Consultants Inc., 1997, app. A). The calibrated 

constant-density, groundwater-flow model was the basis for 

the SEAWAT variable-density, contaminant-transport models 

applied to the transport domain surrounding OU-1, discussed 

further on.

Model Design

The design of the 3D groundwater-flow model was based 

on the 3D hydrogeologic model presented earlier. The model 

domain encompasses 5.9 mi2 and is divided into a variably 

spaced grid with 120 rows and 88 columns; cell widths ranged 

from 25 ft near OU-1 to 500 ft near Model Boundaries (fig. 1). 

The 15 model layers contain a total of 135,477 active cells. 

The bottom four model layers correspond to the bottom four 

hydrogeologic units represented in the geologic model [the 

Clover Park Aquitard (Unit 7; QC1) to the sea-level aquifer 

(Unit 10; QA1), table 1], while the upper six hydrogeologic 

units that represent Quaternary Vashon Stade (Qv) or younger 

deposits (referred to as Qv sediments herein) were subdivided 

into 10 model layers to better represent flow through the 

surficial sediments beneath OU-1. The top model layer (1) 

represents surface water, including the shallow lagoon, the 

marsh pond and its outlet, the tide flats, and the open marine 

waters of Dogfish and Liberty Bays and Port Orchard. The 

active model domain excludes the land area north of Liberty 

Bay, with the exception of the bottom model layer (15), which 

represents the sea-level aquifer (Unit 10; QA1) and extends 

throughout the model domain.

Table 2. Groundwater level response to tidal fluctuations in 

monitoring wells at Operable Unit 1, Naval Base Kitsap, Keyport, 

Washington.

[Abbreviations: ft, foot; ft2/d, square feet per day]

Well identifier Groundwater level 

response (ft)

Hydraulic diffusivity 

(105 ft2/d)

From Opatz and Dinicola (2019)

1MW-1 0.9 5.4

1MW-4 1.0 3.6

MW1-2 2.7 9.9

MW1-9 0.3 3.0

MW1-25 4.1 12.0

MW1-29 3.4 18.0

MW1-38 6.8 2.2

MW1-39 7.3 2.5

MW1-43 0.6 5.7

MW1-45 0.3 5.9

MW1-47 0.4 6.1

MW1-50 0.1 -

MW1-51 0.1 -

MW1-60 0.2 4.0

From URS Consultants, Inc. (1997)

MW1-18 1.4 -

MW1-21 0.5 -

MW1-22 0.9 -

MW1-23 2.2 -

MW1-28 3.5 -

MW1-30 5 -

MW1-32 0.5 -

MW1-34 0.6 -

MW1-36 5.1 -

MW1-37 5.1 -

MW1-40 1.5 -

PW1535HWY308 0.5 -
16 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

Hydrogeologic-Unit Flow (HUF2) Package

The HUF2 package (Anderman and Hill, 2000) was used 

for representing the thin discontinuous hydrogeologic units 

that are common beneath OU-1. The method entails incorporating hydraulic properties for the upper six hydrogeologic 

units described in the previous section (table 1) into flowmodel layers 2–11. The six units do not extend throughout the 

model domain and are discontinuous at various locations. The 

HUF2 package assigns hydraulic properties associated with 

each hydrogeologic unit to the appropriate flow-model layer 

based on the depth interval occupied by the unit. The HUF2 

package then determines the number of hydrogeologic units 

that intersect each flow-model layer and computes an average property value for each model cell. An alternative method 

for representing discontinuous hydrogeologic units would 

have entailed assigning each unit to a separate flow-model 

layer. Very large vertical hydraulic-conductivity values and 

small layer thicknesses would then be specified in model cells 

where the hydrogeologic units are missing to minimize the 

loss of hydraulic head while simulating flow across these cells. 

Application of this alternative method would have required 

the use of very small time steps in the variable-density flow 

and contaminant transport simulations described further 

on, however, and would have resulted in unreasonably long 

simulation times. One limitation of the HUF2 package is that 

its application is limited to MODFLOW-2000-based models, 

such as SEAWAT Version 4, because the package is not currently supported in more recent versions of MODFLOW, such 

as MODFLOW-NWT.

-8

-6

-4

-2

0

2

4

6

8

10

12

12 13 14 15 16 17 18 19 20 21 22

Dogfish Bay stage, in feet above the North American Vertical Datum of 1988

Water-level elevation in wells, in feet above the 

North American Vertical Datum of 1988

2

3

4

5

6

7

8

July 2018

MW1-25

MW1-1

MW1-45

Dogfish Bay

EXPLANATION

Figure 9. Water-level response to tidal fluctuations at selected monitoring wells from July 12 to 22, 2018, at Operable Unit 1, Naval 

Base Kitsap, Keyport, Washington.
Simulation of Constant-Density Groundwater Flow 17

303

308

Liberty 

Bay

Port

Orchard

Dogfish Bay

Tide

Flats

Shallow

Lagoon

BROWNSVILLE H GI HWAY

Royal Valley Creek

Crouch Creek

Bjorgen

Creek

Hoot Creek

0 0.5 1 MILES

0 0.5 1 KILOMETERS

122°39' 122°38' 122°37' 122°36'

47°

42'

47°

41'

T

26

N

T

25

N

Base map modified from U.S. Geological Survey and other digital R 01 E

data at 1:26,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

U.S. Geological Survey model extent

Former landfill - Operable Unit 1 (OU-1)

Not estimated

0.1 to 15

15.1 to 20

20.1 to 25

25.1 to 30

More than 30

Average annual recharge from 

precipitation, in inches per year

Naval Base Kitsap (NBK)

Figure 10. Distribution of average annual recharge from precipitation based on the Soil-Water Balance model near Keyport, 

Washington.
18 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

The thicknesses of the 10 model layers representing 

Qv sediments, which are equal within each vertical stack of 

model cells, were obtained by dividing the total thickness of 

Qv sediments at each grid cell by 10. This method resulted in 

a minimum model-layer thickness of 3 ft beneath OU-1. The 

maximum thickness of these flow-model layers ranged from 

17 ft beneath the Keyport peninsula to 18 ft thick beneath 

the uplands to the south. Hydraulic properties associated 

with each hydrogeologic unit were estimated during model 

calibration.

Surface Water

The groundwater-flow simulation consists of a steadystate period followed by a 31.05-day transient period divided 

into 30 daily tidal cycles. Each semi-diurnal tidal cycle was 

simplified to a diurnal tidal cycle and employed five stress 

periods to represent the approximate 25-hr tidal cycle in the 

open and connected waters of Dogfish and Liberty Bays and 

Port Orchard. Tidal fluctuations were approximated by a 

trapezoidal function with a period of 1.035 days (24.84 hr) 

that ranged from a single low-tide stage to a high-tide stage 

(fig. 11). The high and low tide stages were computed as the 

respective means of the high-high and low-low tides predicted 

for Liberty Bay at Poulsbo, Washington, in 2017, as presented 

earlier. This approximation ignores the intermediate high 

and low tides that are observed in the data, but allowed the 

groundwater-flow and contaminant-transport simulations to be 

conducted in a reasonable amount of time.

Water levels were specified for Dogfish and Liberty Bays 

and Port Orchard in model layer 1 at the beginning and end 

of each stress period as either a low, high, or mean (average 

of high and low) tide using the Time-Variant Specified-Head 

(CHD) package (Harbaugh, 2005) and corrected for the 

specific gravity of seawater (1.022) (Collias and others, 1974) 

(fig. 12A). Stress periods were divided into either two or four 

time steps, and water levels at each time step were interpolated 

by the CHD package from the specified water levels. This 

approximation overestimates the stage in Dogfish Bay during 

the tidal cycle but simplifies the specification of this boundary and permits the completion of a 50-yr transient transport 

simulation discussed further on. The mean tide predicted for 

Liberty Bay at Poulsbo, Washington, in 2017 was specified for 

Dogfish and Liberty Bays and Port Orchard during the steadystate period.

Mean and low tides were lower than the sediment surface 

at the top of the hydrogeologic model for some model cells 

that are exposed at the lower tide stages (for example, in mud 

flats surrounding Dogfish Bay). The minimum water level 

in these cells was specified as equal to the sediment-surface 

elevation and was assigned in model layer 2. The tide flats are 

an extension of Dogfish Bay and experience the same mean 

and high tides. At low tide, however, the water level in the tide 

flats was specified as 0.3 ft, the elevation of the sill connecting the two water bodies. Water levels in the marsh pond and 

its outlet were specified as the average high, low, and mean 

tides measured during July and August 2018 (table 3) using 

the Flow and Head Boundary (FHB) package. The water level 

in the shallow lagoon was specified as a constant value of 7.2 

ft with the CHD package, based on measurements by Opatz 

and Dinicola (2019). The FHB and CHD packages function 

similarly, but the application of the CHD package allowed 

constant-concentration boundaries to be specified for several 

water bodies in the variable-density and contaminant transport 

simulations described further on.

-8

-6

-4

-2

0

2

4

6

8

10

12

19 20 21 22

Water-level elevation, in ft above the

North American Vertical Datum of 1988

July 2018

High stage

0.512 days

1.035 days

Tidal cycle

Mean

stage

Mean

stage

Low

stage

Water level specified

in model simulations

Figure 11. Tidal fluctuation in Dogfish Bay from July 19 to 22, 2018, and tidal 

fluctuation in Keyport groundwater-flow model.
Simulation of Constant-Density Groundwater Flow 19

303

Royal Valley

Creek

Crouch Creek

Marsh

Creek

Seasonal Creek

Shapely

Creek

BROWNSVILLE H GI HWAY

Hoot Creek

Liberty 

Bay

Port

Orchard

Inactive cells

308

0 0.5 1 MILES

0 0.5 1 KILOMETERS

122°39' 122°38' 122°37' 122°36'

47°

42'

47°

41'

Base map modified from U.S. Geological Survey and other digital

data at 1:26,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

U.S. Geological Survey model extent

Former landfill - Operable Unit 1 (OU-1)

Specified head cells in Qv units

Surface water

Mud flat

Water table

Drains

Streams

Seeps

Wet areas

Naval Base Kitsap (NBK)

A

Figure 12. Boundaries specified within Keyport, Washington, groundwater-flow model domain: A, Surface water and Vashon 

Stade (Qv) units, and B, sea-level aquifer (QA1) unit.
20 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

303

308

Seasonal Creek

Liberty 

Bay

Port

Orchard

Dogfish Bay

Tide

Flats

Shallow

Lagoon

BROWNSVILLE H GI HWAY

Royal Valley Creek

Crouch Creek

Bjorgen

Creek

Hoot Creek

0 0.5 1 MILES

0 0.5 1 KILOMETERS

122°39' 122°38' 122°37' 122°36'

47°

42'

47°

41'

T

26

N

T

25

N

Base map modified from U.S. Geological Survey and other digital R 01 E

data at 1:26,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

Flow model domain

Former landfill - Operable Unit 1 (OU-1)

Model boundaries in QA units

Wells

Naval Base Kitsap (NBK)

B

Figure 12.—Continued
Simulation of Constant-Density Groundwater Flow 21

Domain Boundaries

Model domain boundaries coincide with the model grid 

of the 3D Kitsap groundwater-flow model of Frans and Olsen 

(2016), but they do not correspond to natural hydrologic 

boundaries. Hydraulic heads and flow rates computed by the 

Kitsap model were used to set head and flow boundaries in 

the Keyport model in two model layers below the top of the 

Clover Park Aquitard (Unit 7: QC1) with the FHB package. 

Hydraulic heads were specified at domain boundaries in model 

layers 13 and 15 representing the permeable interbeds (Unit 

8; QC1pi) within the Clover Park Aquitard and the sea-level 

aquifer (Unit 10; QA1), respectively (fig. 12B). Vertical flows 

exchanged between the sea-level aquifer and deeper model 

layers in the Kitsap model were also specified with the FHB 

package in model layer 15. Hydraulic heads computed by the 

Kitsap model were adjusted during calibration of the Keyport 

model to better match observed water levels in wells screened 

in the sea-level aquifer.

The Kitsap model was calibrated to water-level data over 

a much larger area than the Keyport model, and it predicts 

water levels above land surface along a large portion of the 

Keyport model boundary, particularly where the uplands slope 

downward to the edge of Dogfish and Liberty Bays and Port 

Orchard. Hydraulic heads were specified along southern and 

western grid boundaries of the Keyport model in a single 

model layer within the Qv sediments. The elevation and head 

at these boundary cells were computed as the top of the Clover 

Park Aquitard (Unit 7; QC1) plus one-third of the thickness 

of the Qv sediments. This boundary maintained the simulated 

water table in the Qv sediments below land surface.

Hydraulic heads were also specified with the CHD 

package along northern and eastern grid boundaries where 

the model domain intersects Liberty Bay and Port Orchard. 

Heads specified at these locations in the top model layer (1) 

correspond to the high, low, and mean tide stages. Heads were 

specified in model layer 12 (Unit 7; QC1) using values computed from the Kitsap model for the upper part of the Clover 

Park Aquitard. The remainder of the model domain boundaries 

were specified as no flow.

Surface Drainage

Groundwater discharge to surface drainage at land 

surface (fig. 12A) was represented by drains using the Drain 

(DRN) package which requires specification of the elevation and hydraulic conductance (discharge per foot of head 

difference between groundwater and the drain, ft2/d) of each 

boundary cell. Perennial streams delineated in the USGS 

National Hydrography Dataset (U.S. Geological Survey, 2014) 

were represented in model layer 2 using the bottom elevation 

of the layer and a conductance value of either 20 or 10 ft2/d for 

the main channels or tributaries, respectively. The small creek 

draining the southern edge of OU-1 and the ditch along the 

Brownsville Highway west of the marsh were also represented 

as drains using elevations scaled from site-scale topographic 

maps and conductance values of 10 ft2/d.

Groundwater seeps in the uplands were represented as 

drains with large conductance values (100 ft2/d) to minimize 

resistance to flow. Elevations of seeps along the bluffs at the 

eastern edge of the uplands above Port Orchard were set at 

the bottom of the Qv sediments in model layer 11 (fig. 12A). 

Drain elevations at land surface were specified in model layer 

2 for two wet areas along the south shore of Dogfish Bay. In 

addition, three lines of seeps were represented as drains in the 

Qv sediments where coarse-grained sediments of the Vashon 

advance aquifer (Unit 4; Qva) pinch out above the fine-grained 

sediments of till (Unit 3; Qvt). The elevations of these drains 

were specified as the top of the Unit 4.

Groundwater Withdrawals

Groundwater withdrawals from Welch and others (2014)

were specified at 20 pumped wells with the Well (WEL) 

package (fig. 12B). Groundwater withdrawals from 16 of the 

wells were specified in model layer 15 (Unit 10; QA1) and 

the remainder in model layer 13 (Unit 8; QC1p). Groundwater 

withdrawals from shallow domestic wells screened in Qv 

sediments and infiltration from domestic septic systems were 

not represented explicitly but are accounted for by representing net outflow (withdrawals minus infiltration) from these 

domestic systems in the estimated recharge values.

Recharge

Recharge computed by the SWB model (app. 1) was 

applied to model layer 2 with the Recharge (RCH) package. 

The model domain was divided into five areas based on the 

surficial geology to enable the adjustment of the recharge 

values during model calibration (fig. 13). Areas where coarsegrained sediments (Units 1 and 4; Qvr and Qva) were present 

at land surface were grouped together, as were areas where 

fine-grained sediments were present (Unit 3; Qvt). The finegrained sediment areas were divided further into two groups 

corresponding to upland and lowland areas. A fourth recharge 

area was delineated where the Clover Park Aquitard (Unit 7; 

Table 3. Water levels specified for surface water bodies at high, 

low, and mean tides.

Surface-water body Tidal stage

High Mean Low

Dogfish and Liberty Bays and Port Orchard 17.94 15.24 1-2.54

Tide flats 17.94 15.24 0.3

Marsh outlet 4.65 3.98 3.3

Marsh pond 6.4 6.1 5.8

1Water levels are corrected for seawater density in model simulation.
22 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

QC1) is present at land surface; recharge in these areas was 

set to the vertical-hydraulic conductivity specified for Unit 7. 

A fifth recharge area was delineated during model calibration 

around a small knoll in the lower part of the uplands southwest 

of Keyport Peninsula. The computed SWB recharge in this 

area was decreased by a factor of 0.1 in this area to limit water 

levels computed by the groundwater-flow simulation to below 

land surface, as described further on.

Model Calibration

The Keyport groundwater-flow model was calibrated by 

adjusting values of hydraulic properties to match 99 observations of four types of data: groundwater levels (69), waterlevel fluctuations in wells that respond to tides (26), vertical 

hydraulic gradients (3) and streamflow (1). Eleven groundwater level observations were in the QCpi and QA1 units, 

while the remainder were in Qv sediments. In addition, pseudo 

water-table altitudes were specified at 24 locations to compensate for the lack of data in the uplands and to ensure the predicted water table was below land surface (fig. 14). Estimation 

of the pseudo water-table altitudes is described below. A total 

of 24 parameter values were specified in the model (table 4), 

15 of which were estimated through nonlinear regression 

using UCODE (Poeter and Hill, 1998), a computer program 

that uses weighted least squares to minimize model residuals 

(difference between observed and simulated data). Model sensitivities corresponding to the remaining nine parameters were 

too small to allow estimation by regression, and values were 

assigned by trial and error, or taken from the literature.

Weights assigned to most observations included in the 

regression were chosen to account for differences in units 

(for example, groundwater levels and streamflow) or ranges 

of measurements (for example, 60 ft for groundwater levels 

and 7 ft for water-level fluctuations) and adjusted such that 

the observations had equal weight in the regression. Waterlevel fluctuations and vertical hydraulic-gradient measurements were weighted slightly larger than groundwater levels, 

and pseudo water-table altitudes were weighted slightly 

smaller. Downward hydraulic gradients were computed from 

the differences in groundwater levels measured in two pairs 

of wells screened in shallow and deep Qv sediments at the 

base of uplands south of the Keyport Peninsula (fig. 14). The 

vertical hydraulic gradient at a third location on the Keyport 

Peninsula was estimated from a water-table map prepared for 

the Qv sediments and results from the Kitsap groundwater 

flow model, which indicated an upward gradient beneath the 

peninsula and Liberty Bay. A single measurement of flow was 

recorded in the perennial stream that drains to the shallow 

lagoon on March 2019 following a period with little surface runoff.

Iterative Calibration Approach

Model calibration involved the application of two 

simulations that were run alternately in an iterative sequence: 

a combined steady-state and transient simulation using 

MODFLOW-2005 and a steady-state simulation using 

MODFLOW-NWT. The use of the HUF2 package necessitated the application of MODFLOW-2005, but that program 

requires the aquifer system to be treated as confined to obtain 

a numerical solution for hydraulic head in the study area, 

where steep slopes in the uplands transition abruptly to the 

level topography of the lowlands. When the aquifer system 

is treated as unconfined, the numerical solvers available in 

MODFLOW-2005 cannot converge to a solution for hydraulic head and dry cells (with heads below the cell bottoms) 

propagate throughout the model domain. Treating the aquifer 

system as confined overestimates underflow from the uplands 

to the lowlands, however, because the saturated thickness 

of each model layer is artificially increased. In contrast, the 

MODFLOW-NWT program can better simulate the transition in hydraulic head from the uplands to the lowlands while 

treating the aquifer system as unconfined. In the unconfined 

MODFLOW-NWT simulation, model layers above the water 

table are dry and the corresponding model cells are inactive. 

Therefore, underflow is computed using a more realistic value 

of saturated thickness. Unfortunately, the numerical solver in 

MODFLOW-NWT is not supported in the SEAWAT Version 4 

program that was later used to simulate variable-density flow 

and contaminant transport, which necessitated the iterative 

procedure described below.

Parameter values were estimated initially from a confined 

transient simulation using MODFLOW-2005 and then updated 

in the unconfined steady-state simulation using MODFLOWNWT. This entailed importing values for horizontal and 

vertical hydraulic-conductivity and recharge directly into the 

MODFLOW-NWT model, which does not support application of the HUF2 package. The water table computed by the 

MODFLOW-NWT simulation was then used to specify as 

inactive the model cells in the MODFLOW-2005 simulation 

that were above the simulated water table. Parameter estimation was then reapplied to the MODFLOW-2005 simulation, 

and the procedure was continued until the volumetric water 

budgets of both simulations were nearly equal. The pseudo 

water-table altitudes used in model calibration were taken 

from the water table computed by the MODFLOW-NWT 

simulation and updated after each iteration of this procedure.
Simulation of Constant-Density Groundwater Flow 23

303

308

Seasonal Creek

Liberty 

Bay

Port

Orchard

Dogfish Bay

Tide

Flats

Shallow

Lagoon

BROWNSVILLE H GI HWAY

Royal Valley Creek

Crouch Creek

Bjorgen

Creek

Hoot Creek

0 0.5 1 MILES

0 0.5 1 KILOMETERS

122°39' 122°38' 122°37' 122°36'

47°

42'

47°

41'

T

26

N

T

25

N

Base map modified from U.S. Geological Survey and other digital R 01 E

data at 1:26,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

U.S. Geological Survey model extent

Former landfill - Operable Unit 1 (OU-1)

Naval Base Kitsap (NBK)

Recharge areas

Units 1, 4 (Qvr, Qva)

Unit 3 (Qvt lowland)

Unit 3 (Qvt upland)

Unit 3 (Qvt upland knoll)

Unit 7 (QC1)

Figure 13. Recharge areas based on surficial geology delineated for model calibration, Naval Base Kitsap, Keyport, Washington.
24 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

303

308

Liberty 

Bay

Port

Orchard

Dogfish Bay

Tide

Flats

Shallow

Lagoon

BROWNSVILLE H GI HWAY

Royal Valley Creek

Crouch Creek

Bjorgen

Creek

Hoot Creek

OU1

0 0.5 1 MILES

0 0.5 1 KILOMETERS

122°39' 122°38' 122°37' 122°36'

47°

42'

47°

41'

T

26

N

T

25

N

Base map modified from U.S. Geological Survey and other digital R 01 E

data at 1:26,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

U.S. Geological Survey model extent

Former landfill - Operable Unit 1 (OU-1)

Gage site

Well type

Groundwater level

Water-level fluctuation

Pseudo water-table elevation

Groundwater level and water-level 

fluctuation

Naval Base Kitsap (NBK)

Figure 14. Observations used in model calibration, Naval Base Kitsap, Keyport, Washington.
Simulation of Constant-Density Groundwater Flow 25

Table 4. Parameter values in Keyport, Washington, groundwater-flow models A and B, specified or 

estimated through model calibration.

[Bold values signify coefficient of variation less than 100 percent. Abbreviations: ft/d, foot per day; –, not applicable]

Parameter Assigned 

value

Model A Model B

Estimated 

value

Coefficient 

of variation 

(percent)

Estimated 

value

Coefficient 

of variation 

(percent)

Hydraulic conductivity (ft/d)

Operable Unit 1(Qvr) – 0.9 – 8.2 83

Operable Unit 2 (Qm):

Qmarsh – 0.18 44 0.25 –

Qestuary – 1 – 19 –

Operable Unit 3(Qvt):

Lowland – 0.66 51 0.07 94

Upland – 0.15 34 0.24 80

Knoll – 10 51 10 3 –

Operable Unit 4(Qva):

Lowland – 0.9 – 0.9 –

Upland – 1.7 – 2 90

Unit 5(Qvsemi) 0.1 – – – –

Unit 6 (Qvinter) – 1 – 3.1 86

Units 7 & 9 (QC1) 0.01 – – – –

Unit 8 (QC1pi) – 1.7 – 1.7 –

Unit 10 (QA1) – 50 70 50 –

Anisotropy, ratio of horizontal to vertical hydraulic conductivity

Coarse-grained units: 2 – – – –

Units 1, 4, 6 & 8

(Qvr, Qva, Qvinter, QC1pi)

Unit 2 (Qm) 5 – – 1.5 –

Unit 3(Qvt) – 26 41 53 61

Unit 5(Qvsemi) 2 – – 4.6 92

Units 7 & 9 (QC1) 100 – – –

Recharge factor

Coarse-grained units:

Units 1 & 4 (Qvr & Qva) – 0.65 53 1.2 –

Unit 3 (Qvt):

Lowland – 0.17 – 0.23 –

Upland – 0.22 52 0.38 99

Knoll 0.1 – – – –

Unit 7 (QC1) 1.0E-04 – – – –

Storage (ft-1)

All units1 2.0E-06 – – – –

Unit 1 (Qvr)2 0.1 – –

1Except Unit 1 in Model B

2Model B only

3Specified in Model B
26 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

Representation of Storage Values

The assumption of confined conditions in the transient 

MODFLOW-2005 simulation does not represent the exchange 

of groundwater to and from storage that would occur near a 

fluctuating water table under unconfined conditions. However, 

hydraulic diffusivity (D) values computed from the observed 

groundwater fluctuations in wells near OU-1 indicate that the 

measured well responses occur under confined conditions. 

Assuming an expected hydraulic conductivity value for the 

aquifer sediments of 1 ft/d, the specific storage (Ss) computed 

from the geometric mean D value of 5.35×105 ft2/d is 1.9×10-6

ft-1, a value within the range of Ss values estimated for unconsolidated material (1.0×10-6 to 2×10-5 ft-1) by several studies 

(Yager and Fountain, 2001).

The ramifications of representing changes in storage 

from the fluctuating water table near OU-1 was assessed using 

two alternative models. In the first model (A), an Ss value of 

2.0×10-6 ft-1 was assigned to all hydrogeologic units. In the 

second model (B), the Ss value of Unit 1 (Qvr sediments), the 

uppermost hydrogeologic unit, was increased to 0.1 to approximate an unconfined storage value associated with specific 

yield. The larger Ss value in Model B overestimates storage 

in areas where the Qvr sediments span multiple model layers, 

however, because this larger value is operative in each layer, 

rather than only in the layer(s) in which the water table is 

fluctuating. The HUF2 package in MODFLOW-2005 does not 

allow the specification of different parameter values within the 

same hydrogeologic unit. A multiplicative scaling factor based 

on the Qvr unit thickness was applied to the Ss values of the 

Qvr sediments to compensate for this problem. The specified 

Ss values were smaller where Unit 1 was thicker and larger 

where the unit was thinner, such that the sum of Ss values in 

each stack of model cells representing Qvr sediments totaled 

to 0.1. This approach resulted in less overestimation of storage, but still overestimated Ss values in the deeper parts of the 

Qvr sediments where Unit 1 spanned multiple model layers.

Simulation Results

Parameter Values

Model A

Coefficients of variation (CV) for 8 of the 15 estimated 

parameters in Model A were less than 100 percent, indicating 

that the regression was relatively sensitive to these parameters 

and that these values are better estimated than those of the 

other parameters (table 4). Estimated hydraulic conductivity 

(K) values for the coarse-grained, permeable Qv sediments 

(Units 1, 4 & 6; Qvr, Qva and Qvinter) range from 0.9 to 1.7 

ft/d, and are lower than the geometric mean values (3.3 to 4.1 

ft/d) computed from 25 slug tests conducted in permeable 

Qv sediments near the OU-1 (URS Consultants, 1998). The 

estimated K value for the Unit 2 (Qestuary) was 1 ft/d; lower 

values had no effect on the regression. Estimated K values for 

the Unit 2 (Qmarsh) and Unit 3 (Qvt) are one order of magnitude lower than for the permeable units, with the exception 

of the till recharge area delineated around a small knoll in the 

lower part of the uplands southwest of the Keyport Peninsula 

(fig. 13). Water levels in two wells screened in shallow Qv 

sediments indicate that the slope of the water-table surface 

changes abruptly and becomes more gradual in this area. The 

hydraulic conductivity of Unit 3 was increased to 10 ft/d and 

recharge was reduced by a factor of 0.1 in this area to match 

the observed groundwater levels. These changes suggest that 

more permeable units occur in this area than those indicated 

by the hydrogeologic model. Values for vertical anisotropy 

(ratio of horizontal to vertical hydraulic conductivity) were 

specified from the literature, with the exception of anisotropy 

for Unit 3, which was estimated as 26:1.

The nonlinear regression decreased the recharge rates 

computed by the SWB model on the coarse-grained units by 

a factor of 0.65, and on Unit 3 (Qvt) by factors of 0.22 and 

0.17 in the uplands and lowlands, respectively. The average 

recharge rate in the flow model domain in Model A is 4.2 in/

yr. This value is much less than the average rate of 16.7 in/yr 

computed with the SWB model, but comparable to recharge 

rates estimated by Welch and others (2014) for the northern 

part of the Kitsap groundwater-flow model domain (4 in/yr). 

The recharge rates estimated for the Keyport groundwaterflow model predict base flow in the perennial stream draining to the shallow lagoon reasonably well, and higher 

recharge rates result in a simulated water table that is above 

land surface.

The lower recharge rates estimated by nonlinear regression for Model A suggests that the SWB model overestimates 

recharge to groundwater in uplands similar to those in the 

Keyport study area, where the topography dips steeply to 

sea level, and where relatively permeable soils at land surface are underlain by much lower permeability sediments at 

depth. Much of the precipitation that percolates through the 

soil root zone in the uplands likely discharges as interflow to 

local drainage features rather than infiltrating downward to 

the groundwater flow system. The recharge factors estimated 

by the nonlinear regression are correlated with the hydraulic 

conductivities estimated for Unit 3 (Qvt) and Unit 4 (Qva) at 

levels greater than 0.85, however, indicating that higher values 

could be specified for both recharge and K in these units with 

little change in model error.

Model B

Coefficients of variation (CV) for 8 of the 16 estimated 

parameters in Model B were less than 100 percent, although 

half of these parameters had CV values greater than or equal to 

90 percent (table 4). Estimated hydraulic conductivity (K) values for some of the coarse-grained, permeable Qv sediments 

(Units 1, 2 & 6; Qvr, Qestuary and Qvinter) were higher than 

in Model A and range from 3.0 to 19 ft/d. Estimated K values 

for lowland till (Unit 3, Qvt) were an order of magnitude less, 

while the remainder of the K values were comparable to those 
Simulation of Constant-Density Groundwater Flow 27

in Model A. Values for vertical anisotropy were estimated 

for Units 2, 3 and 5 (Qm, Qvt and Qvsemi, respectively), and 

were higher for Units 3 and 5, and lower for Unit 2 than in 

Model A. The higher K values estimated in Model B for the 

coarse-grained, permeable Qv sediments compensate for the 

higher storage values specified in Unit 2 (Qvr), so the hydraulic diffusivity (D) value remained unchanged at 5.35×105 ft2/d.

Estimated recharge rates were higher in Model B than 

in Model A. SWB recharge rates on the coarse-grained units 

were increased by a factor of 1.2 and decreased on Unit 3 

(Qvt) by factors of 0.38 and 0.23 in the uplands and lowlands, 

respectively. The average recharge rate in the flow model 

domain in Model B is 7.3 in/yr, slightly less than half the 

SWB recharge rate. The larger recharge factors estimated in 

Model B are correlated in the nonlinear regression with the 

higher hydraulic conductivities, as noted above.

Model Fit

Errors associated with groundwater-flow Model A for 

groundwater levels, water-level fluctuations and pseudo watertable altitudes are depicted in a series of residual plots (figs. 

15 and 16). The standard error in groundwater levels is 7.8 ft, 

which is 13 percent of the 60-ft measurement range, while the 

standard error in water-level fluctuations is 1.0 ft or 14 percent 

of the 7.3-ft measurement range. Errors in groundwater levels 

for the deeper hydrogeologic units below the Clover Park 

Aquitard (Units 7 & 9; QC1) show more scatter than for the 

Qv sediments. Groundwater levels in the vicinity of OU-1 are 

generally less than 20 ft where the agreement is quite good, 

although groundwater levels are overestimated in some areas 

at the base of the uplands south of the Keyport Peninsula. The 

standard error in pseudo water-table altitudes is 11.3 ft, which 

is 4 percent of the 320-ft measurement range. This indicates 

close agreement between the water tables computed by the 

MODFLOW-2005 and MODFLOW-NWT simulations.

Residual plots for Model A indicate little bias in model 

results for water-level fluctuations (fig. 16B), but there is 

a bias towards overestimating measured water levels and 

underestimating pseudo water-table altitudes (figs. 16A and 

C). These biases are in opposite directions and limit the ability of the nonlinear regression to further reduce model error. 

Additional groundwater-level data in the uplands would be 

required to better calibrate the groundwater-flow model.

Errors associated with groundwater-flow Model B for 

groundwater levels and pseudo water-table altitudes are 

similar to those with Model A, but the errors in water-level 

fluctuations are worse (fig. 15). The standard error in groundwater levels and pseudo water-table altitudes are 7.6 and 

11.2 ft, respectively, while the standard error in water-level 

fluctuations is 1.3 ft. The biases in water levels and pseudo 

water-table altitudes are also similar in Model B to those 

observed in Model A, but there is a distinct bias towards 

underestimating water-level fluctuations (fig. 16). The higher 

storage values in Model B dampen the water-level fluctuations, although the higher vertical anisotropy values in Model 

B (especially in Unit 5, Qvsemi) ameliorate this condition to 

some extent by decreasing the hydraulic connection between 

the upper and lower parts of the Qv sediments.

The spatial distribution of residuals indicates that the 

water-table altitude is underpredicted in Model A in the 

uplands and on Keyport Peninsula, but overpredicted in 

lowland areas near the base of the uplands (fig. 17A). Waterlevel fluctuations are generally underpredicted near OU-1 

but overpredicted near the tide flats. The spatial distribution 

of residuals in Model B is similar, but the water-table altitude is better predicted in lowland areas near the base of the 

uplands, and water-level fluctuations are underpredicted near 

OU-1. The directions of the three vertical hydraulic-gradient 

observations were simulated correctly in both models, downward beneath the uplands and upward beneath the Keyport 

Peninsula (fig. 17B). Errors in the simulated gradients ranged 

from 25 to 66 percent of the observed value. The simulated 

value of streamflow in the perennial stream that discharges to 

the shallow lagoon in Model A (0.18 ft3/s) was 70 percent less 

than the measured value of 0.6 ft3/s, while streamflow was 40 

percent less than the measured value in Model B (0.33 ft3/s).

Simulated Groundwater Flow

The water table simulated by groundwater-flow Model A 

for the steady-state period indicates that groundwater enters 

the model domain from the western boundary and flows 

northward to Dogfish and Liberty Bays and southward to the 

southern boundary (fig. 17B). There are two high points on 

the water-table surface that are separated by the watershed 

of the perennial stream that drains to the shallow lagoon. 

Groundwater also flows eastward and discharges along the 

steep bluffs that border Port Orchard. Groundwater flows 

northeastward from the uplands to the base of the Keyport 

Peninsula, where it discharges to the tide flats and the marsh 

and its outlet creek (fig. 17A). Another high point is evident on 

water-table surface on the western side of the peninsula. Some 

groundwater on the peninsula flows south and southwestward 

beneath OU-1, but the hydraulic gradient is gradual. The simulated water table indicates that OU-1 is located adjacent to a 

regional discharge area that prevents the flow of groundwater 

from the landfill toward areas south and west of the tide flats 

and the marsh. The water table simulated by Model B is similar, but the maximum predicted altitudes are 10 ft higher under 

Virginia Point and 20 ft under the Keyport Peninsula, because 

of increased recharge.
28 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

A

B

C

0

10

20

30

40

50

60

70

0 10 20 30 40 50 60 70

Observed head, in feet

Simulated head, in feet

0

1

2

3

4

5

6

7

8

0 1 2 3 4 5 6 7 8

Observed fluctuation, in feet

Simulated fluctuation, in feet

0

50

100

150

200

250

300

350

0 50 100 150 200 250 300 350

Observed head, in feet

Simulated head, in feet

Qv sediments:

Model A

Model B

QA/QC sediments

EXPLANATION

Figure 15. Residual plots showing relations between observed and simulated 

values in models A and B: A, water levels, B, water level fluctuations, and C, pseudo 

water-table altitudes.
Simulation of Constant-Density Groundwater Flow 29

-40

-30

-20

−10

0

10

20

30

40

0 10 20 30 40 50 60 70

Simulated head, in feet

-5

-4

-3

-2

-1

0

1

2

3

4

5

0 1 2 3 4 5 6 7 8

Residual, in feet

Simulated fluctuation, in feet

-40

-30

-20

-10

0

10

20

30

40

0 50 100 150 200 250 300 350

Simulated head, in feet

A

B

C

Qv sediments:

Model A

Model B

QA/QCp sediments

EXPLANATION

Figure 16. Residual plots showing relations between simulated values and residuals 

in models A and B: A, water levels, B, water level fluctuations, and C, pseudo 

water-table altitudes.
30 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

303

308

Liberty 

Bay

Port

Orchard

Dogfish Bay

Tide

Flats

Shallow

Lagoon

Brownsville Highway

Royal Valley Creek

Crouch Creek

Bjorgen

Creek

Hoot Creek

OU1

0 0.5 1 MILES

0 0.5 1 KILOMETERS

122°39' 122°38' 122°37' 122°36'

47°

42'

47°

41'

T

26

N

T

25

N

Base map modified from U.S. Geological Survey and other digital R 01 E

data at 1:26,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

U.S. Geological Survey model extent

Former landfill - Operable Unit 1 (OU-1)

Residuals, in feet—Groundwater levels (circles), 

water-level fluctuations (triangles), pseudo water-table 

elevations (squares) & vertical hydraulic gradients (stars), 

positive is under predicted, negative is over predicted. 

(Note: water-level fluctuation residuals are multiplied 

by ten, hollow circles are groundwater levels in QC1pi 

and QA1 units)

More than 10

5 to 10

0 to 5

−5 to 0

−10 to −5

Less than −10

Naval Base Kitsap (NBK)

A

Figure 17. Keyport groundwater-flow model results: A, distribution of residuals in groundwater-level, water-level fluctuations, 

pseudo water-table altitudes and vertical hydraulic-gradient observations; and B, simulated water table.
Simulation of Constant-Density Groundwater Flow 31

303

308

330

320

310

300

290

280 270

260 250

240

230

220

210

200

180

160 150

140 130

100

80

50

310

300

290

280

270

260

250

240

230

220

210

200

190

180

180

160

150

100

50

30

220

120

100

80

50

140

130

30

20

190

120

110

80

120

50

60

40

120

60

3020

40

20

10

210

200

50

40

130

120 30

20

30

20

190

180

170 170

90 80

70

70

40

30

20

20

90

40

20

200

140

10

10

1 0

1 0

230

1 6 0 140

60

60

60

60

20

Seasonal Creek

Liberty 

Bay

Port

Orchard

Dogfish Bay

Tide

Flats

Shallow

Lagoon

BROWNSVILLE H GI HWAY

Royal Valley Creek

Crouch Creek

Bjorgen

Creek

Hoot Creek

0 0.5 1 MILES

0 0.5 1 KILOMETERS

122°39' 122°38' 122°37' 122°36'

47°

42'

47°

41'

T

26

N

T

25

N

Base map modified from U.S. Geological Survey and other digital R 01 E

data at 1:26,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

U.S. Geological Survey model extent

Former landfill - Operable Unit 1 (OU-1)

Naval Base Kitsap (NBK)

Simulated water table elevation, in feet

10 to 60

70 to 120

130 to 180

190 to 250

260 to 330

B

Figure 17.—Continued
32 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

The Clover Park Aquitard (Unit 7, QC1), which serves 

as a hydraulic barrier in both models A and B, limits the flow 

of groundwater from the Qv sediments to and from the deeper 

hydrogeologic units. Vertical groundwater flow is generally at 

least an order of magnitude less than horizontal flow, except 

along the western boundary of the model domain where the 

flows are nearly equal. Groundwater flow from the Qv sediments is downward beneath the uplands and upward beneath 

Keyport Peninsula and Dogfish and Liberty Bays and Port 

Orchard in Model A (fig. 18). A transition from downward to 

upward flow occurs southeast of the tide flats and marsh at the 

base of the peninsula. A similar distribution of flow is evident in Model B. The simulated distributions of vertical flow 

indicate that upward flow beneath the landfill would prevent 

the migration of dissolved contaminants from OU-1 toward 

deeper hydrogeologic units.

Water budgets simulated for the steady-state period in 

Model A by the confined MODFLOW-2005 and unconfined 

MODFLOW-NWT simulations for the groundwater-flow and 

transport model domains are shown in table 5. The water budgets generally agree closely, although there is more underflow 

from Model Boundaries in the confined MODFLOW-2005 

simulation that overestimates the saturated thickness of Qv 

sediments. Most of the inflow to the groundwater-flow model 

domain is from recharge and underflow along the western 

model boundary. Most of the outflow is discharge to streams 

and underflow along the southern model boundary. The vertical flow to and from deeper hydrogeologic units is specified 

at the bottom model boundary in Unit 10 (QA1; sea-level 

aquifer). The pattern is similar for the transport model domain 

surrounding the Keyport Peninsula, but underflow from the 

uplands to the southern transport boundary is a smaller component of the overall budget. The vertical flow is through the 

Clover Park Aquitard (Unit 7, QC1) which underlies the transport model domain. Underflow from the uplands to transport 

model domain is overestimated by 35 percent by the confined 

MODFLOW-2005 simulation, and the simulated water table is 

as much as 8 ft above land surface in a 4-acre area south of the 

marsh. As a result, the magnitude of groundwater flow from 

the uplands to the marsh is overestimated in this area.

Most water-budget components for Model B are similar to those for Model A, although there is 75 percent more 

recharge and streamflow predicted in Model B. There is correspondingly more underflow predicted to the transport domain 

in Model B (9,360 ft3/d) than in Model A (4,814 ft3/d). The 

agreement between the confined MODFLOW-2005 simulation 

and the unconfined MODFLOW-NWT simulation with Model 

B is comparable to that of Model A. The water table simulated 

by Model B is at or below land surface in the 4-acre area south 

of the marsh, where Model A simulated the water table above 

land surface, suggesting a better model fit in this area.
Simulation of Constant-Density Groundwater Flow 33

303

308

Liberty 

Bay

Port

Orchard

Dogfish Bay

Tide

Flats

Shallow

Lagoon

BROWNSVILLE H GI HWAY

Royal Valley Creek

Crouch Creek

Bjorgen

Creek

Hoot Creek

0 0.5 1 MILES

0 0.5 1 KILOMETERS

122°39' 122°38' 122°37' 122°36'

47°

42'

47°

41'

T

26

N

T

25

N

Base map modified from U.S. Geological Survey and other digital R 01 E

data at 1:26,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

U.S. Geological Survey model extent

Former landfill - Operable Unit 1 (OU-1)

Naval Base Kitsap (NBK)

Vertical flow

Upward

Downward

Figure 18. Vertical flow from (positive) to (negative) Qv sediments simulated by Keyport groundwater-flow Model A, Operable Unit 

1, Naval Base Kitsap, Keyport, Washington.
34 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington Table 5. Water budgets for the groundwater-flow and transport model domains simulated for the steady-state period with Model A by the confined MODFLOW-2005 and unconfined MODFLOW-NWT simulations. MODFLOW-2005 simulation MODFLOW-NWT simulation Groundwater flow domain Groundwater flow domain Source Inflow volume (ft3/d) Discharge Outflow volume (ft3/d) Source Inflow volume (ft3/d) Discharge Outflow volume (ft3/d) Recharge from pre- cipitation 103,810 Streams 88,080 Recharge from precipitation 102,280 Streams 83,348 Underflow from up- gradient areas 85,175 Underflow to downgradi- ent areas 111,740 Underflow from upgradient areas 71,370 Underflow to down- gradient areas 99,472 Vertical flow from deeper units 27,275 Vertical flow to deeper units 4,406 Vertical flow from deeper units 27,275 Vertical flow to deeper units 4,406 Seawater1 2,276 Seawater1 13,328 Seawater1 95 Seawater1 12,830 Wells 982 Wells 982 Total 218,536 Total 218,536 Total 201,020 Total 201,038 Transport domain Transport domain Source Inflow vol- ume (ft3/d) Discharge Outflow volume (ft3/d) Source Inflow volume (ft3/d) Discharge Outflow volume (ft3/d) Recharge from pre- cipitation 18,340 Streams 13,209 Recharge from precipitation 18,308 Streams 12,367 Underflow from up- gradient areas 4,814 Underflow to downgradi- ent areas 272 Underflow from upgradient areas 3,545 Underflow to down- gradient areas 337 Vertical flow from deeper units 677 Vertical flow to deeper units 228 Vertical flow from deeper units 649 Vertical flow to deeper units 232 Seawater2 91 Seawater2 10,215 Seawater2 67 Seawater2 9,644 Total 23,922 Total 23,924 Total 22,570 Total 22,579

1Dogfish and Liberty Bays and Port Orchard.2Dogfish and Liberty Bays.
Simulation of Variable-Density Flow and Transport of Chlorinated Ethenes 35

Simulation of Variable-Density Flow 

and Transport of Chlorinated Ethenes

Two SEAWAT simulations of steady-state, variabledensity flow and transient solute transport were conducted 

for the 1.3 mi2 transport domain surrounding the Keyport 

Peninsula using models A and B. The first 5,000-yr simulation 

used a single solute (seawater) to establish the initial distribution of seawater for the second 50-yr simulation that consisted 

of two solutes, seawater, and CVOCs. Both SEAWAT simulations used the same model design and parameter distributions as in the MODFLOW-2005 groundwater-flow models A 

and B discussed previously but required the specification of 

additional model features to represent transport as described 

further on. The model input, output, and executable files are 

available from Yager (2020).

In the variable-density flow and transport simulations 

described herein, the flow and transport equations were explicitly coupled using a one time-step lag and solved alternately 

with a mass-balance error less than 0.01 percent. An implicit 

finite-difference (FD) method with upstream weighting was 

used to solve the advection-dispersion equation in the singlesolute simulation of seawater. The selection of the FD method 

allowed the use of relatively large time steps (up to 100 days) 

in the 5,000-year simulation and resulted in run times of 4 

hours on a 3.2-GHz Zeon processor. The total-variationdiminishing (TVD) method with a Courant number of 1.0 was 

used to solve the advection equation in the dual-solute simulation of seawater and CVOCs in order to minimize numerical 

dispersion. The TVD method required much smaller time steps 

(0.5 days) in the 50-yr simulation and resulted in run times of 

seven hours. The potential effect of numerical dispersion on 

model results was explored through alternative simulations 

discussed later in this report.

Model Design

The transport-model domain (fig. 19) extends from the 

top of the groundwater-flow model (model layer 1) to the 

upper part of the Clover Park Aquitard (Unit 7, QC1, model 

layer 12) (table 1). Solute concentrations predicted by transport simulations for Unit 7 are disregarded, however, because 

the thickness of model layer 12 is over 100 ft beneath OU-1, 

greater than the numerical dispersion criteria allows. Unit 

7 would need to be discretized into several model layers in 

order to obtain accurate solute concentrations from transport 

simulations. Four additional model parameters (porosity and 

longitudinal, horizontal-transverse, and vertical-transverse 

dispersivities) were specified in both transport simulations. 

A porosity value of 0.2 was specified for all hydrogeologic 

units and a porosity of 1.0 was specified for surface water. 

Longitudinal dispersivity parallel to the direction of groundwater flow was specified as proportional to the grid cell size, 

in accordance with the observation that measured values of 

dispersivity increase with the scale of observation (Zheng and 

Bennett, 2002). Longitudinal dispersivity αl was computed for 

each grid cell as:

αl = 

(_

Δr*Δc)1/2

2 (4)

where

Δr and Δc are the row and column dimensions (L) of 

the cell.

The maximum αl value was limited to 50 ft. The computed αl

 value (12.5 ft) for the smallest grid cells (25 × 25 ft) 

is within the high portion of the range of reported field values 

in Gelhar and others (1992) and therefore conservative, as 

dispersion could be over-estimated. The computed αl value (50 

ft) for the largest cells (351 × 500 ft) falls within the trend of 

the data in Gelhar and others (1992). The largest grid cells are 

located at the eastern edge of the Keyport Peninsula and do 

not affect the simulated migration of CVOCs at OU-1. Values 

of horizontal-tranverse and vertical-tranverse dispersivity perpendicular to the direction of groundwater flow were specified 

as ratios to the longitudinal dispersivity. Ratios of 0.1 and 0.01 

for the horizontal-tranverse and vertical-tranverse dispersivity, 

respectively, were specified in accordance with recommended 

values in Zheng and Bennett (2002).

Seawater Simulation

The transport equations solved by the SEAWAT simulations are written in terms of concentration (ML-3), and 

additional boundary conditions regarding solute concentrations are specified for the transport-model domain. Seawater 

concentrations in both the single- and dual-solute simulations 

ranged from zero percent for fresh water with a specific gravity of 1.0 to 100 percent for seawater in surface-water bodies with a specific gravity of 1.022. Constant-concentration 

boundaries of 100 percent were specified in model layer 1 to 

represent seawater in Dogfish and Liberty Bays, the tide flats 

and Port Orchard (fig. 19). Constant-concentration boundaries 

of 50 percent were specified in model layer 2 to represent pore 

waters in mud flats surrounding Dogfish Bay. The initial seawater concentration in all other grid cells within the transport 

domain was specified as zero in the single-solute, seawater 

simulation. The computed seawater concentrations in model 

cells beneath the surface-water bodies increased during the 

5,000-year simulation, and the computed distribution of seawater concentration at the end of the simulation served as the 

initial seawater-concentration distribution for the dual-solute, 

CVOC simulation.
36 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

0 750 1,500 FEET

0 150 300 METERS

Base map modified from U.S. Geological Survey and other digital

data at 1:12,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

122°38'00" 122°37'30" 122°37'00"

47°

41'

30"

47°

42'

00"

EXPLANATION

Chlorinated volatile organic 

compound mass loading

Constant concentration, in percent

Dogfish Bay (100%)

Tide flats (100%)

Mud flats (50%)

Former landfill - Operable Unit 1 (OU-1)

Transport model domain

Naval Base Kitsap (NBK)

Figure 19. Domain and boundary conditions in transport simulations, Operable Unit 1, Naval Base Kitsap, Keyport, Washington.
Simulation of Variable-Density Flow and Transport of Chlorinated Ethenes 37

CVOC Simulation

Inflow of CVOCs in the former landfill at OU-1 was 

represented though mass-loading rates (MT-1) that remained 

constant throughout the dual-solute simulation. Concentration 

units of mg/L were chosen to represent CVOCs in the dualsolute simulation in order to facilitate comparison of model 

results with existing maps of the CVOCs distribution at OU-1 

that use these concentration units. Units of mass in SEAWAT 

are defined as the product of concentration units and water 

volume. If inconsistent units are specified in the model (for 

example, concentration units of mg/L and length units of 

ft), then mass units are given in the unconventional units of 

milligram-cubic feet per liter [mg-ft3/L]. A conversion factor 

of 2.832×10-5 milligram-cubic feet per liter (kg-L/mg-ft3) 

is applied to convert the mass units computed by model 

simulations to kg (Christopher Neville, S.S. Papadopulos & 

Associates, Inc., written commun, 2009).

Mass-loading rates were specified at model cells identified through analysis of CVOCs concentration data collected 

from sediment cores and groundwater wells within the former 

landfill at OU-1 from 2015 to 2017. Molar concentrations of 

trichloroethene (TCE), 1-2-cis-dichloethene (DCE) and vinyl 

chloride (VC) were summed to obtain total CVOC concentrations at different depths at seven persistent contaminant 

sources in OU-1, including model cells corresponding to 

monitoring wells MW1-50, -56, -57-, and -58 and in the southern part of the landfill, and model cells corresponding to monitoring wells MW1-46, -47-, and -48 in the central part of the 

landfill (fig. 3). Mass-loading rates were adjusted in transport 

simulations to reproduce the maximum CVOC concentrations 

observed at these seven locations.

Complete, sequential reductive dechlorination of CVOCs 

has been demonstrated at OU-1 through the identification of 

redox conditions and the presence of ethene (Dinicola, 2006), 

and mass—degradation rates of total CVOC concentrations 

have been estimated from changes in concentration data along 

groundwater flow paths. Individual degradation rates between 

pairs of parent and daughter compounds in the degradation 

chain have not been estimated, however. For this reason, the 

expected degradation of CVOCs was neglected in the initial 

transport simulation described below, with the understanding 

that a considerable portion of the CVOCs generated through 

mass loading would be biodegraded to ethene and not discharge with groundwater at downgradient boundaries, as 

predicted by the simulation. The goal of the initial transport 

simulations was to reproduce the observed configuration of the 

CVOCs plume at OU-1, and thereby estimate both the massloading rates of CVOCs from persistent contaminant sources 

and the total mass of CVOCs in the contaminant plume. 

An additional simulation described further on addresses 

the potential effects of biodegradation of CVOCs at OU-1 

using Model A.

Simulation Results

Seawater Simulation

The mass of seawater in the transport domain increases 

rapidly at the beginning of the 5,000-yr, single-solute transport 

simulation with Model A, but then increases gradually near the 

end of the simulation as it asymptotically approaches a limit 

of slightly less than 3×109 kg (fig. 20). The mass of seawater 

in the transport simulation with Model B is about 12 percent 

larger than with Model A. The increase in seawater mass at 

the end of the transport simulations results from the continued 

influx of higher-density seawater that displaces lower-density 

fresh water in low permeability sediments of the Clover Park 

Aquitard underlying Dogfish and Liberty Bays. A 3D depiction of the computed seawater distribution (fig. 21) indicates 

that seawater concentrations greater than 1 percent (about 

340 mg/L) extend 100 ft beyond the tide flats, although the 

predicted penetration of seawater is farther on the western end 

of Keyport Peninsula where the model cells are larger.

A vertical section from the marsh pond north through 

the tide flats and the western arm of the Keyport Peninsula 

shows the extent of the transition zone between fresh water 

and seawater near OU-1 with Model A (fig. 22). The transition zone is about 100 ft wide and extends to the bottom of 

the Qv sediments. Saline water overlies fresh water at shallow 

depths south of the tide flats where groundwater in the deeper, 

permeable Qv sediments discharges into seawater. The transition zone with Model B is about 25 ft wider than with Model 

A as a result of the higher K values in Units 1 and 2 (Qvr and 

Qestuary). The computed distribution of saline water extends 

into the Clover Park Aquitard under the deeper part of the tide 

flats and Liberty Bay.

Chlorinated Volatile Organic Compounds 

Simulation

The mass of CVOCs beneath OU-1 increases rapidly at 

the beginning of the 50-yr, dual-solute transport simulation 

with Model A, but then increases gradually and approaches a 

dynamic equilibrium when the contaminant mass in groundwater reaches about 256 kg (fig. 23). Almost all the contaminant mass (134 g/d) discharges to small creeks draining 

the southern edge of OU-1, while most of the remainder 

discharges to the marsh pond and creek (1.7 g/d) and the tide 

flats (49 mg/d). The CVOC mass-loading rate is 150 g/d, 

of which about 95 percent sustains the southern plume. The 

mass-loading rate is within the range of biodegradation rates 

of 130–170 g/d estimated by Dinicola (2006) using CVOC 

concentration data from 1999 to 2000 and 2004, respectively. 

The actual mass-loading rate should be greater than the 

biodegradation rate, as a portion of CVOC mass discharges to 

surface waters adjacent to OU-1. Dinicola (2006) estimated 

that this mass flux ranged from 7 to 9 g/d in 2000 and 2004, 

respectively.
38 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 4,500 5,000

Years

0.0

5.0

1.0

1.5

2.0

2.5

3.0

Seawater mass, x 109 kilograms

Model A

Model B

EXPLANATION

Figure 20. Computed mass of seawater during 5,000-year, single-solute simulations 

with models A and B, Operable Unit 1, Naval Base Kitsap, Keyport, Washington.

C

C'

Landfill

Marsh

 pond

Dogfish

Bay

Shallow

lagoon

OBLIQUE VIEW, NOT TO SCALE

VERTICAL EXAGGERATION x2

Seawater concentration, 

in percent

100

10

1

EXPLANATION

Figure 21. Computed seawater concentrations in the Qv sediments at the end of 5,000-year, single-solute simulation with Model A, 

Operable Unit 1, Naval Base Kitsap, Keyport, Washington.
Simulation of Variable-Density Flow and Transport of Chlorinated Ethenes 39

C C' Tide flats

Transition

zone

Model layers 

1 to 11

0 100 FEET

VERTICAL EXAGGERATION x2

Seawater concentration, 

in percent

100

10

1

EXPLANATION

Figure 22. Vertical section C–C’ through tide flats showing computed seawater concentrations and transition zone between fresh 

water and seawater at the end of 5,000-year, single-solute simulation with Model A, Operable Unit 1, Naval Base Kitsap, Keyport, 

Washington.

0

50

100

150

200

250

300

350

0 5 10 15 20 25 30 35 40 45 50

CVOC mass, in kilograms

Years

Total-variation-diminishing solver

Finite-difference solver

Model A

Model B

EXPLANATION

Figure 23. Computed mass of chlorinated volatile organic compounds (CVOCs) during 50-year, dual-solute 

simulations using total-variation-diminishing and implicit finite-difference solvers with Model A, and TVD 

solver with Model B, Operable Unit 1, Naval Base Kitsap, Keyport, Washington.
40 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

The mass of CVOCs beneath OU-1 is slightly greater 

with Model B (310 kg), but the CVOC mass-loading rate is 

only 73 g/d. The lower mass-loading rate results from the 

lower volume of groundwater flow in the deeper Qv sediments 

beneath OU-1 in Model B, in which the vertical anisotropy 

in Unit 5 (Qvsemi) is higher than in Model A. As a result, 

less contaminant mass (48 g/d) discharges to the small creeks 

draining the southern edge of OU-1 than in Model A. More 

CVOC mass discharges to the marsh pond and creek (7.8 g/d) 

and the tide flats (160 mg/d) in Model B, however, because the 

contaminant plumes are diverted westward to Dogfish Bay by 

the higher K values in Unit 2 (Qestuary).

Two distinct contaminant plumes emanate from the 

persistent sources in the transport simulation with Model 

A, a northern plume that migrates toward the marsh pond 

and creek, and a southern plume that migrates toward small 

creeks south of OU-1 (fig. 24). The extents of the simulated 

CVOC plumes agree favorably with those delineated from 

observed data (fig. 3), and the maximum CVOC concentrations for the northern and southern plumes (120 and 630 mg/L, 

respectively) are close to those estimated for the persistent 

contaminant sources. The southern plume accounts for about 

70 percent of the total contaminant mass. Both CVOC plumes 

extend from the water table to the bottom of the Qv sediments. 

The CVOC plumes migrate downward to the base of the Qv 

sediments and then laterally before migrating upward to discharge in the marsh and creeks.

The two contaminant plumes emanating from the 

persistent sources intermingle in the transport simulation 

with Model B and become indistinguishable downgradient 

(fig. 25). Unfortunately, no contaminant data are available 

in locations between the two distinct plumes predicted by 

Model A, so the possibility that the Model B simulation is a 

more accurate depiction of actual CVOC migration cannot be 

either verified or ruled out on the basis of available data. The 

maximum CVOC concentrations at the northern and southern persistent sources with Model B are the same as those 

predicted by Model A (120 and 630 mg/L, respectively). As 

indicated above, the combined plumes in the Model B simulation migrate toward Dogfish Bay and less contaminant mass 

discharges to the small creeks at the southern edge of OU-1.

A vertical section from the marsh pond north through 

the tide flats shows the extent of the northern CVOC plume 

beneath Marsh Creek that discharges to the tide flats in the 

Model A simulation (fig. 26). The northern plume extends 

beneath the tide flats where it mixes with higher-density seawater. The seawater presents a hydraulic barrier that limits the 

upward discharge of CVOCs to the tide flats. As a result, most 

of the CVOCs in this area are forced to discharge to Marsh 

Creek. A vertical section through the CVOC contaminant distribution simulated with Model B is similar, but the combined 

plumes underlie the tide flats to the east of the plume simulated with Model A.

Application of the FD solver in place of the TVD solver 

in the 50-yr, dual-solute transport simulation with Model 

A produces a contaminant mass (304 kg) that is about 19 

percent larger as a result of numerical dispersion (fig. 23). 

The increased dispersion generates larger CVOC plumes that 

overlap beneath the pond in the marsh (fig. 27). The difference 

in the extents of the CVOC plumes predicted by the TVD and 

FD solvers indicates that the magnitude of numerical dispersion associated with the FD solver ranges from 50 to 100 ft, or 

two to four grid-cell widths.
Simulation of Variable-Density Flow and Transport of Chlorinated Ethenes 41

Shapely

Creek

Marsh Creek

Seasonal Creek

Tide

flats

Marsh 

Pond

KEYS RD

BRADLEY RD

BROWNSVILLE HWY NE

SHARPLEY RD

NE TAGHOLM RD

NE STATE HWY 308

GADBERRY ST

TORPEDO RD

D

D'

0 125 250 FEET

0 25 50 METERS

122°37'45" 122°37'40" 122°37'35" 122°37'30"

47°

41'

55"

47°

41'

50"

47°

41'

45"

Base map modified from U.S. Geological Survey and other digital

data at 1:2,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

Former landfill - Operable Unit 1 (OU-1)

Cross section trace

Naval Base Kitsap (NBK)

Chlorinated volatile organic compound 

concentration, in milligrams per liter

0.01 to 0.5

0.51 to 1.0

1.01 to 10

10.01 to 100

100.01 to 300

300.01 to 650

A

D D'

Figure 24. Computed chlorinated volatile organic compound (CVOC) concentrations in the Qv sediments at the end of 50-year, 

dual-solute simulation with Model A: A, plan view in Unit 6 (Qvinter), B, three-dimensional perspective view.
42 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

D D' D"

Landfill

Marsh

 Pond

Northern

Southern Plume

Plume

Tide flats

OBLIQUE VIEW, NOT TO SCALE

VERTICAL EXAGGERATION x2

B

Chlorinated volatile organic 

compounds concentration, 

in milligrams per liter

1,000

10

1

0.01

EXPLANATION

Figure 24.—Continued
Simulation of Variable-Density Flow and Transport of Chlorinated Ethenes 43

Shapely

Creek

Marsh Creek

Seasonal Creek

Tide

flats

Marsh 

Pond

KEYS RD

BRADLEY RD

BROWNSVILLE HWY NE

SHARPLEY RD

NE TAGHOLM RD

NE STATE HWY 308

GADBERRY ST

TORPEDO RD

0 125 250 FEET

0 25 50 METERS

122°37'45" 122°37'40" 122°37'35" 122°37'30"

47°

41'

55"

47°

41'

50"

47°

41'

45"

Base map modified from U.S. Geological Survey and other digital

data at 1:2,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

Former landfill - Operable Unit 1 (OU-1)

Naval Base Kitsap (NBK)

Chlorinated volatile organic compound 

concentration, in milligrams per liter

0.01 to 0.5

0.51 to 1.0

1.01 to 10

10.01 to 100

100.01 to 300

300.01 to 650

A

Figure 25. Computed chlorinated volatile organic compound (CVOC) concentrations in the Qv sediments at the end of 50-year, 

dual-solute simulation with Model B: A, plan view in Unit 6 (Qvinter), B, three-dimensional perspective view.
44 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

D D' D"

Landfill

Marsh

 Pond

Northern

plume Southern

plume

Tide flats

OBLIQUE VIEW, NOT TO SCALE

VERTICAL EXAGGERATION x2

B

Chlorinated volatile organic 

compounds concentration,

in milligrams per liter

1,000

10

1

0.01

EXPLANATION

Figure 25.—Continued

D'

D D'

D"

Tide flats

Model layers 

1 to 11

0 100 FEET

VERTICAL EXAGGERATION x2 Chlorinated volatile organic 

compounds concentration, 

in milligrams per liter

1,000

10

1

0.01

Seawater concentration, 

in percent

100

10

1

EXPLANATION

Figure 26. Vertical section from marsh pond through tide flats showing computed chlorinated volatile organic compound (CVOC) and 

seawater concentrations at the end of 50-year, dual-solute simulation with Model A.
Simulation of Variable-Density Flow and Transport of Chlorinated Ethenes 45

Shapely

Creek

Marsh Creek

Seasonal Creek

Tide

flats

Marsh 

Pond

KEYS RD

BRADLEY RD

BROWNSVILLE HWY NE

SHARPLEY RD

NE TAGHOLM RD

NE STATE HWY 308

GADBERRY ST

TORPEDO RD

0 125 250 FEET

0 25 50 METERS

122°37'45" 122°37'40" 122°37'35" 122°37'30"

47°

41'

55"

47°

41'

50"

47°

41'

45"

Base map modified from U.S. Geological Survey and other digital

data at 1:2,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

Former landfill - Operable Unit 1 (OU-1)

Numerical dispersion with finite-difference 

solver

Naval Base Kitsap (NBK)

Chlorinated volatile organic compounds 

concentration, in milligrams per liter

0.01 to 0.5

0.51 to 1.0

1.01 to 10

10.01 to 100

100.01 to 300

300.01 to 650

Figure 27. Comparison of the extents of chlorinated volatile organic compound (CVOC) migration with total-variation-diminishing 

and implicit finite-difference solvers with Model A showing effect of numerical dispersion.
46 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

Discussion of Simulation Results

Containment of Chlorinated Volatile Organic 

Compounds Contaminants

Results of the groundwater-flow and contaminanttransport simulations indicate that both geologic and hydraulic barriers limit the migration of CVOCs from the former 

landfill at OU-1. The CVOC contaminant plumes at OU-1 in 

both models A and B are largely contained within Naval Base 

Kitsap, although a small quantity of CVOCs discharges to the 

tide flats. The Clover Park Aquitard, a low-permeability geologic barrier that underlies the Keyport Peninsula and adjacent 

uplands, limits vertical groundwater flow between shallow, 

contaminated Qv sediments at OU-1 and deeper geologic 

units. Moreover, evidence from measured groundwater levels 

and simulated groundwater flow in both the regional Kitsap 

Model And Keyport model indicate that the hydraulic gradient is upward beneath the Keyport Peninsula, and presents an 

additional hydraulic barrier to the potential downward migration of CVOCs from OU-1 through the Clover Park Aquitard.

Underflow from upland areas south of the Keyport 

Peninsula discharges to the marsh at the base of the peninsula 

and presents an additional hydraulic barrier, limiting CVOC 

migration to the south. Contaminant transport simulations 

indicate that the CVOC plumes discharge mainly to the marsh 

pond and creeks. Only a very small portion of the simulated 

contaminant mass (less than 0.1 percent in Model A and 0.3 

percent in Model B) discharges into seawater in the tide flats. 

It should be noted that neither of these simulations include the 

effects of biodegradation, which would decrease the predicted 

CVOC mass reaching the tide flats. The saline water that 

underlies seawater in the tide flats and in Dogfish and Liberty 

Bays west of OU-1 has a higher density than fresh groundwater and limits upward discharge of CVOCs into seawater. This 

saline water presents an additional hydraulic barrier that prevents westward CVOC migration beneath Dogfish Bay toward 

residential wells on Virginia Point.

Model Uncertainty

There are no shallow pumped wells on the Keyport 

Peninsula and, consequently, no data are available concerning the rate of groundwater flow through the Qv sediments 

in the area surrounding OU-1. Recharge to groundwater has 

been estimated by the SWB model described earlier, but the 

SWB recharge rates were adjusted downward during calibration of both models A and B to better match the observed 

base flow in the perennial stream that drains to the shallow 

lagoon, and to maintain the simulated water table below land 

surface. Estimated recharge rates are correlated with estimated 

hydraulic conductivities for hydrogeologic units in the Qv 

sediments, however. A combination of higher recharge rates 

and larger hydraulic conductivities would therefore yield 

higher groundwater flow rates with little change in model 

error. Uncertainty in the rates of recharge and groundwater 

flow produces uncertainty in the estimated mass-loading rate 

of CVOCs.

At higher rates of recharge and groundwater flow, the 

CVOCs entering the flow system are transported faster to 

discharge boundaries, requiring larger mass-loading rates to 

maintain a stable configuration of the contaminant plume. The 

CVOC mass-loading rate with Model B (73 g/d) was lower, 

however, because most of the increased recharge in that model 

flows through the upper part of the Qv sediments, due to the 

increased vertical anisotropy in Unit 5 (Qvsemi). Data on the 

rate of groundwater flow on the Keyport Peninsula and under 

OU-1, collected through cross-hole aquifer tests or additional 

measurements of base flow in perennial streams, could be used 

to narrow the 73–240 g/d range of uncertainty associated with 

the estimated CVOC mass-loading rate at the OU-1 site.

Neither Model A nor B predicted groundwater flow or 

contaminant transport from OU-1 to wells MW1-38 and -39 

located on the State Route 308 causeway to the northwest of 

the tide flats because of the hydraulic barriers identified above. 

The contaminants 1,4-dioxane and vinyl chloride have been 

detected in one or both of those wells at low concentrations 

(maximum values of 4.1 and 2.7 μg/l). Dinicola and others 

(2002) could not determine whether groundwater at these 

causeway wells originated from OU-1, based on multiple lines 

of evidence from chemical isotope and common-ion geochemical data. There are no wells in the tide flats to provide data 

on the continuity of hydrogeologic units between the landfill 

and the causeway, so there is uncertainty as to whether the 

contaminants measured in those wells originated in the landfill 

or elsewhere.

Effect of Biodegradation

Previous studies (Dinicola and others, 2002; Dinicola, 

2006) have established that reductive dechlorination of 

CVOCs is occurring at OU-1, and Dinicola (2006) estimated 

CVOC biodegradation rates ranging from 130 to 170 g/d. A 

first-order, removal rate-constant λ for CVOCs at OU-1 was 

computed from these rate estimates using the CVOC mass in 

the contaminant plume obtained from the 50-yr, dual-solute 

transport simulation with Model A. The rate constant λ (d-1) is 

related to the reaction half-life t1/2 as:

λ =

ln 2 _

t_1

2

 (2)

where

t1/2 is the time required for the CVOC mass to 

decrease to half the original mass under the 

influence of biodegradation alone.

The half-life of the 256-kg contaminant plume predicted 

by the simulation is about 850 days at a biodegradation rate of 

150 g/d, which is equivalent to a rate constant of about 8×10-4

d-1. The comparable half-life and rate constant computed with 
Discussion of Simulation Results 47

the CVOC mass predicted with transport simulation with 

Model B are 1,000 days and 7×10-4 d-1. These estimated halflives are within the range of t1/2 values (35 to 3,600 days) for 

TCE, cis-DCE and VC reported for field sites by Wilson and 

others (1996).

An additional 50-yr, dual solute contaminant transport 

simulation with Model A was conducted in which CVOC 

biodegradation was represented as a first-order reaction with 

a rate constant of 10-3 d-1. In this simulation the CVOC plume 

reached a stable mass of 77 kg within 10 years of migration, although the CVOC mass-loading rates and maximum 

CVOC concentrations were equivalent to those in the original contaminant-transport simulation described previously 

(fig. 28). About one half the CVOC mass was lost through 

biodegradation at a rate of 75 g/d, and nearly all the remainder (71 g/d) was discharged by the southern plume to a small 

creek south of OU-1. The CVOC mass discharged from the 

northern plume to the marsh and tide flats was 0.95 mg/d. The 

areas of the northern and southern plumes (fig. 29) are 67 percent and 32 percent less, respectively, than the areas computed 

by the original transport simulation (fig. 24), and the northern 

plume does not extend past the landfill boundary at concentrations exceeding 10 μg/L.

Refined estimates of CVOC mass-loading rates and 

biodegradation rates could be obtained in future transport 

simulations by matching CVOC concentrations at different 

locations within the contaminant plumes. Alternatively, individual biodegradation rates for components in the reductive 

dechlorination chain could be estimated by matching TCE, 

cis-DCE, and VC concentrations in a transport simulation that 

represents coupled sequential reactions. These simulations 

would also better estimate the mass-loading rates of CVOCs 

at OU-1. This effort was beyond the scope of the current study 

and would necessitate another variable-density, contaminanttransport program as SEAWAT Version-4 does not support 

coupled sequential reactions.

0

50

100

150

200

250

300

0 5 10 15 20 25 30 35 40 45 50

Chlorinated volatile organic compound mass, in kilograms

Years

Model A

With biodegradation

EXPLANATION

Figure 28. Computed mass of chlorinated volatile organic compounds (CVOCs) during original 50-year, 

dual-solute simulation with Model A, and alternate simulation with biodegradation.
48 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

Shapely

Creek

Marsh Creek

Seasonal Creek

Tide

flats

Marsh 

Pond

KEYS RD

BRADLEY RD

BROWNSVILLE HWY NE

SHARPLEY RD

NE TAGHOLM RD

NE STATE HWY 308

GADBERRY ST

TORPEDO RD

0 125 250 FEET

0 25 50 METERS

122°37'45" 122°37'40" 122°37'35" 122°37'30"

47°

41'

55"

47°

41'

50"

47°

41'

45"

Base map modified from U.S. Geological Survey and other digital

data at 1:2,000. Coordinate Reference System State Plane

Washington North. Horizontal datum is North American Datum

of 1983

EXPLANATION

Former landfill - Operable Unit 1 (OU-1)

Naval Base Kitsap (NBK)

Chlorinated volatile organic compounds 

concentration, in milligrams per liter

0.01 to 0.5

0.51 to 1.0

1.01 to 10

10.01 to 100

100.01 to 300

300.01 to 650

A

Figure 29. Computed chlorinated volatile organic compound (CVOC) concentrations in the Qv sediments at the end of 50-year, 

dual-solute simulation with Model A including biodegradation: A, plan view, B, three-dimensional perspective view.
Discussion of Simulation Results 49

Alternative Models

Two models (A and B) of groundwater flow and contaminant transport are presented that each provide an adequate 

fit to hydraulic and contaminant data. As indicated above, 

the determination as to which model more accurately predicts actual contaminant migration at OU-1 cannot be made 

without additional contaminant data. Both models predict 

similar contaminant masses and biodegradation rate constants 

and also predict that contaminant migration is limited by 

the same geologic and hydraulic barriers. Model predictions 

differ, however, in terms of the hydraulic parameters used in 

model simulations (recharge rates, and horizontal and vertical 

hydraulic conductivities) and in the CVOC mass-loading rates 

and contaminant plume configurations. The groundwater flow 

system simulated with Model B better represents the effects 

of tidal fluctuations on unconfined storage than does Model A 

and also uses recharge rates that are closer to those estimated 

with the SWB Model And predicts water table altitudes near 

OU-1 that are below land surface. Model B overestimates the 

effects of unconfined storage, however, and underpredicts the 

response of groundwater levels to tidal fluctuations and also 

overpredicts measured stream flow. In the following sections, Model A, which better represents the groundwater-level 

response to tidal fluctuations, is used to demonstrate the effects 

of the representation of tidal fluctuations and potential remedial activities. The effects illustrated using Model A would be 

similar if Model B were chosen instead, however.

A version of MODFLOW 6 will soon be available that 

supports variable-density flow and contaminant transport 

(C.D. Langevin, U.S. Geological Survey, written commun., 

2019). Application of this MODFLOW 6 version could represent the distribution of discontinuous hydrogeologic units in 

the Keyport model domain directly (without using the HUF2 

package), permit the use of the NWT solver and, therefore, 

better represent unconfined and confined storage. A single 

model could then be developed to simulate groundwater flow 

and CVOC contaminant migration at OU-1 and greatly simplify the modeling approach described in the current study.

Effect of Tidal Fluctuations

The 50-yr, steady-state, contaminant-transport simulations of CVOC migration with Model A discussed previously represent the mean (average of high and low tide) tidal 

conditions in Dogfish and Liberty Bays and Port Orchard 

based on predicted high and low tides in 2017. The depiction 

of mean tidal conditions was necessary in order to complete 

the simulations within a reasonable time, as determined by 

the small time increments required by the TVD solver. Two 

alternate simulations with Model A were conducted to assess 

the sensitivity of predicted CVOC migration from OU-1 to the 

depiction of tidal fluctuations. In the first simulation (C), the 

mean tide was computed for Dogfish and Liberty Bays, Port 

Orchard, and the tide flats and Marsh Creek, as the average of 

Landfill

Marsh

 Pond

Northern

Southern plume

plume

Tide flats

OBLIQUE VIEW, NOT TO SCALE.

VERTICAL EXAGGERATION x2

B

Chlorinated volatile organic 

compounds concentration, 

in milligrams per liter

1,000

10

1

0.01

EXPLANATION

Figure 29.—Continued
50 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

water levels measured during the 31-day period in July and 

August 2018. The mean water levels in steady-state simulation 

A were lowered by 1.2 ft in Dogfish and Liberty Bays, 0.87 

ft in the tide flats, and 0.1 ft in Marsh Creek. In the second 

simulation (D) the FD solver was applied in a 51-yr transient 

simulation that consisted of 600 repetitions of the 31.05-day 

transient periods used for calibration of the groundwaterflow model. In transient simulation B, tidal fluctuations were 

represented as shown in figure 11. The initial conditions for 

the simulation D were the same as those used in 50-yr, steadystate transport simulation computed with the TVD solver with 

Model A and required several days to complete.

The results of the steady-state TVD simulation C and 

the transient FD simulation D are compared with the original 

steady-state TVD and FD simulations with Model A depicted 

previously in figure 24. Under steady-state simulation C the 

CVOC mass in the contaminant plumes increased slightly 

faster than in the original TVD simulation, and the CVOC 

mass after 50 years (272 kg) was 6 percent higher (fig. 30). 

The simulated northern CVOC plume extended about 100 ft 

farther under the tide flats, and CVOC discharge to the tide 

flats increased from 49 mg/d under the original TVD simulation to 890 mg/d. CVOC discharge to Marsh Creek and pond 

were reduced by 46 percent from 1,700 to 770 mg/d. Under 

transient simulation D the CVOC mass in the contaminant 

plumes increased slightly slower than in the original FD 

simulation, and the CVOC mass after 50 years (297 kg) was 

3 percent less (fig. 30). The simulated northern CVOC plume 

extended about 75 ft farther under the tide flats, and CVOC 

discharge to the tide flats increased from 67 mg/d under the 

original FD simulation to 75 mg/d. The extent of the transition 

zone between fresh water and seawater north of the marsh was 

increased from 100 ft in the original steady-state FD simulation to 125 ft in transient simulation B.

The results of these pairs of simulations indicate that 

the representation of tidal fluctuations has a small effect 

on the sizes and masses of the simulated CVOC plumes. 

Comparisons of the TVD and FD simulations must be considered separately because the effects of numerical dispersion are 

greater than the effects of tidal fluctuations on the simulated 

CVOC plumes. Representing the mean tidal stages in steadystate simulation C as the averages of water levels measured 

during the 31-day period in July and August 2018 predicts 

larger CVOC plumes than in the original TVD simulation, 

in which the mean tidal stages computed with the average of 

high and low tide stages are higher. The larger northern CVOC 

plume results in more CVOC discharge to the tide flats and 

less discharge to the marsh. Data required to compute the averages of water levels in all surface-water bodies are only available for the 31-day period in July and August 2018, however, 

which may not be representative of annual tidal fluctuations.

In contrast, representation of tidal fluctuations in transient simulation D results in smaller CVOC plumes and only 

slightly more CVOC discharge to the tide flats than in the 

original steady-state FD simulation. Simulated tidal stages 

in transient simulation D are more frequently higher than the 

mean stage, despite the representation of the low tide period. 

In both pairs of simulations, the increase in CVOC discharge 

from the northern CVOC plume to the tide flats is offset by a 

decrease in CVOC discharge to the marsh. It is important to 

note that neither of the pairs of simulations represents biodegradation, which was previously shown to reduce CVOC 

discharge to the tide flats and the marsh by a factor of 50.

The small effect of tidal fluctuations on the sizes and 

masses of the simulated CVOC plumes is expected because 

the period of the fluctuations (12 hr) is rapid compared with 

the groundwater velocity (20 ft/yr) beneath OU-1. The rapid 

tidal fluctuations slightly alter groundwater-flow paths and 

the representation of lower tidal stages causes the simulated 

northern CVOC plume to lengthen slightly under the tide 

flats. These results are consistent with those of La Licata and 

others (2007), who simulated the effect of tidal fluctuations on 

variable-density groundwater flow and contaminant transport 

in a coastal aquifer in Italy. They found that representing tidal 

fluctuations had little effect on the contaminant mass discharged to seawater. The principal effect of tidal fluctuations 

in their simulations was to increase the width of the transition 

zone between fresh water and seawater, thereby increasing dispersion of contaminants near the seashore. An increase in the 

transition zone was also noted in transient simulation D, indicating that tidal fluctuations could affect CVOC concentrations 

near the tide flats. These effects are likely small compared with 

the larger effects produced by biodegradation, however.

Implications for Remediation of Contaminant 

Sources

The U.S. Navy has conducted a comprehensive characterization program at the OU-1 site to delineate the locations of persistent contaminant sources in the former landfill 

(Battelle Memorial Institute, 2018). Current studies at the site 

are focused on potential remediation strategies to attenuate 

or isolate these contaminant sources. Three additional 50-yr, 

steady-state contaminant-transport simulations with Model A 

were conducted as a complement to these studies to predict the 

time required for potential source-control measures to diminish the mass of CVOCs in groundwater at OU-1. The initial 

conditions for the contaminant-transport simulations were 

taken from the end of the original 50-yr, steady-state simulation of CVOC migration with Model A without biodegradation 

discussed above. The mass-loading rates of CVOCs in the 

three simulations were reduced by factors of 50, 90 and 99 

percent, respectively.

In each of the three contaminant source-control scenarios 

the mass of CVOCs decreased rapidly at first and then slowly 

approaches a new dynamic equilibrium (fig. 31). The reduction in CVOC mass in groundwater after 50 years in each 

scenario is slightly less than the magnitude of source reduction 

(table 6). The reductions of CVOC mass in groundwater after 

50 years for the three scenarios are 46, 86, and 95 percent, 

respectively. The sizes of the contaminant plumes are largely 

unchanged from the original transport simulation (fig. 24), 

although the CVOC concentrations in each are reduced.
Discussion of Simulation Results 51

0

50

100

150

200

250

300

350

0 5 10 15 20 25 30 35 40 45 50

Chlorinated volatile organic compound mass,

in kilograms

Years

Finite-difference solver, Model A

Total-variation-diminishing solver, Model A

Finite-difference solver, simulation D

Total-variation-diminishing solver, simulation C

EXPLANATION

Figure 30. Computed mass of chlorinated volatile organic compounds (CVOCs) during original 50-year, 

dual-solute simulations using total-variation-diminishing and implicit finite-difference solvers with Model A, and 

alternate steady-state simulation C with lower tidal stages and alternate transient simulation D with transient 

representation of tidal fluctuations.

EXPLANATION

0

50

100

150

200

250

300

0 5 10 15 20 25 30 35 40 45 50

Mass, in kilograms

Years

90

99

50

Reduction in mass-loading rates,

in percent

Figure 31. Reduction in chlorinated volatile organic compound (CVOC) mass with time resulting from three 

contaminant source-control scenarios with Model A that reduce CVOC mass-loading rates.
52 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

Summary

Chlorinated volatile organic compounds (CVOCs) have 

migrated to groundwater beneath a former 9-acre landfill 

at Operable Unit 1 (OU-1) of the Naval Base Kitsap on the 

Keyport Peninsula, in Kitsap County, Washington. Previous 

studies demonstrated that biodegradation of CVOCs in shallow groundwater at OU-1 is substantial and prevents most 

of the mass of dissolved-phase CVOCs in groundwater from 

discharging to surface water, but contaminant concentrations 

up to 630 milligrams per liter persist in localized areas, likely 

from the dissolution of residual, non-aqueous phase liquids. 

Variable-density groundwater-flow and contaminant-transport 

models were developed using SEAWAT-Version 4 to simulate 

the direction and rate of groundwater flow in a 5.9 mi2 area 

surrounding the Keyport Peninsula, estimate the CVOC mass 

in groundwater and the rate of mass loading, and to assess 

possible remedial activities at OU-1.

The landfill at OU-1 is on a narrow strip of land connecting the Keyport Peninsula to the mainland and is adjacent to 

tide flats extending from Dogfish Bay, a segment of Liberty 

Bay connected to Puget Sound. The study area is underlain by 

Quaternary deposits consisting of alternating glacial and interglacial sediments from 500 to 1,500 ft thick. A sequence of six 

units was delineated in a hydrogeologic model within a relatively thin package (less than 100 ft) of recent sediments from 

the Vashon Stade (Qv) and younger deposition beneath the 

Keyport Peninsula. These Qv sediments are underlain by the 

much thicker (more than 300 ft) Clover Park Aquitard (QC1), 

which in turn overlies a confined, sea-level aquifer (QA1).

A shallow lagoon, tide flats, marsh (including a pond 

and outlet creek), and Dogfish Bay are the principal areas 

of groundwater discharge near OU-1. Water levels in these 

surface-water bodies fluctuate in response to tidal cycles, but 

manmade structures limit the extent of water-level fluctuations in the tide flats, marsh, and shallow lagoon. Groundwater 

flow at OU-1 is derived from lateral flows from the Keyport 

Peninsula to the northeast and from the Kitsap Peninsula 

uplands to the southwest and from vertical flow from the 

underlying sea-level aquifer (QA1).

Daily tidal fluctuations in Liberty Bay follow a semidiurnal pattern that includes two high tides and two low tides 

each day, with a higher and lower high tide and a higher and 

lower low tide. The mean high-high tide in Liberty Bay in 

2017 was 7.94 ft, and the mean low-low tide was -2.54 ft 

(based on tides predicted by NOAA at Poulsbo, Washington), 

resulting in a mean range of 10.53 ft. The daily tidal cycle 

had a mean length of 24.84 hr. The daily range in water-level 

fluctuations measured in July and August 2018 in the tide flats, 

marsh outlet creek and marsh pond were 8, 1.5, and 0.5 ft, 

respectively. Groundwater-level response to tidal fluctuations 

in monitoring wells in and near OU-1 was as much as 7.3 ft. 

The hydraulic diffusivity was estimated as 5.4×105 ft2/d from 

the attenuation of tidal fluctuations measured in 12 monitoring wells.

The magnitude and distribution of groundwater recharge 

in the 5.9 mi2 study area was computed using national spatial 

datasets and the SWB model using a modified Thornthwaite 

Mather, soil-water-balance approach at daily time steps for the 

period 1980 to 2015. The daily recharge rates were calculated 

on a 25-ft grid, summed and averaged to obtain spatially 

distributed, mean-annual recharge rates that reflect the spatial 

distribution of land cover and soil properties in the study area. 

The average recharge rate computed for the study area was 

16.7 in/yr.

Constant-Density Groundwater-Flow Simulation

Three-dimensional (3D), constant-density groundwater flow in the 5.9 mi2 study area was simulated using 

MODFLOW-2005 and MODFLOW-NWT models run 

alternately in an iterative sequence: a combined steady-state 

and transient simulation using MODFLOW-2005 followed 

by a steady-state simulation using MODFLOW-NWT. The 

HUF2 package was employed to represent the discontinuous hydrogeologic units beneath OU-1 in MODFLOW-2005, 

a program that requires the aquifer system to be treated as 

confined to obtain a numerical solution for hydraulic head in 

the study area, where steep slopes in the uplands transition 

abruptly to the level topography of the lowlands. In contrast, 

the MODFLOW-NWT program (which does not support the 

HUF2 package) can better simulate the transition in hydraulic head from the uplands to the lowlands while treating the 

aquifer system as unconfined and computes underflow from 

the uplands using a more realistic value of saturated thickness 

than does MODFLOW-2005. The water table computed by the 

MODFLOW-NWT simulation was used to specify as inactive 

the model cells in the MODFLOW-2005 simulation that were 

above the simulated water table, thereby reducing the simulated saturated thickness in the confined simulation.

Table 6. Reduction in chlorinated volatile organic compounds 

mass after 50 years of contaminant source-control measures 

with Model A, Operable Unit 1, Naval Base Kitsap, Keyport, 

Washington.

[Abbreviations: g/d, gram per day; CVOC, chlorinated volatile organic 

compounds]

Scenario

Percentage

reduction in

mass-loading 

rate

Massloading

rate 

(g/d)

Percentage 

reduction in 

CVOC mass in 

plumes

1 50 74 46

2 90 24 86

3 99 2.4 96
Summary 53

The 3D groundwater-flow model encompasses a 5.9 

mi2 area that is divided into a variably spaced grid with 15 

model layers, 120 rows, and 88 columns; cell widths ranged 

from 25 ft near OU-1 to 500 ft near model boundaries. The 

bottom 4 model layers correspond to 4 hydrogeologic units 

from the Clover Park Aquitard (Unit 7; QC1) to the sea-level 

aquifer (Unit 10; QA1), while 10 model layers represent the 

Quaternary Vashon Stade (Qv) or younger deposits. These 

Qv sediments were subdivided into 10 model layers to better 

represent flow through the surficial sediments beneath OU-1. 

The top model layer represents surface water. The thicknesses 

of the 10 model layers representing Qv sediments are equal 

within each vertical stack of model cells, and the mapping of 

hydraulic properties of the hydrogeologic units to model layers 

was computed using the HUF2 package.

The 3D groundwater-flow simulation included a steadystate period that provided initial conditions to a 31-day 

transient period that represented 30 tidal cycles. Tidal fluctuations in Dogfish and Liberty Bays and Port Orchard were 

approximated by a trapezoidal function with a period of 1.035 

days that ranged from a single low-tide stage to a high-tide 

stage (the respective means of the low-low and high-high 

tides predicted for Liberty Bay at Poulsbo, Washington, in 

2017), corrected for the specific gravity of seawater. The mean 

tide predicted for 2017 was specified for the model boundary 

during the steady-state period. Water levels specified in the 

tide flats were the same as in Dogfish Bay, with the exception 

of low tide when the water level was specified as the elevation of the sill connecting the two water bodies (0.3 ft). Water 

levels in the marsh pond and its outlet were specified as the 

average high, low and mean levels measured during July and 

August 2018.

Keyport model boundaries coincide with the model 

grid of the 3D Kitsap groundwater-flow model of Frans and 

Olsen (2016), and hydraulic heads and flow rates computed 

by the Kitsap model were used to set head and flow boundaries in the Keyport model. Hydraulic heads were specified at 

model boundaries in model layers 13 and 15 representing the 

permeable interbeds within the Clover Park Aquitard (Unit 8; 

QC1pi) and the sea-level aquifer (Unit 10; QA1), respectively, 

and vertical flows exchanged between the sea-level aquifer 

and deeper model layers in the Kitsap model were specified in 

model layer 15. Hydraulic heads were specified along southern 

and western boundaries in a single model layer within the Qv 

sediments. The elevation and head were computed as the top 

of the Clover Park Aquitard (Unit 7; QC1) plus one-third of 

the thickness of the Qv sediments to maintain the simulated 

water table below land surface. Hydraulic heads corresponding 

to the high, low and mean tide stages were also specified in 

model layer 1 along the northern and eastern boundaries where 

the model domain intersects Liberty Bay and Port Orchard.

The Drain (DRN) package was used to represent groundwater discharge to surface water at land surface, including 

perennial streams, seeps in the uplands and wet areas in the 

lowlands. Groundwater withdrawals from wells in model 

layers 15 (Unit 10; QA1) and 13 (Unit 8; QC1pi) were 

represented with the Well (WEL) package. Groundwater 

withdrawals from shallow domestic wells screened in Qv 

sediments and infiltration from domestic septic systems were 

accounted for by representing the net outflow (withdrawals 

minus infiltration) in the estimated recharge values. Recharge 

computed by the SWB model was applied to model layer 2 

with the Recharge (RCH) package, and the model domain was 

divided into five areas to adjust the recharge values during 

model calibration, including a coarse-grained area, three finegrained (Unit 3; Qvt) areas, and an area where the Clover Park 

Aquitard (Unit 7; QC1) is present at land surface.

Two alternative models (A and B) with differing storage properties were constructed to assess the ramifications 

of representing changes in storage from the fluctuating water 

table near OU-1 with a confined numerical model. An Ss

value of 2.0×10-6 ft-1 was assigned to all hydrogeologic units 

in Model A, while in Model B, the Ss value of the uppermost 

Unit 1 (Qvr sediments) was increased to 0.1 to approximate 

an unconfined storage value associated with specific yield. A 

multiplicative scaling factor based on the Qvr unit thickness 

was applied in Model B to reduce the Ss values where Unit 1 

spanned multiple model layers. Unconfined storage near the 

water table at OU-1 was underestimated in Model A and overestimated in Model B.

Both groundwater-flow models were calibrated by 

adjusting values of hydraulic properties and recharge to 

match 99 observations of four types of data: groundwater 

levels (69), water-level fluctuations in wells in response to 

tides (26), vertical hydraulic gradients (3), and streamflow 

(1). In addition, pseudo water-table altitudes computed by 

the MODFLOW-NWT simulation were specified as observations at 24 locations in the uplands. A subset of 15 of the 

24 parameter values specified in the model were estimated 

through nonlinear regression using UCODE, while values for 

the remaining 9 parameters were assigned by trial and error 

or taken from the literature. Weights assigned to observations 

included in the regression were chosen to account for differences in units or ranges of measurements and adjusted such 

that the observations were weighted equally. Parameter values 

were initially estimated in a confined transient simulation 

using MODFLOW-2005 and then updated in an unconfined 

steady-state simulation using MODFLOW-NWT. The water 

table computed by the MODFLOW-NWT simulation was then 

used to update inactive model cells in the MODFLOW-2005 

simulation, and parameter estimation was then reapplied to the 

MODFLOW-2005 simulation. This procedure was continued 

until the volumetric water budgets of both simulations were 

nearly equal.

Coefficients of variation (CV) for 8 of 15 estimated 

parameters were less than 100 percent in Model A, indicating 

that the regression was relatively sensitive to these parameters. Estimated hydraulic conductivity values for the coarsegrained, permeable Qv sediments (Units 1, 4, and 6; Qvr, Qva 

and Qvinter) ranged from 0.9 to 1.7 ft/d, slightly lower than 

the geometric mean values (3.3–4.1 ft/d) computed from 25 

slug tests. Estimated hydraulic conductivity values for Unit 2 
54 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

(Qmarsh) and Unit 3 (Qvt) were one order of magnitude lower 

than for the permeable units, and the estimated hydraulic 

conductivity for Unit 2 (Qestuary) was 1 ft/d. The hydraulic 

conductivity in the till recharge area delineated around a small 

knoll in the lower part of the uplands southwest of the Keyport 

Peninsula was estimated as 10 ft/d, indicating that more 

permeable units occur in this area than those mapped by the 

hydrogeologic model.

SWB recharge rates were reduced during nonlinear 

regression by factors of 0.65 for coarse-grained sediments, 

0.22 and 0.17 for Unit 3 (Qvt) in the uplands and lowlands, 

respectively, and 0.1 for the small knoll in the uplands. The 

average recharge rate in the model domain was 4.2 in/yr, 

much less than the average rate of 16.7 in/yr computed with 

the SWB Model But comparable to recharge rates (4 in/yr) 

estimated for the northern part of the Kitsap groundwaterflow model domain. The lower recharge rates estimated by 

nonlinear regression suggest that the SWB model overestimates recharge to the groundwater in uplands where relatively 

permeable soils at land surface are underlain by much lower 

permeability sediments at depth. In these areas infiltration 

through the soil profile likely discharges as interflow to local 

drainage, rather than to the groundwater-flow system.

The CVs for 8 of the 16 parameters estimated with Model 

B were less than 100 percent, but half of these parameters had 

CV values greater than or equal to 90 percent. Estimated K

values for Units 1, 2, and 6 (Qvr, Qestuary and Qvinter) were 

higher than in Model A (ranging from 3.0 to 19 ft/d), while K

values Unit 3 (Qvt) were an order of magnitude less. Values 

for vertical anisotropy were estimated for Units 2, 3, and 5 

(Qm, Qvt and Qvsemi) and were higher for Units 3 and 5, 

and lower for Unit 2 than in Model A. The average recharge 

rate in Model B was 7.3 in/yr, slightly less than half the SWB 

recharge rate, and 60 percent higher than in Model A.

The standard error in model fit in Model A is 7.8 ft for 

groundwater levels (13 percent of the 60-ft measurement 

range), 1.0 ft for water-level fluctuations (14 percent of the 

7.3-ft measurement range), and 11.3 ft for pseudo water-table 

altitudes (4 percent of the 320-ft measurement range). The 

corresponding numbers for Model B are 7.6 ft and 11.2 ft for 

groundwater levels and pseudo water-table altitudes and 1.3 ft 

in water-level fluctuations because the higher storage values 

dampen the response to tidal fluctuations. Residual plots 

indicate little bias in water-level fluctuations in Model A, but 

offsetting biases in water levels (overestimated) and pseudo 

water-table altitudes (underestimated). Biases in Model B are 

similar, but there is a distinct bias towards under estimating 

water-level fluctuations The directions of the three vertical 

hydraulic-gradient observations were simulated correctly in 

both models, with errors ranging from 25 to 66 percent, while 

the simulated streamflow for a perennial stream (0.18 ft3/s) 

was 70 percent less than the measured value in Model A and 

40 percent less in Model B.

The hydraulic-head distribution simulated by both 

groundwater-flow models for the steady-state period indicates 

that OU-1 is just north of a regional discharge area that limits 

groundwater flow from the landfill toward uplands south and 

west of the tide flats and the marsh, and that the Clover Park 

Aquitard (Unit 7, QC1) is a hydraulic barrier limiting flow 

from the Qv sediments to and from the deeper hydrogeologic 

units. Vertical groundwater flow is generally at least an order 

of magnitude less than horizontal flow and downward beneath 

the uplands and upward beneath Keyport Peninsula, Dogfish 

and Liberty Bays, and Port Orchard.

Variable-Density Groundwater Flow and 

Contaminant-Transport Simulations

Two SEAWAT simulations of steady-state, variabledensity flow and transient solute transport were conducted 

for models A and B in a 1.3 mi2 transport domain surrounding the Keyport Peninsula. The first 5,000-yr simulation used 

a single solute (seawater) to establish the initial distribution 

of seawater in the second 50-yr simulation that consisted of 

two solutes, seawater and CVOCs. Both simulations used 

the same model design and parameter distributions as in the 

MODFLOW-2005 groundwater-flow model. The transportmodel domain extends from surface water (model layer 1) 

to the upper part of the Clover Park Aquitard (Unit 7, QC1, 

model layer 12). A porosity value of 0.2 was specified for all 

hydrogeologic units. Longitudinal dispersivity was spatially distributed and proportional to the grid cell size (with 

a maximum value of 50 ft), and horizontal-transverse and 

vertical-transverse dispersivities were specified as ratios to the 

longitudinal dispersivity (0.1 and 0.01, respectively).

Seawater concentrations in both single- and dual-solute 

simulations ranged from 0 percent for fresh water (specific 

gravity 1.0) to 100 percent for seawater (specific gravity 

1.022) in surface-water bodies. Constant-concentration boundaries of 100 percent were specified in model layer 1 to represent Dogfish and Liberty Bays, the tide flats and Port Orchard. 

An implicit finite-difference (FD) method with upstream 

weighting was used to solve the advection-dispersion equation in the single-solute simulation to allow the use of large 

time steps (up to 100 days) in the 5,000-year simulation. The 

mass of seawater increased rapidly at the beginning of both 

single-solute transport simulations, but then gradually near 

the end and asymptotically approached limits of 2.8×109 kg 

in Model A and 3.1×109 kg in Model B. The transition zone 

between fresh water and seawater near the tide flats west of 

OU-1 was about 100 ft wide in Model A and 125 ft wide in 

Model B and extended to the bottom of the Qv sediments. 

Saline water overlies fresh water at shallow depths in this area 

where groundwater in the deeper, permeable Qv sediments 

discharges into seawater.

Inflow of CVOCs within the former landfill at OU-1 was 

represented in the dual-solute simulations by using massloading rates (MT-1) to depict total CVOC concentrations 

(TCE, cis-DCE, and VC) at different depths at seven persistent 

contaminant sources. The mass-loading rates were adjusted to 

reproduce the maximum observed CVOC concentrations in 
Summary 55

2017. The expected biodegradation of CVOCs was initially 

neglected in the transport simulations, with the understanding that a considerable portion of the CVOCs generated 

through mass loading would be biodegraded to ethene and 

not discharge with groundwater at downgradient boundaries, 

as predicted by the simulations. The goals of the transport 

simulations were to reproduce the observed configuration of 

the CVOCs plume and to estimate the mass-loading rates of 

CVOCs and the total mass of CVOCs in groundwater. The 

total-variation-diminishing (TVD) method with a Courant 

number of 1.0 was used to solve the advection-dispersion 

equation in the dual-solute simulation to minimize numerical 

dispersion while using small time steps (0.5 days) in the 50-yr 

simulation.

The mass of CVOCs beneath OU-1 increased rapidly 

at the beginning of the dual-solute simulation with Model A, 

then increased gradually and approached a dynamic equilibrium with a contaminant mass of 256 kg. Almost all the contaminant mass (134 g/d) discharged to small creeks draining 

the southern edge of OU-1, while the remainder discharged 

to the marsh pond and creek (1.7 g/d) and the tide flats (49 

mg/d). The CVOC mass-loading rate was 150 g/d, within the 

range of biodegradation rates of 130 to 170 g/d estimated by 

Dinicola (2006). The CVOC mass was 310 kg in Model B, but 

the mass-loading rate is only 73 g/d because there is a lower 

volume of groundwater flow in the deeper Qv sediments than 

in Model A, due to the higher vertical anisotropy in Unit 5 

(Qvsemi). As a result, less CVOC mass (48 g/d) discharges 

to the small creeks draining the southern edge of OU-1. More 

CVOC mass discharges to the marsh pond and creek (7.8 g/d) 

and the tide flats (160 mg/d), however, because contaminant 

plumes are diverted westward in Model B to Dogfish Bay by 

the higher K values in Unit 2 (Qestuary).

Two distinct contaminant plumes were simulated in 

Model A: a northern plume with 30 percent of the contaminant 

mass that migrates toward the marsh pond and creek and a 

southern plume with 70 percent of the contaminant mass that 

migrates toward small creeks south of OU-1. Both CVOC 

plumes extend from the water table to the bottom of the Qv 

sediments where they migrate laterally and upward to discharge in the marsh and creeks. The northern plume extends 

beneath the tide flats where it mixes with higher-density 

saline water, which presents a hydraulic barrier that limits the 

upward discharge of CVOCs to the tide flats. The contaminant 

plumes emanating from two persistent sources intermingle 

in Model B and become indistinguishable downgradient, but 

contaminant data are lacking to determine whether either the 

two distinct plumes predicted by Model A, or the single plume 

predicted by Model B more accurately depicts actual CVOC 

migration at OU-1.

Discussion of Simulation Results

Results of groundwater-flow and contaminant-transport 

simulations with both models A and B indicate that geologic 

and hydraulic barriers limit the migration of CVOCs from 

the former landfill at OU-1 and contain the CVOC plumes 

on Naval Base Kitsap. The Clover Park Aquitard is a lowpermeability geologic barrier that limits vertical groundwater 

flow between shallow, contaminated Qv sediments at OU-1 

and deeper hydrogeologic units. In addition, hydraulic barriers 

limit downward migration of CVOCs through the Clover Park 

Aquitard and lateral migration of CVOCs into upland areas 

south of the Keyport Peninsula. Finally, saline water beneath 

the tide flats and in Dogfish and Liberty Bays limits westward 

CVOC migration toward residential wells on Virginia Point.

Uncertainty in the rates of recharge and groundwater 

flow produces uncertainty in the estimated mass-loading rate 

of CVOCs at OU-1. Higher rates of recharge and groundwater 

flow require larger mass-loading rates to maintain a stable 

configuration of the contaminant plume. The lower massloading rate in Model B (73 g/d), however, indicates that 

higher vertical anisotropy can restrict groundwater flow rates 

through deeper Qv sediments beneath OU-1. The estimated 

CVOC mass-loading rate, therefore, ranges from 73 to 240 

g/d. Although the contaminants 1,4-dioxane and vinyl chloride 

have been detected at low concentrations on the State Route 

308 causeway northwest of the tide flats, neither Model A 

nor B predicted groundwater flow or contaminant transport 

from OU-1 to this area. There is uncertainty as to whether the 

contaminants measured in this area originated in the landfill or 

elsewhere.

Previous studies established that reductive dechlorination of CVOCs is occurring at OU-1 with biodegradation rates 

ranging from 130 to 170 g/d. A first-order, removal rateconstant of 10-3 d-1 was computed from CVOC contaminant 

masses computed by the transport simulations of models A and 

B and the biodegradation rate (150 g/d) estimated by Dinicola 

(2006). An alternate 50-yr, dual-solute contaminant-transport 

simulation with Model A and a first-order reaction predicted 

that a stable 77-kg CVOC plume formed within 10 years of 

migration, using the same mass-loading rates as in the original transport simulation. One half the CVOC mass was lost 

through biodegradation (75 g/d), and nearly all the remainder 

(71 g/d) was discharged by the southern plume to a small 

creek south of OU-1. The CVOC mass discharged from the 

northern plume to the marsh and tide flats was 0.95 mg/d. The 

areas of the northern and southern plumes were 67 percent and 

32 percent less, respectively, than the areas computed by the 

original transport simulation.
56 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

Both models A and B provide an adequate fit to hydraulic 

and contaminant data, but additional data are required to determine which model more accurately predicts actual contaminant migration at OU-1. Both models predict similar contaminant masses and biodegradation rate constants and also predict 

that contaminant migration is limited by the same geologic 

and hydraulic barriers. Model predictions differ in terms of 

the hydraulic parameters (recharge rates, and horizontal and 

vertical hydraulic conductivities), and in the CVOC massloading rates and contaminant plume configurations. Model B 

better represents the effects of tidal fluctuations on unconfined 

storage, uses recharge rates closer to the SWB estimated rates, 

and predicts water-table altitudes below land surface near 

OU-1. Model A better predicts the groundwater-level response 

to tidal fluctuations, however, and was chosen to demonstrate 

the effects representation of tidal fluctuations and potential 

remedial activities, although the effects would be the same if 

Model B were chosen instead. A soon-to-be-released version 

of MODFLOW 6 that supports variable-density flow and 

contaminant transport, and the distribution of discontinuous 

hydrogeologic units with the NWT solver would better represent unconfined and confined storage and greatly simplify the 

simulation of groundwater and contaminant transport at OU-1.

Two alternate 50-yr, contaminant-transport simulations 

assessed the sensitivity of predicted CVOC migration from 

OU-1 to the depiction of tidal fluctuations in the original transport simulation with Model A. In steady-state simulation C, 

which used the TVD solver, mean tides were lowered by 1.2 ft 

in Dogfish and Liberty Bays, 0.87 ft in the tide flats, and 0.1 ft 

in Marsh Creek (the averages of water levels measured during 

the 31-day period in July and August 2018). In transient simulation D, the FD solver was applied to compute 600 repetitions 

of the 31.05-day transient periods that represented diurnal 

tidal fluctuations used for calibration of the groundwaterflow model.

Both alternative representations of tidal fluctuations 

in simulations C and D had a small effect on the sizes and 

masses of the simulated CVOC plumes, when compared with 

the results of the original, steady-state TVD and FD simulations with Model A. The lower high and low tide stages in 

simulation C predicted larger CVOC plumes (272 kg) than 

in the original TVD simulation. The northern CVOC plume 

extended about 100 ft farther under the tide flats and CVOC 

discharge to the tide flats increased (from 49 to 890 mg/d) 

and CVOC discharge to Marsh Creek and pond decreased 

(from 1,700 to 770 mg/d). Representation of transient tidal 

fluctuations in simulation D predicted smaller CVOC plumes 

(297 kg), and the northern CVOC plume extended about 75 

ft farther under the tide flats. The CVOC discharge to the tide 

flats increased (from 67 mg/d to 75 mg/d), and the width of the 

transition zone between fresh water and seawater north of the 

marsh also increased (from 100 ft to 125 ft). The effect of tidal 

fluctuations on the predicted CVOC plumes is small because 

the period of the fluctuations (6 hr) is rapid compared with 

the groundwater velocity (20 ft/yr) beneath OU-1, and tidal 

fluctuations only slightly alter groundwater-flow paths and 

lengthen the northern CVOC plume under the tide flats.

Three 50-yr, steady-state contaminant-transport simulations with Model A predicted the effects of potential sourcecontrol scenarios on the contaminant plumes by reducing 

CVOC mass-loading rates by factors of 50, 90, and 99 percent, 

respectively. Under each scenario, the mass of CVOCs 

decreased rapidly and then slowly approached a new dynamic 

equilibrium. The reductions of CVOC mass in groundwater 

were 46, 86, and 95 percent, respectively, but sizes of the 

contaminant plumes were largely unchanged.

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. [Also available at https://doi.org/10.3133/ofr00342.]

Battelle Memorial Institute, 2018, 2017 Site recharacterization, Phase II, Operable Unit 1, Naval Base Kitsap, Keyport, 

Washington: Prepared by Battelle Memorial Institute, 

Columbus, Ohio for Naval Facilities Engineering Command 

Northwest, Silverdale, Wash., variously paged.

Collias, E.E., McGary, Noel, and Barnes, C.A., 1974, Atlas 

of physical and chemical properties of Puget Sound and 

its approaches: Seattle, Wash. and London, Washington 

Sea Grant Publication distributed by the University of 

Washington Press, 235 p.

Department of the Navy, 1993, Remedial Investigation Report, 

NUWC Keyport: Prepared by URS Consultants and Science 

Applications International Corporation for EFA NW under 

CLEAN Contract No. N62474-89-D-9295, CTO 10.

Department of the Navy, 2015, Fourth Five-Year Review, 

Naval Base Kitsap, Keyport, Washington—Prepared 

by Naval Facilities Engineering Command Northwest: 

Washington, Silverdale, 444 p.

Dinicola, R.S., Cox, S.E., Landmeyer, J.E., and Bradley, 

P.M., 2002, Natural attenuation of chlorinated volatile 

organic compounds in ground water at Operable Unit 

1, Naval Undersea Warfare Center, Division Keyport, 

Washington: U.S. Geological Survey Water-Resources 

Investigations Report 2002–4119, 115 p. [Also available at 

https://pubs.er.usgs.gov/publication/wri024119.]
References Cited 57

Dinicola, R.S., 2006, Continued biodegradation of chloroethene compounds in ground water at Operable Unit 1, Naval 

Undersea Warfare Center, Division Keyport, Washington: 

U.S. Geological Survey Scientific Investigations Report 

2006–5056, 42 p. [Also available at https://doi.org/10.3133/

sir20065056]

Frans, L.M., and Olsen, T.D., 2016, Numerical simulation 

of the groundwater-flow system of the Kitsap Peninsula, 

west-central Washington (ver. 1.1, October 2016): U.S. 

Geological Survey Scientific Investigations Report 

2016–5052, 63 p. [Also available at https://pubs.er.usgs.gov/

publication/sir20165052.]

Garling, M.E., Molenaar, Dee, and others, 1965, Water 

resources and geology of the Kitsap Peninsula and certain adjacent islands: Olympia, Washington, Washington 

Department of Conservation, Washington Division of Water 

Resources Water-Supply Bulletin 18, 309 p.

Gelhar, L.W., Welty, C., and Rehfeldt, K.W., 1992, A critical review of data on field-scale dispersion in aquifers: 

Water Resources Research, v. 28, no. 7, p. 1955–1974, 

https://doi.org/10.1029/92WR00607.

Gendaszek, A.S., and Welch, W.B., 2018, Water budget of 

the upper Chehalis River Basin, southwestern Washington: 

U.S. Geological Survey Scientific Investigations Report 

2018–5084, 17 p. [Also available at https://doi.org/10.3133/

sir20185084.]

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. [Also 

available at https://doi.org/10.3133/ofr200092.]

Harbaugh, A.W., 2005, MODFLOW-2005, The U.S. 

Geological Survey modular ground-water model—the 

Ground-Water Flow Process: U.S. Geological Survey 

Techniques and Methods 6–A16, variously paged. [Also 

available at https://pubs.usgs.gov/tm/2005/tm6A16/.]

Hargreaves, G.H., and Samani, Z.A., and the George H. 

Hargreaves, and the Zohrab A. Samani, 1985, Reference 

crop evapotranspiration from temperature: Applied 

Engineering in Agriculture, v. 1, no. 2, p. 96–99, 

https://doi.org/10.13031/2013.26773.

Headman, A.O, 2020, Soil water balance (SWB) model of 

Keyport, Washington: U.S. Geological Survey data release, 

https://doi.org/10.5066/P95WQ7TM.

Jones, M.A., 1999, Geologic framework of the Puget Sound 

aquifer system, Washington and British Columbia: U.S. 

Geological Survey Professional Paper 1424–C, 31 p. [Also 

available at https://doi.org/10.3133/pp1424C.]

La Licata, I., Langevin, C.D., and Dausman, A.M., 2007, 

Effect of tidal fluctuations on contaminant transfer to 

the ocean: A New Focus on Groundwater–Seawater 

Interactions, Proceedings of Symposium HS1001 at 

IUGG2007, Perugia, July 2007, IAHS Publ. 312, 2007.

Langevin, C.D., Thorne, D.T., Jr., Dausman, A.M., Sukop, 

M.C., and Guo, Weixing, 2007, SEAWAT Version 4—A 

computer program for simulation of multi-species solute 

and heat transport: U.S. Geological Survey Techniques and 

Methods, book 6, chap. A22, 39 p. [Also available at https://

pubs.er.usgs.gov/publication/tm6A22.]

Keller, C.K., van der Kamp, G., and Cherry, J.A., 1989, A 

multiscale study of the permeability of a thick clayey till: 

Water Resources Research, v. 25, no. 11, p. 2299–2317, 

https://doi.org/10.1029/WR025i011p02299.

Miller, D.A., and White, R.A., 1998, A conterminous United 

States multilayer soil characteristics dataset for regional 

climate and hydrology modeling: Earth Interactions, v. 2, 

no. 2, p. 1–26, https://doi.org/10.1175/1087-3562(1998)002

<0001:ACUSMS>2.3.CO;2.

Niswonger, R.G., Panday, Sorab, and Ibaraki, Motomu, 

2011, MODFLOW-NWT, A Newton formulation for MODFLOW-2005: U.S. Geological Survey 

Techniques and Methods 6–A37, 44 p. [Also available at 

https://pubs.usgs.gov/tm/tm6a37/.]

National Oceanic and Atmospheric Administration, 2018, 

Tides and currents: National Oceanic and Atmospheric 

Administration, web, accessed August 2018 at 

https://tidesandcurrents.noaa.gov/noaatidepredictions.html?

id=9445719&units=standard&bdate=20180712&edate=

20180808&timezone=LST/LDT&clock=12hour&datum=

MLLW&interval=hilo&action=data.

Opatz, C.C., and Dinicola, R.S., 2019, Analysis of groundwater response to tidal fluctuations, Operable Unit 1, Naval 

Base Kitsap, Keyport, Washington: U.S. Geological Survey 

Open-File Report 2019–1098, 36 p. [Also available at 

https://doi.org/10.3133/ofr20191098.]

Poeter, E.P., and Hill, M.C., 1998, Documentation of UCODE, 

a computer code for universal inverse modeling: U.S. 

Geological Survey Water-Resources Investigations Report 

98–4080, 116 p.
58 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

U.S. Geological Survey, 2014, Hydrography—National 

hydrography dataset, watershed boundary dataset: U.S. 

Geological Survey Web site, accessed August 27, 2018, at 

https://nhd.usgs.gov/index.html

URS Consultants Inc, 1997, Final summary data assessment report for Operable Unit 1 Naval Undersea Warfare 

Center Division Keyport, Washington—Prepared by URS 

Consultants, Seattle, Washington, for Engineering Field 

Activity 3 vols: Poulsbo, WA, Northwest, Naval Facilities 

Engineering Command.

URS Consultants Inc, 1998, Final record of decision for 

Operable Unit 1 Naval Undersea Warfare Center Division 

Keyport, Washington—Prepared by URS Consultants, 

Seattle, Wash., for Engineering Field Activity: Poulsbo, 

Wash., Northwest, Naval Facilities Engineering 

Command, 111 p.

Vaccaro, J.J., Hansen, A.J., Jr., and Jones, M.A., 1998, 

Hydrogeologic framework of the Puget Sound aquifer 

system, Washington and British Columbia: U.S. Geological 

Survey Professional Paper 1424–D, 77 p. [Also available at 

https://doi.org/10.3133/pp1424D

Welch, W.B., Frans, L.M., and Olsen, T.D., 2014, 

Hydrogeologic framework, groundwater movement, and water budget of the Kitsap Peninsula, westcentral Washington: U.S. Geological Survey Scientific 

Investigations Report 2014–5106, 44 p. [Also available at 

https://doi.org/10.3133/sir20145106.]

Westenbroek, S.M., Kelson, V.A., Dripps, W.R., Hunt, 

R.J., and Bradbury, K.R., 2010, SWB—A modified 

Thornthwaite-Mather soil-water balance code for estimating 

groundwater recharge: U.S. Geological Survey Techniques 

and Methods, book 6, chap. A31, 60 p. [Also available at 

https://pubs.usgs.gov/tm/tm6-a31/.]

Wilson, J.T., Kampbell, D.H., and Weaver, J.W., 1996: 

Environmental chemistry and the kinetics of biotransformation of chlorinated organic compounds in ground 

water, in U.S. Environmental Protection Agency, 

Symposium in natural attenuation of chlorinated organics 

in ground waters, September 11–13, 1996: Dallas, Texas, 

EPA/540/R-96/509.

Yager, R.M., and Fountain, J.C., 2001, Effect of natural 

gas exsolution on specific storage in a confined aquifer undergoing water level decline: Ground Water, v. 

39, no. 4, p. 517–525. https://doi.org/10.1111/j.1745-

6584.2001.tb02340.x.

Yager, R.M., and Kappel, W.M., 1998, Infiltration and hydraulic connections from the Niagara River to a fractureddolomite aquifer in Niagara Falls, New York: Journal 

of Hydrology (Amsterdam), v. 206, no. 1-2, p. 84–97, 

https://doi.org/10.1016/S0022-1694(98)00088-2.

Yager, R.M., 2020, MODFLOW-2005, MODFLOW-NWT, 

and S EAWAT V.4 models used to simulate variable-density 

groundwater flow and contaminant transport at Naval Base 

Kitsap, Keyport, Washington: U.S. Geological Survey data 

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

Zheng, Chunmiao, and Bennett, G.D., 2002, Applied contaminant transport modeling: John Wiley and Sons, Inc., New 

York, 621 p.

Zheng, C., and Wang, P.P., 1999, MT3DMS—A modular 

three-dimensional multispecies transport model for simulation of advection, dispersion and chemical reactions of 

contaminants in ground-water systems–Documentation 

and user’s guide: U.S. Army Corps of Engineers Contract 

Report SERDP-99-1.
Appendix 1. Soil-Water Balance (SWB) Model 59

Appendix 1. Soil-Water Balance (SWB) Model

The Soil-Water-Balance (SWB) computer program 

(Westenbroek and others, 2010) provides estimates for 

groundwater recharge by calculating water-balance components at daily time steps for each model cell using a modified Thornthwaite Mather soil-water-balance approach. The 

model documented in this report utilizes SWB ver1.2. The 

model input, output, and executable files are available from 

Headman (2020).

Soil-Water Balance (SWB) Model 

Spatially Distributed Datasets

Five spatially distributed (gridded) datasets were used 

in the SWB recharge simulations. All datasets were in ASCII 

format (with the exception of the climate data, which was in 

NetCDF format), and cover the 5.9 mi2 study area with a uniform cell size of 25 ft. The input datasets include daily climate 

data (precipitation, maximum and minimum air temperature), 

land cover, flow direction, hydrologic soil group and available 

soil-water capacity (table 1.1).

Daily Climate Data

Daily climate data required for the SWB recharge 

calculation include precipitation (inches) and maximum and 

minimum temperature (Fahrenheit). Daily evapotranspiration is computed using the method of Hargreaves and Samani 

(1985). Daily climate data were obtained from the Daymet 

dataset available from Oak Ridge National Laboratory 

(http://daymet.ornl.gov). The Daymet climate data is computed through interpolation from spatially referenced observation stations and available in 2-degree areas at a resolution 

of 0.62 mi for all of North America. Daily Daymet climate 

data were obtained for the study area for the period 1980 

through 2015.

Land Cover

The 2011 National Land Cover Database (NLCD; 

Homer and others, 2015) was used to characterize land 

cover for the study area (table 1.2). The NLCD is produced 

by the Multi-Resolution Land Characteristics consortium 

(http://www.mrlc.gov) to describe land-cover characteristics across the United States using a 16-class classification 

schema at a resolution of 100 ft. The NLCD was obtained 

from https://www.mrlc.gov/nlcd11_data.php using an Albers 

Conical Equal Area projection (NAD83). The NLCD data 

was resampled to the 25-ft cell size of the SWB model using a 

Table 1.1. Spatial data assigned to parameters in the soil-water 

balance model for Keyport study area.

[NLCD, National Land Cover Database; STATSGO, State Soil Geographic 

Database; ft, feet; --, not applicable (vector data)]

Type Data source Resolution Soil-water balance 

parameters

Land Cover NLCD 20111 100 ft Land Cover

Soils

gSSURGO2 or 

STATSGO23 –

Runoff Curve 

number

Root depth

Maximum Recharge

Available Water 

Capacity

Flow 

Direction

Puget Sound 

LiDAR 6 ft D8 Flow Direction

Climate Daymet 0.62 mi Temperature

Precipitation

Evapotranspiration

1Homer and others (2015).

2National Resources Conservation Service, 2018

3Miller and White (1998).

Table 1.2. Distribution of land cover in the Keyport study area in 

2011.

[mi2, square mile]

Land 

use 

code

Description Area

(mi2)

Percentage 

of study 

area

11 Open water 1.85 31.3

21 Developed, open space 1.06 17.8

22 Developed, low intensity 0.64 10.9

23 Developed, medium intensity 0.16 2.7

24 Developed, high intensity 0.10 1.7

31 Barren land 0.15 2.5

41 Deciduous forest 0.49 8.3

42 Evergreen forest 0.45 7.6

43 Mixed forest 0.56 9.5

52 Shrub and scrub 0.07 1.1

71 Grassland/herbaceous 0.04 0.7

81 Pasture/hay 0.04 0.7

90 Woody wetlands 0.14 2.4

95 Emergent herbaceous wetlands 0.15 2.6
60 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

majority sampling technique in ArcGIS. A comparison of the 

two gridded datasets indicates that the spatial pattern for the 

principal land-cover classes in the study area is similar. Most 

of the land area is either developed or covered by forest.

Flow Direction

A flow-direction grid to route surface runoff 

between SWB model cells used the D8 flow-routing 

convention of Lehner and others (2008). The flowdirection grid was generated from Lidar elevation data 

at a 6 ft resolution obtained from Puget Sound Lidar 

(http://pugetsoundlidar.ess.washington.edu/). Individual Lidar 

files were combined and resampled to the 25-ft cell size of the 

SWB model using a majority sampling technique in ArcGIS. 

The resulting elevation grid was processed to smooth over 

small sinks that were artifacts of data processing and the D8 

flow-direction grid was then generated with the ArcGIS “flow 

direction” tool.

Hydrologic Soil Group and Available Soil-Water 

Capacity

Two gridded datasets of soil properties, hydrologic 

soil group (HSG) and available water capacity (AWC), 

are required for the SWB model. HSG and AWC data are 

compiled by the U.S. Department of Agriculture Natural 

Resource Conservation Service (NRCS) for the entire United 

States. NRCS classifies HSGs in four groups (A, B, C and 

D) on a continuum from high infiltration capacity with low 

runoff potential to low infiltration capacity with high runoff 

potential. AWC is the amount of water a particular soil column 

can hold at various depths; the AWC for the top meter of soil 

was used for the SWB model described in this report. Data for 

HSG and AWC were mostly obtained from the NRCS Gridded 

Soil Survey Geographic (gSSURGO) database (National 

Resources Conservation Service, 2018). Some areas for which 

gSSURGO data were not available used the NRCS State Soil 

Geographic (STATSGO2) database (Miller and White, 1998). 

HSG and AWC grids were resampled to 25-ft cells using the 

majority technique for HSG data and the mean sampling technique for AWC data.

Tabular Input Datasets

Tabular information, which includes runoff curve numbers, vegetation routing depths, maximum infiltration rates, 

and maximum daily recharge rates, are also required for SWB 

model for each combination of HSG and Land Cover type. 

The table was populated with values used by Gendaszek 

and Welch (2018), Tillman (2015) and Westenbroek (2010) 

(table 1.3).
Appendix 1. Soil-Water Balance (SWB) Model 61

Table 1.3. Tabular data for soil-water balance (SWB) model of Keyport study area.

Land 

use 

code

Description

Impervious 

area

(percent)

Curve numbers Maximum recharge rate (inches/day) Interception Root zone depth (ft)

A1 B1 C1 D1 A1 B1 C1 D1 Growing 

season

Nongrowing 

season

A1 B1 C1 D1

11 Open Water 0 100 100 100 100 2 0.6 0.12 0.24 0 0 0 0 0 0

21 Developed, Open Space 0 49 69 79 84 2 0.6 0.12 0.24 0.0835 0 8.53 8.53 8.53 8.5322 Developed, Low Intensity 30 77 86 91 94 2 0.6 0.12 0.24 0.0835 0 8.53 8.53 8.53 8.5323 Developed, Medium Intensity 60 89 92 94 95 2 0.6 0.12 0.24 0.0835 0 8.53 8.53 8.53 8.5324 Developed, High Intensity 90 98 98 98 98 2 0.6 0.12 0.24 0.0835 0 8.53 8.53 8.53 8.5331 Barren Land 0 77 86 91 94 2 0.6 0.12 0.24 0 0 1 1 1 141 Deciduous Forest 0 32 48 57 63 2 0.6 0.12 0.24 0.0835 0 9.5 9.5 9.5 9.542 Evergreen Forest 0 39 58 73 80 2 0.6 0.12 0.24 0.0835 0 12.8 12.8 12.8 12.843 Mixed Forest 0 46 60 68 74 2 0.6 0.12 0.24 0.0835 0 11.15 11.15 11.15 11.1552 Shrub/Scrub 0 49 68 79 84 2 0.6 0.12 0.24 0.0835 0 3.5 3.5 3.5 3.571 Grassland/Herbaceous arid 0 64 71 81 89 2 0.6 0.12 0.24 0.0835 0 8.53 8.53 8.53 8.5381 Pasture/Hay fair 0 49 69 79 84 2 0.6 0.12 0.24 0.0835 0 8.53 8.53 8.53 8.5382 Cultivated Crops SR+CR 

poor

0 71 80 87 90 2 0.6 0.12 0.24 0.0835 0 2 2 2 2

90 Woody Wetlands 0 88 89 90 91 2 0.6 0.12 0.24 0.0835 0 4.5 4.5 4.5 4.595 Emergent Herbaceous Wetlands 0 89 90 91 92 2 0.6 0.12 0.24 0 0 4.5 4.5 4.5 4.5

1Hydrologic soil group (HSG).
62 Variable-Density Groundwater Flow and Contaminant Transport, Operable Unit 1, Naval Base Kitsap, Keyport, Washington

References Cited

Carslaw, H.S., and Jaeger, J.C., 1959, Conduction of heat in 

solids (2nd ed.): New York, Oxford University Press, 510 p.

Gendaszek, A.S., and Welch, W.B., 2018, Water budget of 

the upper Chehalis River Basin, southwestern Washington: 

U.S. Geological Survey Scientific Investigations Report 

2018–5084, 17 p. [Also available at https://doi.org/10.3133/

sir20185084.]

Hargreaves, G.H., and Samani, Z.A., and the George H. 

Hargreaves, and the Zohrab A. Samani, 1985, Reference 

crop evapotranspiration from temperature: Applied 

Engineering in Agriculture, v. 1, no. 2, p. 96–99, 

https://doi.org/10.13031/2013.26773.

Headman, A.O, 2020, Soil water balance (SWB) model of 

Keyport, Washington: U.S. Geological Survey data release, 

https://doi.org/10.5066/P95WQ7TM.

Homer, C.G., Dewitz, J.A., Yang, L., Jin, S., Danielson, P., 

Xian, G., Coulston, J., Herold, N.D., Wickham, J.D., and 

Megown, K., 2015, Completion of the 2011 National Land 

Cover Database for the conterminous United States—

Representing a decade of land cover change information: 

Photogrammetric Engineering and Remote Sensing, v. 81, 

no. 5, p. 345–354.

Lehner, B., Verdin, K., and Jarvis, A., 2008, New global 

hydrography derived from spaceborne elevation data: Eos 

(Washington, D.C.), v. 89, no. 10, p. 93–94, https://doi.org/

10.1029/2008EO100001.

Miller, D.A., and White, R.A., 1998, A conterminous United 

States multilayer soil characteristics dataset for regional 

climate and hydrology modeling: Earth Interactions, v. 2, 

no. 2, p. 1–26, https://doi.org/10.1175/1087-3562(1998)002

<0001:ACUSMS>2.3.CO;2.

National Oceanic and Atmospheric Administration, 2020, 

Tides and currents: National Oceanic and Atmospheric 

Administration, web.

National Resources Conservation Service, 2018, 

Web Soil Survey: United States Department of 

Agriculture, web, accessed March 19, 2018, at 

https://websoilsurvey.sc.egov.usda.gov/.

Tillman, F.D., 2015, Documentation of input datasets for 

the soil-water balance groundwater recharge model of 

the Upper Colorado River Basin: U.S. Geological Survey 

Open-File Report 2015–1160, 17 p. [Also available at 

https://doi.org/10.3133/ofr20151160.]

Westenbroek, S.M., Kelson, V.A., Dripps, W.R., Hunt, 

R.J., and Bradbury, K.R., 2010, SWB—A modified 

Thornthwaite-Mather soil-water balance code for estimating groundwater recharge: U.S. Geological Survey 

Techniques and Methods, book 6, chap. A31, 60 p., 

https://pubs.usgs.gov/tm/tm6-a31/
Publishing support provided by the U.S. Geological SurveyScience 

Publishing Network, Tacoma Publishing Service Center

For more information concerning the research in this report, contact the 

Director, Washington Water Science Center

U.S. Geological Survey

934 Broadway, Suite 300

Tacoma, Washington 98402

https://wa.water.usgs.gov
Yager and others—Variable-Density Groundwater Flow and Contaminant Transport, Naval Base Kitsap, Keyport, Washington—SIR 2020–5066

ISSN 2328-0328 (online)

https://doi.org/10.3133/sir20205066