Use of Inverse-Modeling Methods to 

Improve Ground-Water-Model Calibration 

and Evaluate Model-Prediction Uncertainty, 

Camp Edwards, Cape Cod, Massachusetts

By Donald A. Walter and Denis R. LeBlanc

Prepared in cooperation with the

Impact Area Groundwater Study Program of the

National Guard Bureau and U.S. Army Environmental Command

Scientific Investigations Report 2007–5257

U.S. Department of the Interior

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

DIRK KEMPTHORNE, Secretary

U.S. Geological Survey

Mark D. Myers, Director

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

For product and ordering information:

World Wide Web: http://www.usgs.gov/pubprod

Telephone: 1-888-ASK-USGS

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

natural hazards, and the environment:

World Wide Web: http://www.usgs.gov

Telephone: 1-888-ASK-USGS

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

U.S. Government.

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

reproduce any copyrighted materials contained within this report.

Suggested citation:

Walter, D.A., and LeBlanc, D.R., 2008, Use of inverse-modeling methods to improve ground-water-model calibration 

and evaluate model-prediction uncertainty, Camp Edwards, Cape Cod, Massachusetts: U.S. Geological Survey 

Scientific Investigations Report 2007–5257, 57 p., available online at http://pubs.usgs.gov/sir/2007/5257.
iii

Contents

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

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

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

Site Description and History ...............................................................................................................4

Hydrogeology of the Camp Edwards Area .......................................................................................4

Geologic Setting...........................................................................................................................4

Hydrologic Setting .......................................................................................................................7

Methods of Analysis.............................................................................................................................8

Numerical Ground-Water-Flow Model.......................................................................................................9

Regional Model Design........................................................................................................................9

Modifications to the Regional Model ..............................................................................................11

Representation of Parameters ................................................................................................13

Incorporation of Observed Heads and Flows .......................................................................16

Use of Inverse Methods to Improve Model Calibration ........................................................................18

Evaluation of Model Sensitivity ........................................................................................................18

Formulation of the Least-Squares Objective Function .................................................................23

Heads and Flows.................................................................................................................................23

Improvement in Calibration......................................................................................................25

Optimal Parameter Estimates ..................................................................................................25

Precision of Parameter Estimates ..........................................................................................29

Modification to the Model Framework...................................................................................29

Advective-Transport Observations ..................................................................................................31

Limitations of Analysis .......................................................................................................................37

Evaluation of Model-Prediction Uncertainty...........................................................................................40

Forward Particle Tracking .................................................................................................................42

Reverse Particle Tracking .................................................................................................................44

Summary and Conclusions.........................................................................................................................47

Acknowledgments.......................................................................................................................................48

References Cited..........................................................................................................................................48
iv

Figures

1–2. Maps showing—

1. Location of Massachusetts Military Reservation and Camp Edwards on 

western Cape Cod, Massachusetts..................................................................................2

2. Simulated water table and flow directions and Camp Edwards training 

installations, Royal Dutch Explosive (RDX) contaminant plumes, and production 

wells, northwestern Cape Cod, Massachusetts.............................................................5

3. Map and diagram showing (A) surficial geology and regional water table, western 

Cape Cod, Massachusetts, and (B) hydrologic section through western Cape Cod........6

4. Map and cross section showing (A) areal extent of regional model grid, modeled 

hydrologic boundaries, and the simulated water table; and vertical grid discretization 

for the linearized model (model 2, table 1) along (B) row 97, and (C) column 107, 

western Cape Cod, Massachusetts.........................................................................................10

5–6. Maps showing—

5. Comparison of simulated water tables for the trial-and-error, linearized, and 

parameterized regional models, western Cape Cod, Massachusetts......................12

6. Distribution of (A) hydraulic conductivity and boundary leakance parameter 

zones in the parameterized regional model (model 3, table 1), and (B) hydraulic 

conductivity in the trial-and-error regional model (model 1, table 1), western 

Cape Cod, Massachusetts ...............................................................................................14

7. Cross sections showing vertical distribution of hydraulic conductivities in the 

trial-and-error regional model (model 1, table 1) along (A) row 97, and (B) column 107; 

and vertical distribution of hydraulic conductivity and boundary leakance parameter 

zones in the parameterized model (model 3, table 1) along (C) row 97, and 

(D) column 107.............................................................................................................................15

8. Graph showing absolute mean head residuals for the trial-and-error and 

parameterized regional models of western Cape Cod, Massachusetts ...........................16

9. Map showing locations of head and streamflow observations, Sagamore flow lens, 

western Cape Cod, Massachusetts.........................................................................................17

10. Graphs showing composite-scaled sensitivities of heads to (A) all parameters 

(horizontal and vertical hydraulic conductivity and boundary leakances), and 

(B) horizontal hydraulic conductivity parameters for vertical group 1..............................20

11. Maps showing (A–F) the 1-percent scaled sensitivities of simulated heads to 

selected parameters (model 3, table 1), western Cape Cod, Massachusetts..................21

12. Map and graphs showing (A) locations of selected observations on western 

Cape Cod, Massachusetts, and (B–J) the 1-percent scaled sensitivities for 

the observation locations..........................................................................................................22

13–15. Graphs showing—

13. Composite-scaled sensitivities of heads (for modified parameter distribution; 

model 4, table 1) to (A) all parameters, and (B) horizontal hydraulic conductivity 

and boundary leakance parameters ..............................................................................24

14. Absolute mean head residuals calculated for the parameterized model 

(model 4, table 1) by using parameter values derived from the trial-and-error 

calibration and optimal parameters from the inverse calibration.............................25

15. Magnitudes of changes calculated in residuals at individual observation 

locations for (A) heads and (B) streamflows for the parameterized model 

(model 4, table 1) by using parameter values derived from the trial-and-error 

calibration and optimal parameter values from the inverse calibration..................26
v

16. Map showing spatial distribution and magnitudes of changes calculated in absolute 

head residuals for the parameterized model (model 4, table 1) by using parameter 

values derived from the trial-and-error calibration and optimal parameter values from 

the inverse calibration, western Cape Cod, Massachusetts ..............................................27

17–18. Graphs showing—

17. (A) Initial and final estimated parameters and precisions of the estimates 

(as standard deviations) for parameters included in the parameter-estimation 

regression, and (B) estimated and specified parameter values for different 

parameter regressions for model 4 and variants 1 and 2 (table 1)............................28

18. The 1-percent sensitivity of advective-transport predictions in the central part 

of the Impact Area to estimated parameters and precision of the parameter 

estimates expressed as standard deviations and coefficients of variation............30

19. Map showing spatial extent of the RDX plume downgradient from Demolition 

Area 1, the advective-transport observation location, and particle tracks and 

advective-transport positions predicted by the parameterized model (model 4, 

variant 2, table 1) by using parameter values derived from the trial-and-error 

calibration and optimal parameter values from the inverse calibration, western 

Cape Cod, Massachusetts ........................................................................................................32

20–21. Graphs showing—

20. Sensitivity of the simulated three-dimensional particle position of the leading 

edge of a contaminant plume at Demolition Area 1 to hydraulic conductivity 

and leakance parameters for (A–C) original (model 4, variant 2, table 1) and 

(D–F) modified (model 5, table 1) aquifer zonation.......................................................34

21. Estimated and specified parameter values for regressions with and without 

inclusion of the advective-transport observation for (A) original aquifer zonation 

(model 4, variant 2, table 1) and (B) modified aquifer zonation (model 5, table 1) 

for various assumed traveltimes.....................................................................................35

22. Map and cross section showing (A) spatial extent of the RDX plume downgradient 

from Demolition Area 1 and particle tracks and advective-transport positions 

predicted by the parameterized models (model 4, variant 2, and model 5; table 1) by 

using parameter values derived from the trial-and-error calibration and optimal 

parameter values from inverse calibration to heads, flows, and the advectivetransport observation; and (B) vertical extent of the RDX plume and particle tracks 

from the inverse models for different assumed traveltimes to the advective-transport 

observation location, western Cape Cod, Massachusetts..................................................36

23–25. Graphs showing—

23. Absolute mean head residuals for versions of the inversely calibrated regional 

model with and without inclusion of the advective-transport observation and 

with different assumed traveltimes to the observation location...............................37

24. Comparisons of (A) weighted observations and weighted simulated equivalents, 

and (B) weighted simulated equivalents and weighted residuals for hydraulic 

head observations .............................................................................................................39

25. Normal-order statistic and weighted residuals for hydraulic head observations .40

26. Schematic diagrams showing uncertainties associated with (A) forward and 

(B) reverse particle tracks represented as confidence intervals on particle positions.41

27. Map and cross section showing (A) forward particle track and advective-transport 

positions near the Central Impact Area and 68-, 95-, and 99-percent confidence 

intervals associated with the predictions; and (B) vertical forward particle track and 

vertical 95-percent confidence intervals associated with the prediction, western 

Cape Cod, Massachusetts ........................................................................................................43
vi

28. Graph showing (A–C) the 1-percent sensitivity of simulated particle positions at 

three locations along the forward particle track to selected model parameters............44

29–30. Maps showing—

29. Simulated recharge locations for eight subsurface contaminant detections in 

and around the Central Impact Area and approximate areas of 68-, 95-, and 

99-percent confidence associated with the predicted recharge locations, 

western Cape Cod, Massachusetts................................................................................45

30. Simulated recharge area to a production well near the Impact Area and 

approximate areas of 68-, 95-, and 99-percent confidence associated with the 

predicted recharge area, western Cape Cod, Massachusetts..................................46

Tables

1. Numerical models developed to represent the aquifer of western Cape Cod, 

Massachusetts, and agreement between observed and simulated heads......................11

2. Model parameters used to represent glacial sediments on western Cape Cod, 

Massachusetts............................................................................................................................13
vii

Conversion Factors, Datums, and Abbreviations

Inch/Pound to SI

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 mile (mi2

) 2.590 square kilometer (km2)

Volume

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

)

cubic foot (ft3

) 28.32 cubic decimeter (dm3)

cubic foot (ft3

) 0.02832 cubic meter (m3)

Flow rate

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)

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

/s)

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

Hydraulic conductivity

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

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 follows:

 °F=(1.8×°C)+32

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

 °C=(°F–32)/1.8

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

(NGVD 29).

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

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

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

times foot of aquifer thickness [(ft3

/d)/ft2]ft. In this report, the mathematically reduced form, 

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

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

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

OTHER ABBREVIATIONS USED IN REPORT

AMEC AMEC Earth & Environmental, Inc.

CIA Central Impact Area

CV coefficient of variation

Demo 1 Demolition Area 1

DNT dinitrotoluene

HK horizontal hydraulic conductivity

HMX Her Majesty’s Explosive

IAGWSP Impact Area Groundwater Study Program

K hydraulic conductivity

MassDEP Massachusetts Department of Environmental Protection

MMR Massachusetts Military Reservation

NGB National Guard Bureau

RDX Royal Dutch Explosive

USEPA U.S. Environmental Protection Agency

USGS U.S. Geological Survey

UXO unexploded ordnance

VK vertical hydraulic conductivity
Use of Inverse-Modeling Methods to Improve GroundWater-Model Calibration and Evaluate Model-Prediction 

Uncertainty, Camp Edwards, Cape Cod, Massachusetts

By Donald A. Walter and Denis R. LeBlanc

Abstract

Historical weapons testing and disposal activities at 

Camp Edwards, which is located on the Massachusetts 

Military Reservation, western Cape Cod, have resulted in the 

release of contaminants into an underlying sand and gravel 

aquifer that is the sole source of potable water to surrounding communities. Ground-water models have been used at the 

site to simulate advective transport in the aquifer in support 

of field investigations. Reasonable models developed by 

different groups and calibrated by trial and error often yield 

different predictions of advective transport, and the predictions 

lack quantitative measures of uncertainty. A recently (2004) 

developed regional model of western Cape Cod, modified to 

include the sensitivity and parameter-estimation capabilities 

of MODFLOW-2000, was used in this report to evaluate the 

utility of inverse (statistical) methods to (1) improve model 

calibration and (2) assess model-prediction uncertainty.

Simulated heads and flows were most sensitive to 

recharge and to the horizontal hydraulic conductivity of the 

Buzzards Bay and Sandwich Moraines and the Buzzards Bay 

and northern parts of the Mashpee outwash plains. Conversely, 

simulated heads and flows were much less sensitive to vertical 

hydraulic conductivity. Parameter estimation (inverse calibration) improved the match to observed heads and flows; 

the absolute mean residual for heads improved by 0.32 feet 

and the absolute mean residual for streamflows improved by 

about 0.2 cubic feet per second. Advective-transport predictions in Camp Edwards generally were most sensitive to the 

parameters with the highest precision (lowest coefficients of 

variation), indicating that the numerical model is adequate 

for evaluating prediction uncertainties in and around Camp 

Edwards. The incorporation of an advective-transport observation, representing the leading edge of a contaminant plume 

that had been difficult to match by using trial-and-error calibration, improved the match between an observed and simulated plume path; however, a modified representation of local 

geology was needed to simultaneously maintain a reasonable 

calibration to heads and flows and to the plume path.

Advective-transport uncertainties were expressed as 

about 68-, 95-, and 99-percent confidence intervals on threedimensional simulated particle positions. The confidence 

intervals can be graphically represented as ellipses around 

individual particle positions in the X-Y (geographic) plane 

and in the X-Z or Y-Z (vertical) planes. The merging of 

individual ellipses allows uncertainties on forward particle 

tracks to be displayed in map or cross-sectional view as a cone 

of uncertainty around a simulated particle path; uncertainties 

on reverse particle-track endpoints—representing simulated 

recharge locations—can be geographically displayed as areas 

at the water table around the discrete particle endpoints. This 

information gives decisionmakers insight into the level of confidence they can have in particle-tracking results and can assist 

them in the efficient use of available field resources.

Introduction

Live-fire training and munitions disposal at Camp 

Edwards on the Massachusetts Military Reservation (MMR), 

Cape Cod, Massachusetts (fig. 1), have resulted in the release 

of contaminants into the environment. The underlying sand 

and gravel aquifer is the sole source of potable water to 

the residents of western Cape Cod, and there is concern 

that migration of contaminants from Camp Edwards could 

adversely affect current and future water-supply resources. In 

1997, the U.S. Environmental Protection Agency (USEPA) 

issued an administrative order that suspended live-fire training 

at Camp Edwards until ground-water contamination at the site 

could be evaluated. As a result, the National Guard Bureau 

(NGB) initiated the Impact Area Groundwater Study Program 

(IAGWSP) to assess ground-water and soil contamination 

at the site. The contaminant of most concern is Royal Dutch 

Explosive (RDX), a potentially carcinogenic compound that 

is transported conservatively in ground water under the oxic 

conditions that are generally observed in the sand and gravel 

aquifer. Other contaminants of concern at the site include 

Her Majesty’s Explosive (HMX), perchlorate, picric acid, and 

dinitrotoluene (DNT).
2 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 1. Location of Massachusetts Military Reservation and Camp Edwards on western Cape Cod, Massachusetts.

Cape Cod Canal

Massachusetts

Military

Reservation

Massachusetts

Military

Reservation

Sandwich

Mashpee

Falmouth

Barnstable

Bourne

 Cape Cod Bay

Vineyard Sound

Buzzards Bay

Snake 

Pond

Mashpee 

 Pond

Santuit River

Pocasset 

 River

Bourne 

 River

Backus 

 River

Squeteague

Brook

Quaker Run

Shawme

 Pond

Little 

River

Mashpee River

Quashnet River

Childs River

Coonamessett River

Marston Mills River

Massachusetts

Military

Reservation

Otis Air

Force Base

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator projection,

Zone 19, 1:24,000, 1991

Wastewater-treatment

facility

Hatchville

precipitation

gage

Impact 

 Area 

Impact 

 Area 

Camp

Edwards 

Massachusetts

Military

Reservation

Otis Air

Force Base

Massachusetts

Cape

Cod

Sagamore flow lens

Extent of map

41o42'30"

41o

32'30"

70o

37'30" 70o27'30"

0

0

5 MILES

5 KILOMETERS
Introduction 3

The U.S. Geological Survey (USGS) has assisted the 

IAGWSP in the investigation since 1997. One role of the 

USGS has been to improve understanding of the regional 

ground-water-flow system and to provide assistance in evaluating ground-water contamination at the site in a hydrologic 

context. USGS activities include the development and application of steady-state and transient numerical models to simulate the advective transport of contaminants in the aquifer; 

advective transport is the movement of a solute at the average 

ground-water velocity. Numerical models have been used at 

the site by USGS and other researchers to provide support for 

IAGWSP field investigations. Model-predicted particle tracking has been used to determine subsurface sampling locations 

and the placement of observation wells downgradient from 

known contaminant sources, as identified by soil sampling, 

and to identify possible source areas for contaminants detected 

in the subsurface during reconnaissance drilling. Numerical 

models also have been used to delineate areas at the water 

table that contribute water to production wells. 

Contaminants emanating from sources at the former Otis 

Air Force Base (since 1973, Otis Air National Guard Base) 

(fig. 1), to the south of Camp Edwards, have formed large, 

distinct, and contiguous plumes of contaminated ground water. 

Source areas in this part of the MMR, which include training areas, industrial sites, and fuel spills, generally are well 

defined and have well-known source histories, and subsurface contamination downgradient from these sources has 

been delineated by detailed ground-water-sampling efforts. 

Conversely, source areas at Camp Edwards, which include 

unexploded ordnance (UXO) and scattered areas of munitions 

disposal, are not as well understood; consequently, plumes of 

contaminated ground water from these sources are less distinct 

and contiguous, particularly in and near the Impact Area, 

an area in Camp Edwards used for live-fire training (fig. 1). 

In these cases, field efforts have relied, to a large extent, on 

numerical models to determine downgradient drilling locations 

and to link sporadic subsurface detections of contaminants to 

possible source areas.

The reliance on numerical models makes it necessary for 

models to be well calibrated to observations of heads, flows, 

and advective transport. Ground-water-flow models generally 

are calibrated by trial and error, which means that initial 

stresses and intrinsic aquifer properties are manually defined, 

and the models are calibrated to hydrologic data by a trialand-error method in which inputs into the model are varied 

manually to achieve a reasonable match to observed data. 

Although this approach can produce reasonably calibrated 

models, the final calibrated input parameters likely do not 

represent optimal values that give the best possible fit to 

the observations. Also, different reasonably calibrated trialand-error models can yield different model predictions. In 

addition to model calibration, it also is necessary to better 

understand the uncertainty of model predictions, particularly 

if the uncertainty relates to prediction of advective transport 

from particle tracking. Currently, no quantitative measure of 

particle-tracking uncertainty is associated with the trial-anderror models. Model-prediction uncertainty is particularly 

important when using particle tracking to guide field activities, 

such as determining locations for observation wells. A particle 

track represents an estimate of advective transport through 

the aquifer that is based on a given set of model-input values. 

Because uncertainty is inherent in the predictions, observation 

wells generally are placed close to the predicted particle track 

and at some distance from the predicted particle track to 

account for the unknown uncertainty. Without a quantitative 

measure of this uncertainty, it is difficult for decisionmakers 

to agree on an appropriate placement of wells that both 

adequately accounts for particle-tracking uncertainty and 

allows for the reasonable use of available resources, and 

thereby sets a technically justifiable limit to the number of 

required wells. There are several components of uncertainty 

associated with both model design and calibration. The work 

presented here addresses the model-calibration component of 

uncertainty; it should be noted, however, that uncertainties 

associated with the conceptual model of the system are 

inherent in any model predictions. 

In 2003, the USGS, in cooperation with the IAGWSP of 

the U.S. Army Environmental Command (USAEC), began an 

evaluation of model calibration and model-prediction uncertainty at Camp Edwards. The overall objective of the investigation was to evaluate the utility of inverse-modeling methods 

to improve model calibration and evaluate model-prediction 

uncertainty. The specific objectives of the investigation were 

to (1) quantitatively evaluate sensitivities of simulated heads, 

flows, and advective transport to input parameters; (2) use 

inverse modeling to statistically estimate optimal model-input 

parameters; and (3) use the resulting measures of statistical 

precision of the estimated parameters to quantify uncertainties 

associated with particle-tracking predictions.

Purpose and Scope

This report documents the use of inverse-modeling 

methods to quantify model sensitivities, improve model 

calibration, and evaluate particle-tracking uncertainties in and 

around Camp Edwards. The report presents an example that 

(1) illustrates the use of these methods to improve diagnostic 

understanding of the sensitivities of simulated hydraulic 

conditions to model-input parameters and (2) shows how 

model-calculated sensitivities can assist in trial-and-error 

calibration and how the sensitivities are used when estimating 

optimal parameter values. The report also presents an 

example of the use of these sensitivities along with parameter 

estimation to improve model calibration. Examples are 

presented of the use of parameter estimation to calibrate to 

heads and flows as well as the incorporation of an observation 

of advective transport into inverse-model calibration. In 

addition, the report illustrates the use of statistical measures 

of precision from estimated input parameters to quantify the 

uncertainty associated with estimates of advective transport 

determined from particle tracking for three types of advective-
4 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

transport predictions: a forward particle-tracking prediction 

from a potential contaminant source, reverse particle tracking 

from the locations of subsurface contaminant detections to 

the water table, and a model-delineated contributing area to a 

production well.

Site Description and History

Military activity at the MMR, a multiuse facility that 

encompasses about 22,000 acres on western Cape Cod, 

began as early as 1911. Camp Edwards encompasses about 

14,000 acres in the north-central part of the MMR and 

consists of the Impact Area, which is about 2,000 acres in 

area, surrounded by several training ranges. The site was 

operated by the U.S. Army until about 1974 and is now used 

as a training facility by the Massachusetts National Guard. 

Since the mid-1930s, the Impact Area has been used for 

live-fire mortar and artillery training, and the surrounding 

training ranges have been used for small-arms training and 

troop maneuvers. Other historical activities at the site include 

ordnance training and the testing and disposal of ordnance 

by the military and military contractors. These activities have 

resulted in the release of contaminants into the underlying 

aquifer. The contaminants of most concern are explosive 

compounds and perchlorate. Explosive compounds (RDX) 

and perchlorate have been detected in a number of areas of 

Camp Edwards (U.S. Army Environmental Command, 2007a; 

2007b). Areas where plumes of subsurface contamination have 

been mapped and remedial actions are planned or underway 

include Demolition Area 1 (Demo 1), the J-Ranges, and the 

Central Impact Area (CIA) (AMEC Earth & Environmental, 

Inc. (AMEC), 2005; 2002; 2001) (fig. 2). 

The J-Ranges were used for mortar, rocket, and smallarms training activities from about 1935 to the mid-1980s. 

From the late 1960s to the mid-1990s, the J-Ranges also were 

used as a test facility by military contractors for testing of 

munitions (AMEC, 2002). Although little is known about the 

exact nature of activities that occurred in this area, a number 

of structures that likely were used to test and dispose of ordnance have been identified. Activities in this area have resulted 

in contamination of the underlying ground water, where 

observed RDX concentrations are as high as 30 μg/L (AMEC, 

2002). Contamination emanating from the J-Ranges has been 

observed as far downgradient as Snake Pond, where RDX has 

been detected beneath the northern part of the pond. Contamination from the J-Ranges has been observed as deep as 120 ft 

below the water table (AMEC, 2002). 

The CIA, located in the central part of the Impact Area, 

has been used as a live-fire mortar and artillery training area 

since about 1940. Inert and live ordnance are commonly 

observed in the CIA, particularly near targets. The release of 

explosive compounds into the environment occurs primarily 

through either low-order (partial) or high-order detonation of 

shells (AMEC, 2001). The pattern of contamination, which 

emanates from a large number of small, isolated sources, 

resembles the release of contamination from a nonpoint source 

and occurs as diffuse, widely distributed contaminants in the 

aquifer. Maximum observed RDX concentrations in the CIA 

are about 45 μg/L, and contaminants emanating from the 

CIA have been observed 3,000 ft downgradient from the CIA 

boundary (AMEC, 2001). Contamination emanating from the 

CIA has a maximum depth of about 100 ft below the water 

table. The exact locations of all potential source areas in the 

J-Ranges and CIA and detailed source histories are not well 

known; general areas of subsurface contamination have been 

delineated by drilling and sampling.

Demo 1 was used for demolition training and the disposal 

of unexploded ordnance by detonation from about 1986 to 

about 1997; prior to that time (from about the mid-1970s), the 

area was used as a training range (AMEC, 2005). These activities resulted in a well-defined plume of contaminated ground 

water that contains RDX and extends about 2,000 ft downgradient from the source. Unlike contamination in the CIA or 

J-Ranges, the location and source history at Demo 1 are fairly 

well known. The maximum measured RDX concentration 

in contaminated ground water at the site is about 250 μg/L 

beneath the source area (AMEC, 2005). 

Hydrogeology of the Camp Edwards Area

The glacial sediments beneath southeastern Massachusetts generally are sandy and highly permeable and can exceed 

500 ft in thickness in the south-central part of Cape Cod. 

Average annual rainfall in the region is about 45 in. (LeBlanc 

and others, 1986), and the unconsolidated sandy sediments 

underlying the region are the sole source of potable water for 

local communities. The geologic history and hydrology of 

Cape Cod have been documented in numerous publications, 

including LeBlanc and others (1986), Masterson and others 

(1997a), Oldale (1992), and Uchupi and others (1996). The 

nomenclature for the unconsolidated deposits on Cape Cod is 

documented in Oldale and Barlow (1986).

Geologic Setting

The unconsolidated sediments underlying Cape Cod 

consist of gravel, sand, silt, and clay that were deposited 

during the Pleistocene Epoch between 15,000 and 

16,000 years ago (Oldale and Barlow, 1986). Sediment was 

deposited at or near the edges of a retreating continental ice 

sheet through direct deposition from the ice and proglacial 

deposition from meltwater. The surficial geology of western 

and central Cape Cod is characterized by broad, gently sloping 

glacial-outwash plains and hummocky terrain associated 

with glacial moraines (fig. 3). A number of different types 

of glacial deposits are found on Cape Cod. Moraines were 

deposited at or near the edge of the ice lobes and are either 

ablation moraines that were deposited in place by melting 

ice, such as the Buzzards Bay Moraine, or tectonic moraines, 

such as the Sandwich Moraine, that consist of reworked 
Introduction 5

Figure 2. Simulated water table and flow directions and Camp Edwards training installations, Royal Dutch Explosive 

(RDX) contaminant plumes, and production wells, northwestern Cape Cod, Massachusetts.

Bourne

Snake

Pond

60

EXPLANATION

Massachusetts Military Reservation

boundary

Simulated line of equal water-table altitude—

Contour interval is 5 feet. Datum is NGVD 29

General direction of ground-water flow

Production well

Impact Area

Ranges

Training areas

Demolition Area 1 RDX

contaminant plume

Other RDX contaminant plumes

Production well #3

J-Ranges

RDX Contamination

40

30

50

60

10

20

Camp Edwards 

Impact 

 Area 

Demolition

 Area 1

 J-Ranges

Top of watertable mound 

Central Impact

 Area

 

Sandwich

Mashpee

RDX contamination

Production well #3

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator projection,

Zone 19, 1:24,000, 1991

Cape Cod Canal

41o45'30"

41o

40'30"

70o

35'30" 70o30'30"

0

0

2 MILES

2 KILOMETERS

Plumes of ground-water contamination mapped

by U.S. Army Environmental Command (2007a,b)
6 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 3. (A) Surficial geology and regional water table, western Cape Cod, Massachusetts, and (B) hydrologic section 

through western Cape Cod.

Top of the

water-table mound 

VERTICAL SCALE GREATLY EXAGGERATED

30

40

50

60

10

20

Buzzards Bay

Location of generalized

hydrologic section 

River

 Camp 

Edwards 

EXPLANATION

A

Monomoy

flow lens Sagamore flow lens

Bass

A

B

Cape Cod Canal

Nantucket Sound

Vineyard Sound

Cape Cod Bay

Mashpee outwash

Barnstable outwash

Harwich outwash

Buzzards Bay outwash

Sandwich Moraine

Buzzards Bay Moraine 

Nantucket ice-contact deposits

Cape Cod Bay lacustrine deposits 

Falmouth ice-contact deposits 

Massachusetts Military

Reservation boundary

Simulated line of equal water-table 

altitude—Contour interval is 5 feet. 

Datum is NGVD 29

General direction of groundwater flow

Top of the

Water-Table Mound 

Saline 

ground 

water

Bedrock

Freshwater-saltwater 

interface

 NGVD 29

-100

-200

-300

 200

 100

North South

Vineyard

 Sound

-400

-500

-600

Fresh ground water

Cape Cod

 Bay

Areal recharge

A A'

A'

Top of the

water-table mound 

2020

ALTITUDE, IN FEET 

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator projection,

Zone 19, 1:24,000, 1991

41o42'30"

41o

32'30"

70o

37'30" 70o

17'30"

0

0

10 MILES

10 KILOMETERS

0

0

10 MILES

10 KILOMETERS
Introduction 7

outwash sediments pushed into place by local readvances of 

the ice lobes (Uchupi and others, 1996; Oldale, 1992; Oldale 

and O’Hara, 1984). Kames are ice-contact deposits that were 

deposited in high-energy meltwater environments in holes 

in the ice sheets. Ice-contact deposits also include sediments 

that were deposited by meltwater in high-energy fluvial 

environments near the ice margin. Outwash sediments, the 

most common unconsolidated sediments underlying Cape 

Cod, were deposited by meltwater streams in depositional 

environments associated with proglacial lake deltas; these 

depositional environments are analogous to those found in 

present-day fluvial deltas (Uchupi and others, 1996; Oldale, 

1992). Deltaic sediments can be divided into three general 

facies: topset, foreset, and bottomset deposits (Masterson 

and others, 1997a). These sediments generally become finegrained with depth; in general, coarse-grained fluvial sand 

and gravel (topset sediments) overlie fine to medium sand 

deposited in nearshore lacustrine environments (foreset beds), 

which in turn are underlain by silty sand and clay deposited in 

offshore lacustrine environments (bottomset beds). Collapse 

structures within the outwash deposits, often associated with 

topographic depressions, formed when buried ice blocks 

melted, causing the collapse of overlying sediments; coarsegrained sediments often extend deeper in the aquifer in these 

areas than elsewhere. The glacial sediments described above 

are underlain by basal till in most places; basal till consists 

primarily of low-permeability clay that was produced by 

mechanical erosion of bedrock by movement of the overlying 

ice sheet. The unconsolidated glacial sediments are underlain 

by relatively impermeable crystalline bedrock.

The sequence of glacial deposits on western Cape Cod 

ranges in thickness from about 70 ft near the Cape Cod Canal 

to more than 500 ft along Nantucket Sound. The lithology 

of the glacial sediments, which include moraines, kames and 

other ice-contact deposits, and outwash, differs according to 

the depositional environment in which the sediments were 

deposited. Outwash sediments generally are characterized 

by two grain-size trends: (1) a fining-downward trend in 

which grain size decreases and silt and clay content increases 

with depth, and (2) a trend in which sediments become finer 

grained and coarse-grained sand and gravel deposits become 

thinner with increasing distance from the sediment sources. 

Most of western Cape Cod is underlain by Mashpee Pitted 

Plain deposits, which generally are coarse-grained, wellsorted glacial outwash deposits (Oldale and Barlow, 1986). 

The sediment source for these deltaic sediments was close to 

the apex of the Sandwich and Buzzards Bay Moraines near 

the Cape Cod Canal (fig. 3A); as a result, sediments in the 

Mashpee Pitted Plain deposits generally become finer grained 

to the south. Moraine sediments, which consist of gravel, 

sand, silt, and clay, were deposited in low-energy depositional 

environments and generally are poorly sorted and have a 

highly variable lithology. The Buzzards Bay and Sandwich 

Moraines are near the shores of Buzzards Bay and Cape Cod 

Bay, respectively (fig. 3A). Ice-contact and kame deposits are 

found along the shore of Nantucket Sound (fig. 3A). Basal 

till, which is fine-grained and consists primarily of clay, was 

deposited at the base of the ice sheets. 

Grain size and degree of sorting determine the watertransmitting properties of aquifer sediments. The trends in 

hydraulic conductivity (K) of outwash sediments generally are 

parallel to the trends in grain size; the hydraulic conductivity 

of sediments generally decreases with depth and with 

increasing distance from sediment sources, or generally 

southward (Masterson and others, 1997a). Aquifer tests during 

previous investigations have identified general relations 

between sediment grain size and hydraulic conductivity 

(Masterson and others, 1997a; Masterson and Barlow, 1997). 

Medium- to coarse-grained sand and gravel deposits have 

hydraulic conductivities that range from 200 to 350 ft/d. Fine 

to medium sands have hydraulic conductivities typically 

ranging from 70 to 200 ft/d. The hydraulic conductivity of 

very fine sand and silt typically ranges from 30 to 70 ft/d; 

silt and clay deposits have hydraulic conductivities of 10 to 

30 ft/d. Ice-contact and kame deposits generally consist of 

medium- to coarse-grained sand and gravel and have hydraulic 

conductivities similar to those of coarse-grained outwash 

deposits. Moraine deposits have a variable lithology—ranging 

from gravel and sand to silt and clay—and generally have 

lower average hydraulic conductivities than outwash deposits. 

Hydraulic conductivity decreases with depth in most areas, 

including the moraines.

Camp Edwards is underlain by moraine and outwash 

deposits (fig. 3A). The Impact Area is on the outwash plain 

near and to the southeast of the sediment source. As a result, 

the area is underlain by coarse-grained sediments characteristic of a high-energy depositional environment. The sediments 

generally consist of medium- to coarse-grained sand and 

gravel with local deposits of silt and fine sand, particularly 

deeper in the aquifer. Moraine deposits north and west of 

the Impact Area consist of gravel, sand, silt, and clay; have a 

more variable lithology; and generally are more fine-grained 

than the adjacent outwash deposits. Saturated thickness in 

the Impact Area ranges from about 100 to 250 ft (Walter and 

Whealan, 2005). 

Hydrologic Setting

Camp Edwards is located on the Sagamore flow lens 

(fig. 1), which is the largest and westernmost lens of six 

separate ground-water-flow lenses that underlie Cape Cod 

(LeBlanc and others, 1986). The unconsolidated glacial 

sediments underlying the region compose an unconfined 

aquifer system that is surrounded by saltwater: Cape Cod 

Bay to the northeast, Buzzards Bay to the west, and Nantucket 

Sound to the south (fig. 3A). The Sagamore flow lens is 

hydraulically separated at its northwestern extent from 

mainland Massachusetts by the Cape Cod Canal and at its 

eastern extent from the adjacent Monomoy flow lens 

(fig. 3A) by the Bass River and associated wetlands. The 

aquifer system is bounded below by impermeable bedrock and 

at the top by the water table, across which recharge enters 
8 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

(fig. 3B). Recharge from precipitation is the sole source of 

water to the aquifer system. About 45 in. of precipitation 

falls annually on Cape Cod; about half of the precipitation 

recharges the aquifer (LeBlanc and others, 1986). The 

remainder is lost to evapotranspiration; surface runoff is 

negligible owing to the sandy soils and low topographic relief 

of the area. 

Ground water flows outward from regional ground-water 

divides towards natural discharge locations at streams, coastal 

estuaries, and the ocean (fig. 3A). Most flow is in shallow 

sediments and discharges to streams and estuaries; ground 

water recharging the aquifer near the central ground-water 

divides flows deeper in the aquifer and discharges to the ocean 

(fig. 3B). About half of the ground-water discharge is into 

saltwater bodies. About 40 percent of ground-water discharge 

is into freshwater streams and wetlands, and a small amount 

(about 6 percent) is removed from the system for water supply 

(Walter and Whealan, 2005). Water-table contours and local 

ground-water-flow patterns are strongly affected by numerous 

kettle-hole ponds. These ponds are flow-through ponds 

characterized by ground-water-flow paths converging in areas 

upgradient of the ponds, where ground water discharges into 

the ponds, and diverging in downgradient areas, where pond 

water recharges the aquifer. Streams generally are areas of 

ground-water discharge (gaining streams) and receive water 

from the aquifer over most of their length. Some stream 

reaches may lose water to the aquifer (losing streams), 

particularly in areas downgradient from pond outflows. 

Regional ground-water-flow patterns on western Cape 

Cod, in the vicinity of Camp Edwards, are dominated by a 

water-table mound, the top of which is southeast of the Impact 

Area (fig. 2) (Walter and Whealan, 2005). Water-level altitudes in this area exceed 65 ft on average. Ground water flows 

radially outward from the mound and, as a result, groundwater-flow directions in Camp Edwards differ depending on 

location relative to the position of the top of the mound 

(fig. 2). This radial flow field has important implications for 

the advective transport of contaminants from Camp Edwards. 

In the Impact Area, ground water flows to the northwest 

towards the Cape Cod Canal, whereas ground water from 

the southern J-Ranges flows southward toward Snake Pond. 

Ground water near Demo 1 flows westward towards Buzzards 

Bay (fig. 2). Ground-water-flow patterns have a stronger component of vertical flow near the top of the water-table mound 

than in areas away from the mound. Ground-water flow 

is nearly vertically downward at the top of the water-table 

mound and nearly horizontal in downgradient areas of the 

aquifer; flow near discharge boundaries has a strong component of upward vertical flow (fig. 3B). Contaminants from the 

J-Ranges, which enter the ground-water-flow system closer to 

the top of the water-table mound, would be expected to move 

deeper in the system relative to horizontal transport distance 

than contaminants from sources farther from the mound, such 

as Demo 1.

Because of high recharge rates (about 27 in/yr) and the 

generally high permeability of the aquifer sediments, advective transport is the dominant component of contaminant transport in the aquifer (LeBlanc, 1984). Ground-water velocities 

of more than 1.5 ft/d have been observed along the southern 

boundary of the MMR downgradient from the wastewatertreatment facility (fig. 1; LeBlanc and others, 1991). At Camp 

Edwards, RDX has been observed in ground water at distances 

as great as 2 mi from likely sources in the Impact Area and at 

depths as great as 100 ft below the water table (AMEC, 2006).

The rate of advective transport in the aquifer is a function 

of the location of source areas relative to regional groundwater divides; water recharged near divides, where groundwater flow is more nearly vertical, moves more slowly than 

water recharged away from divides, where horizontal flow 

predominates (fig. 3B) (Walter and Masterson, 2003; Walter 

and others, 2004). Traveltimes, defined as the total time it 

takes water to move from recharge locations at the water table 

to natural discharge locations range from essentially zero 

adjacent to discharge boundaries to hundreds of years near 

ground-water divides (Walter and others, 2004).

Methods of Analysis

A numerical ground-water-flow model that represents the 

hydrologic system of the Sagamore flow lens on western Cape 

Cod was developed in 2004 by the USGS in cooperation with 

the Massachusetts Department of Environmental Protection 

(MassDEP). This model, which forms the basis of the work 

described in this report, uses the USGS three-dimensional 

ground-water-flow model MODFLOW-2000 (Harbaugh 

and others, 2000) to simulate regional ground-water flow. 

A detailed documentation of the model, including model 

design, hydrologic boundaries, hydrologic stresses, aquifer 

characteristics, and model calibration, is presented in Walter 

and Whealan (2005). The advective transport of contaminants 

in the aquifer was simulated by using the USGS particletracking software algorithm MODPATH (Pollock, 1994). 

MODFLOW-2000 and its suite of supporting programs 

allowed for the incorporation of parameter estimation into 

the modeling analysis. Specifically, the Observation (OBS), 

Sensitivity (SEN), and Parameter Estimation (PES) Processes 

(Hill and others, 2000) were used to incorporate observed 

hydraulic data into the model and to use these observations to 

evaluate model sensitivities and estimate optimal parameter 

values. The 68-, 95-, and 99-percent confidence intervals on 

predicted particle tracks based on optimal parameter values 

were estimated by using the program YCINT-2000 (Hill and 

others, 2000).
Numerical Ground-Water-Flow Model 9

Numerical Ground-Water-Flow Model

The regional model of the Sagamore flow lens that 

forms the basis of the analyses documented in this report was 

developed as part of a recent (2001–04) investigation into the 

regional hydrology of central and western Cape Cod (Walter 

and Whealan, 2005). The purpose of this previous modeling 

effort, which was done in cooperation with MassDEP, was 

to evaluate (1) the areas contributing recharge to production 

wells and natural receptors (ponds, streams, and estuaries) in 

the region, and (2) the effects of time-varying recharge and 

pumping on surface water in the region. The regional model 

of the Sagamore flow lens (as well as a regional model of 

the adjacent Monomoy flow lens) is documented in detail in 

Walter and Whealan (2005, app. 1). 

Regional Model Design

The regional model of the Sagamore flow lens extends 

from the Cape Cod Canal in Bourne to the Bass River in 

Yarmouth (fig. 3A). The finite-difference ground-watermodeling program MODFLOW-2000 (Harbaugh and others, 

2000; McDonald and Harbaugh, 1988) was used to simulate 

the ground-water system of the Sagamore lens. The model has 

246 rows and 365 columns with a total active modeled area of 

252 mi2

; the horizontal discretization is a uniform 400 by 

400 ft (fig. 4A). The model has 20 layers with thicknesses 

ranging from 10 ft in the top 170 ft of saturated thickness to 

more than 200 ft in the deepest layer (figs. 4B and C). 

Estuaries, open coastal waters, and streams are 

represented as head-dependent hydraulic boundaries in the 

model (fig. 4A). Open coastal waters, saltwater estuaries, 

and some freshwater wetlands were simulated by using the 

General Head Boundary (GHB) and Drain (DRN) Packages 

in MODFLOW-2000 (Harbaugh and others, 2000; McDonald 

and Harbaugh, 1988). The saltwater-boundary altitudes 

surrounding the aquifer range from 0.5 to 1.9 ft and are based 

on either tidal-gage information from the National Oceanic 

and Atmospheric Administration, where available, or on 

freshwater-equivalent heads estimated from bathymetry. 

Freshwater streams and some freshwater wetlands were 

simulated by using the Streamflow Routing Package (STR) 

Package (Prudic, 1989), which allows for an accounting of 

streamflow and pumping-induced streamflow depletion. Ponds 

were simulated as active parts of the aquifer but were assigned 

a very large hydraulic conductivity (100,000 ft/d) so that there 

was effectively no resistance to flow through the ponds.

The only input of water into the model is from 

recharge. A natural recharge rate onto aquifer sediments of 

27.1 in/yr was specified; this value was obtained from longterm precipitation records from Hatchville, Mass. (fig. 1), 

was adjusted during model calibration, and is consistent with 

previous investigations on Cape Cod and southeastern New 

England (Barlow and Dickerman, 2001; DeSimone and others, 

2002; Masterson and others, 1997b; Walter and Masterson, 

2003). Natural recharge onto surface-water bodies varied: 

16 in/yr onto ponds, representing net recharge after pan 

evaporation, and no recharge onto wetlands, which likely are 

areas of net ground-water discharge. Recharge was further 

adjusted to account for septic-system return flow in residential 

areas. Parcel-scale water-use data were not available, so 

the volume of generated wastewater was determined from 

the total volume of ground water withdrawn by each town; 

a consumptive loss of 15 percent was assumed, and the 

remaining volume was evenly distributed across model cells 

representing non-sewered residential areas. Recharge also was 

adjusted in areas of existing wastewater-disposal facilities 

in Barnstable, Falmouth, and at the MMR (fig. 1). In 2004, 

about 96 production wells were operating on western Cape 

Cod (Walter and Whealan, 2005) and withdrew a total of 

about 17.4 Mgal/d of water from the aquifer. The four towns 

bordering the MMR (Bourne, Falmouth, Mashpee, and 

Sandwich) withdrew a total of about 7.7 Mgal/d of water. 

Withdrawal of water at production wells was simulated using 

the WEL Package (McDonald and Harbaugh, 1988). 

Aquifer properties were estimated from lithologic logs 

and from previously developed depositional models of western Cape Cod (Masterson and others, 1997b; Byron Stone, 

U.S. Geological Survey, written commun., 2002). Hydraulic 

conductivities varied spatially and with depth and ranged from 

350 ft/d for coarse sand and gravel to 10 ft/d for silt and clay 

(Walter and Whealan, 2005). The steady-state regional model 

was calibrated by using measured water levels and streamflows that represent long-term average hydraulic conditions, 

and delineated contaminant plumes as indicators of advective 

transport. Initial estimates of recharge and intrinsic aquifer 

properties were adjusted during the trial-and-error calibration 

process to achieve an acceptable match between observed 

and simulated hydrologic conditions at the calibration points 

(Walter and Whealan, 2005). 

The regional model is based on the assumption of steadystate hydraulic stresses (recharge and pumping). Although 

these stresses vary over time and versions of the regional 

model that incorporate these time-varying stresses were 

developed, steady-state models are adequate for the simulation 

of advective transport in most areas, including the CIA and 

around Demo 1, because the time scale of the variability of 

these stresses is small compared to the time scale of transport in the aquifer. Walter and Masterson (2003) compared 

advective-transport paths computed from a transient flow field 

in which recharge varied yearly (between 14 and 41 in/yr) to 

flow paths computed with a steady-state, average recharge 

rate. The results showed that computed flow paths in most 

areas of Camp Edwards, including Demo 1 and the CIA, were 

nearly identical. This similarity indicates that the effect of 

transient recharge on advective-transport paths was negligible 

in most areas and that the assumption of steady-state conditions in simulating advective transport generally is valid.
10 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 4. (A) Areal extent of regional model grid, modeled hydrologic boundaries, and the simulated water table; and vertical grid 

discretization for the linearized model (model 2, table 1) along (B) row 97, and (C) column 107, western Cape Cod, Massachusetts. 

In the original regional model (model 1, table 1), layers 1–7 have a uniform thickness of 10 feet.

Bedrock

Bedrock

Aquifer

Aquifer

VERTICAL SCALE GREATLY EXAGGERATED

-100

 NGVD 29

-140

-240

 60

Layer 19

Layer 20

Layer 18

Layers 8–17

Layers 1–7

-100

-140

-240

 60

Layer 19

Layer 20

Layer 18

Layers 8–17

Layers 1–7

60

50

40

10

20

30

Row 97

Column 107

Row 97

Column 107

A

B

C

ALTITUDE, IN FEET

ALTITUDE, IN FEET NGVD 29

41o42'30"

41o

32'30"

70o37'30" 70o

17'30"

0

0

10 MILES

10 KILOMETERS

202

Coastal boundary (GHB)

Estuary boundary (DRN)

Saltwater boundary

Pond

Stream (STR)

Freshwater boundary

EXPLANATION

Simulated line of equal

water-table altitude—

Contour interval is 5 feet. 

Datum is NGVD 29

Production well
Numerical Ground-Water-Flow Model 11

The aquifer is unconfined, and the water table is 

simulated as a free surface. The bottoms of the top six layers 

of the model are uniform and are above sea level. As a result, 

model layers are free to go dry and rewet during the model 

simulation. In areas with simulated heads of less than 10 ft, 

the top six layers are dry; no model layers are dry in areas with 

simulated heads of more than 60 ft. 

Modifications to the Regional Model

The regional model described above was modified to 

facilitate the implementation of parameter estimation in the 

analysis of model calibration and uncertainty. Before performing inverse-model calibration, modifications were made that 

included changes in (1) model layering to minimize possible 

nonlinearities and (2) data structure to allow for the representation of model inputs as parameters. Descriptions of the trialand-error and inversely calibrated models developed as part of 

this investigation and the match to observed heads, as represented by the absolute mean head residual, are summarized in 

table 1. In this report, the original regional model is referred to 

as the trial-and-error model (model 1, table 1).

The drying and rewetting of model cells introduce nonlinearities into the model solution that may interfere with the 

parameter-estimation process. To minimize this problem, the 

top and bottom altitudes of the upper layers—greater than an 

altitude of 0 ft—were modified to ensure that no model cells 

went dry during the simulation. This was accomplished by 

setting the bottom altitude of layer 1 to 2 ft below the calibrated water-table altitude for a given row and column, and 

evenly spacing the bottoms of layers 2–6 between the bottoms 

of layer 1 and layer 7 (altitude of 0 ft) for that row and column 

(figs. 4B and C). This allowed the uppermost model layer 

(layer 1) to remain saturated during the simulation, resulting in 

a more linear and efficient solution. The simulated water table 

produced from the linearized model (model 2, table 1) was 

nearly identical to that of the model with the original layering 

scheme (fig. 5) and resulted in a similar match to observed 

heads (table 1); this result would be expected because model 

inputs (aquifer properties, boundary conditions, and hydraulic stresses) were unchanged. In addition to these changes, 

Table 1. Numerical models developed to represent the aquifer of western Cape Cod, Massachusetts, and agreement between 

observed and simulated heads.

[--, no value]

Model Model description

Absolute mean 

head residual, 

in feet

Parameter values from trial-and-error calibration, model 1

Trial-and-error model (1) 1

Original model, trial-and-error calibration 1.63

Linearized model (2) Model (1) modified to remove nonlinearities 1.62

Parameterized model (3) 2

Model (2) modified to represent hydraulic conductivity and boundary leakance 

parameters for simplified geology

1.68

Parameterized model (4) Model (3) simplified by lumping parameters from vertical groups 1 and 2 1.73

Parameter values from inverse calibration: original aquifer zonation

Parameterized model (4) Model (4) calibrated to head and flow observations 1.41

Variant 1: BUZZMS_12 extended to bedrock and FALMIC_3 specified at 

110 feet per day

--

Variant 2: BUZZMS_12 and FALMIC_3 specified at 350 and 110 feet per day, 

respectively

1.50

Parameterized model (4) Model (4), variant 2, calibrated to head, flow, and advective-transport observation; 

assumed traveltime of 25 years

2.13

Parameter values from inverse calibration: modified aquifer zonation

Parameterized model (5) Model (4), variant 2, with modified zonation of the Buzzards Bay Moraine calibrated 

to head, flow, and advective-transport observation; assumed traveltime of 25 years

1.58

Parameterized model (5) Model (5) calibrated to an advective-transport traveltime of 20 years 1.75

Parameterized model (5) Model (5) calibrated to an advective-transport traveltime of 15 years 1.81

1 Documented in Walter and Whealan (2005).

2 Shown in figures 7C and D.
12 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

41o42'30"

41o

32'30"

70o37'30" 70o

17'30"

30

40

50

60

10

20

202

202

202

202

Coastal boundary (GHB)

Estuary boundary (DRN)

Saltwater boundary

Pond

Stream (STR)

Freshwater boundary

EXPLANATION

Production well

Simulated line of equal water-table

altitude—Contour interval is 5 feet.

Datum is NGVD 29

 Trial-and-error model (1)

 Linearized model (2)

 Parameterized model (3)

 Parameterized model (4): Vertical

 groups 1 and 2 combined

Figure 5. Comparison of simulated water tables for the trial-and-error, linearized, and parameterized regional models, western 

Cape Cod, Massachusetts. Parameter values derived from trial-and-error calibration. Model numbers shown in table 1.
Numerical Ground-Water-Flow Model 13

the structure of the model was changed to allow for model 

inputs to be represented as parameters in the model instead of 

specified values. These parameters included intrinsic aquifer 

properties and boundary leakances. In addition, input files 

were created for the OBS, SEN, and PES Processes necessary 

for the analysis.

Representation of Parameters

In a trial-and-error model, all active model cells are 

assigned specified values of horizontal and vertical hydraulic 

conductivity (K); these values are determined during the trialand-error calibration process. In a parameterized model, these 

properties are represented as parameter values assigned to 

specified regions of the aquifer. Although several contiguous 

cells within the domain of a trial-and-error model can share 

a common value of hydraulic conductivity, input is required 

for each individual model cell in that region of the model. 

Conversely, parameters allow for a single value to be assigned 

to multiple cells within a region of the model. This allows for 

more efficient management of input data, and the parameters 

can vary during the estimation of optimal parameter values. 

In using parameters to define input properties, zones of the 

aquifer (representing the distribution of lithologic units) are 

defined, and each zone is assigned a common parameter value. 

The parameterized model (model 3, table 1) has a total 

of 34 separate hydraulic conductivity zones; as a result, there 

are a total of 34 horizontal hydraulic conductivity (HK) 

and 34 vertical hydraulic conductivity (VK) parameters in 

the model (table 2). Some simplification of model inputs 

was done during conversion to a parameterized model. The 

simplification was implemented to (1) remove unnecessary 

complexities from the trial-and-error model and (2) allow 

for a more efficient parameterization process. The degree of 

complexity in a parameterized model should be a function 

of the available observation data such that there is sufficient 

information in the calibration data to reliably estimate 

parameter values. The spatial distribution of hydraulic 

conductivity zones in the shallow parts of the aquifer 

(altitude greater than -10 ft) (fig. 6A) is similar to that in the 

original trial-and-error model (fig. 6B), because only minor 

simplifications were implemented. The parameterized model, 

to a similar extent as the trial-and-error model, represents the 

depositional model on which the aquifer zones are based. 

As an example, the shallowest part of the trial-and-error 

model (altitude greater than 60 ft) has a total of 16 zones 

(HK and VK), whereas the shallowest part of the parameterized model has 13 zones (HK and VK). Specifically, 

the Harwich and Barnstable outwash plain deposits were 

each simulated as a single hydraulic conductivity zone, and 

the Mashpee outwash deposits were subdivided into three 

hydraulic conductivity zones instead of the four zones in the 

trial-and-error model (fig. 6B).

The vertical distribution of hydraulic conductivity in the 

trial-and-error model is a function of model layering. Each of 

the 20 layers has a different distribution of hydraulic conductivity zonation reflecting the three-dimensionality of the 

depositional model (fig. 7). Therefore, as many as 20 separate 

zones can be present in any vertical profile within the model. 

This level of complexity was implemented to represent vertical trends in grain size in the aquifer and to represent local 

confining layers in some areas. The spatial complexity of 

hydraulic conductivity (along an X-Y plane) generally was 

preserved during conversion to the parameterized model 

(fig. 6), whereas the vertical distribution of hydraulic conductivity was simplified considerably (fig. 7). The 20 layers were 

Table 2. Model parameters used to represent glacial sediments 

on western Cape Cod, Massachusetts.

[Altitudes of vertical groups are as follows: Group 1 (70 to -10 feet); 

Group 2 (-10 to -100 feet); Group 3 (-100 to -240 feet); and Group 4 (deeper 

than -240 feet). Altitudes are relative to National Geodetic Vertical Datum 

of 1929. Parameter names shown are for horizontal hydraulic conductivity. 

Parameter names for vertical hydraulic conductivity include the letter “V” 

preceding the underscore.]

Parameter Description

Vertical extent 

of parameter 

zones

MASH1_1,2,3 Northern part of the Mashpee 

Pitted Plain

Groups 1, 2, 

and 3

MASH2_1,2,3 Central part of the Mashpee 

Pitted Plain

Groups 1, 2, 

and 3

MASH3_1 Southern part of the Mashpee 

Pitted Plain

Group 1

KAME_1,2 Nantucket ice-contact deposits Groups 1 and 2

FALMIC_1,2,3 Falmouth ice-contact deposits Groups 1, 2, 

and 3

BUZZMS_1,2,3 Southern part of the Buzzards 

Bay Moraine

Groups 1, 2, 

and 3

BUZZMN_1,2 Northern part of the Buzzards 

Bay Moraine

Groups 1 and 2

SAND_1,2 Western part of the Sandwich 

Moraine

Groups 1 and 2

BUZZP_1,2,3 Buzzards Bay Plain deposits Groups 1, 2, 

and 3

CCLAK_1,2 Cape Cod Bay lacustrine 

deposits

Groups 1 and 2

SAND2_1,2 Eastern part of the Sandwich 

Moraine

Groups 1 and 2

BARN_1,2 Barnstable Plain deposits Groups 1 and 2

HARW_1,2 Harwich Plain deposits Groups 1 and 2

MORCC_3 Moraine and lacustrine deposits Group 3

BHP_3 Barnstable and Harwich Plain 

deposits

Group 3

BASAL_4 Basal lacustrine deposits Group 4

COLL Collapsed outwash sediments Groups 1, 2, 3, 

and 4
14 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 6. Distribution of (A) hydraulic conductivity and boundary leakance parameter zones in the parameterized regional 

model (model 3, table 1), and (B) hydraulic conductivity in the trial-and-error regional model (model 1, table 1), western Cape 

Cod, Massachusetts.

MASH1_1

MASH3_1

MASH2_1

SAND2_1

SAND_1

BUZZMS_1

BUZZMN_1

HARW_1

BARN_1

SAND2_1

CCLAK_1

BUZZP_1

FALMIC_1

KAME_1

COAST2

COAST1

STREAM

BUZZM

Row 97

Column 107

Hydraulic conductivity—Group 1: altitude

greater than -10 feet

 

EXPLANATION

Modified hydraulic conductivity

(model 5, table 1)

BUZZMC_12

BUZZMC_12

Parameter type

Pond: K = 100,000 feet per day

Falmouth 

ice contact

deposits

210

230

250

220

240

100

100

300

300

300

200

240

300

300

350

130

200

Hydraulic conductivity from trial-and-error

model: altitude greater than 60 feet

Sandwich Moraine

Buzzards Bay Moraine

(South)

(North)

Mashpee

Pitted

Plain

Barnstable Plain Harwich

 Plain

Buzzards Bay Plain

deposits Cape

Bay lacustrine Cod

Falmouth 

ice-contact

deposits

A

B

Coastal boundary

Estuary boundary

Stream

Leakance

41o42'30"

41o

32'30"

70o37'30" 70o

17'30"

0

0

10 MILES

10 KILOMETERS
Numerical Ground-Water-Flow Model 15

Coastal

Leakance Hydraulic conductivity

Parameter zones represented in parameterized model

10–90

EXPLANATION

Hydraulic conductivity in trial-and-error model, in feet per day

100–190 200–290 300–350

Pond

100,000

Estuary

Column 107

Group 1 COAST1

MASH1_1

SAND_2

MASH1_3

MASH1_2

SAND_1

SAND_3

MASH3_1

MASH2_1

MASH2_2

MASH2_3

BASAL_4

COLL

COAST1

COAST2

COAST2

Group 4

Group 3

Group 2

 -10

-100

-240

BEDROCK

D

(each block represents a discrete value)

Column 107

 NGVD 29

-100

-240

65

-140

B

NGVD 29

-100

-240

65

-140

Row 97

A

Row 97

MASH1_1

MASH2_1

MASH1_3

MASH1_2

MASH3_1

SAND_123

BASAL_4

BUZZMN_1 BARN_1 HARW_1

BARN_2 HARW_2

BHP_3

COLL

COLL

COLL

COLL

MASH2_3

BUZZMN_2 MASH2_2

MORCC_3

BUZZP_1

BUZZP_2

COAST1 Group 1 COAST1

Group 4

Group 3

Group 2

 -10

-100

-240

BEDROCK

WATER TABLE

C

ALTITUDE, IN FEET

ALTITUDE, IN FEET

ALTITUDE, IN FEET

ALTITUDE, IN FEET

WATER TABLE

0

0

10 MILES

10 KILOMETERS

VERTICAL SCALE GREATLY EXAGGERATED

Figure 7. Vertical distribution of hydraulic conductivities in the trial-and-error regional model (model 1, table 1) along 

(A) row 97, and (B) column 107; and vertical distribution of hydraulic conductivity and boundary leakance parameter zones in 

the parameterized model (model 3, table 1) along (C) row 97, and (D) column 107.
16 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

incorporated into four vertical groups. Each group has a 

different spatial distribution of hydraulic conductivity parameters; group 1 (altitude greater than -10 ft) (figs. 6A and 7C–D) 

has 14 different hydraulic conductivity parameter zones, 

whereas group 4 (altitude less than -240 ft) (figs. 7C–D) consists of 2 zones. The parameter COLL, which is included in all 

four groups (table 2), represents sediments within collapse 

structures; these sediments generally occur beneath ponds. 

Hydraulic conductivities from the trial-and-error calibration 

for all model cells within each parameter zone were averaged 

to estimate initial parameter values.

In addition to hydraulic conductivity parameters, 

leakances at hydrologic boundaries also are represented as 

parameters. Leakances at freshwater streams are represented 

by parameter STREAM (fig. 6). Coastal boundary leakances 

are represented by parameters COAST1 and COAST2 

(fig. 6). Parameter COAST1 represents leakances in open 

coastal waters, whereas COAST2 represents leakances in 

inland estuaries. Recharge was included as a parameter; 

however, the parameter value was specified as the value 

(27.1 in/yr) resulting from trial-and-error calibration. A specified recharge value was used because inclusion of recharge as 

an estimated parameter resulted in non-convergence of initial 

parameter-estimation regressions.

The simulated water table produced from the 

parameterized model with initial values derived from the 

trial-and-error calibration was similar, but not identical, to 

that of the original trial-and-error model (fig. 5). The largest 

discrepancies (generally less than 1 ft) were near groundwater divides (or farthest from hydrologic boundaries) where 

hydraulic gradients are smallest. Bulk hydraulic conductivities 

are preserved in the parameterized model because the initial 

value for each parameter zone is the average of calibrated 

values for all model cells within the zone; however, the 

vertical simplification would be expected to result in a slightly 

different head distribution because flows through the aquifer 

to hydrologic boundaries, particularly in shallower parts of the 

aquifer, would differ. Simulated heads from both the trial-anderror and parameterized models reasonably match observed 

heads (table 1, fig. 8); the absolute mean residuals for the two 

models were 1.63 and 1.68 ft, respectively.

Incorporation of Observed Heads and Flows

The OBS Process was used to incorporate observations 

of heads and flows into the analysis of model calibration 

and uncertainty. This package allows the user to specify the 

values and locations of observations of hydrologic conditions 

within the model domain as well as the uncertainty associated 

with the observations. Sources of observation uncertainty, 

which are additive, could include measurement and survey 

errors. For observations representing estimates of steady-state 

hydrologic conditions, an additional component of uncertainty 

relates to how well the estimates represent actual long-term 

hydraulic conditions. Uncertainty is represented as a measure 

of spread, such as a variance, standard deviation, or coefficient 

of variation. 

Water-level data were assembled from four sources: 

(1) a network of 51 long-term observation wells jointly managed by USGS and the Cape Cod Commission, (2) synoptic 

measurements collected in March 1993, (3) measurements 

collected in and around the MMR in June–December 2000, 

and (4) measurements collected in and around Camp Edwards 

in June 2001. Historical data from long-term observation wells 

were used to estimate long-term mean water levels at the 51 

observation wells. The long-term data indicate that the synoptic measurements (March 1993, June–December 2000, and 

June 2001) were made during periods when conditions were 

near average. A total of 424 water-level observations—409 

wells and 15 ponds—were assembled and used for trial-anderror calibration during development of the regional model. A 

subset of 128 water-level observations, including observations 

from all 51 long-term observation wells, was used for purposes of the analyses of model calibration and uncertainty 

(fig. 9). This simplification was implemented to eliminate 

Figure 8. Absolute mean head residuals for the trial-anderror and parameterized regional models of western Cape 

Cod, Massachusetts. Parameter values derived from trialand-error calibration. Model numbers shown in table 1. ABSOLUTE MEAN HEAD RESIDUAL, IN FEET

MODEL

1.0

0.5

0

1.5

2.0

Trial-and-error model (1)

Linearized model (2)

Parameterized model (3)

 Parameterized model 

(groups 1 and 2 combined) (4)
Numerical Ground-Water-Flow Model 17

Figure 9. Locations of head and streamflow observations, Sagamore flow lens, western Cape Cod, Massachusetts.

EXPLANATION

11058837

11058837

Buzzards Bay

Nantucket Sound

Vineyard Sound

Cape Cod Bay Cape Cod Canal

Ground water

 Long-term observation well

 Synoptic: March 1993

 

 Synoptic: June–December 2000

 Synoptic: June 2001

 Pond

 Long-term streamflow-

 gaging station and number

 Partial-record streamflow-

 measurement site, May 2002

Impact Area

Massachusetts 

Military 

Reservation

Surface water

41o42'30"

41o

32'30"

70o37'30" 70o

17'30"

0

0

10 MILES

10 KILOMETERS
18 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

redundant data and allow for a more efficient parameterestimation solution. Water-level observations were eliminated 

when more than one observation applied to a single model 

cell or when a number of wells with similar water levels were 

located in close proximity. Because the aquifer is unconfined 

and vertical hydraulic gradients generally are small, multiple 

observations at different depths in the aquifer also were eliminated by including only the shallowest water-level observation 

(closest to the water table). 

Streamflows were measured at 19 sites on the Sagamore 

flow lens, including at one long-term streamflow-monitoring 

site at the Quashnet River (11058837) (fig. 9), in May 2002. 

Measurements at the long-term monitoring site and at a 

former long-term monitoring site on the Herring River on the 

adjacent Monomoy flow lens suggested that streamflows in 

May 2002 were close to long-term averages. All 19 observations of streamflow were included in the analysis of model 

calibration and uncertainty. 

The OBS Process allows for different weights, which 

are indicators of the degree of confidence in the observations, 

to be assigned to observations from different sources. 

Weights are calculated as the reciprocal of the measure of 

the variability of a population of observations expressed as 

a standard deviation, variance, or coefficient of variation. 

Streamflows can respond quickly to precipitation events 

and are variable over time. Single, partial-record streamflow 

observations may not be reliable measures of long-term 

flows, whereas data from long-term monitoring sites likely 

are more reasonable measures of long-term streamflows. 

As a result, the observation from the long-term observation 

site on the Quashnet River was given a higher weight than 

were observations from the 18 partial-record sites; a weight 

corresponding to a standard deviation (σ) of 1.0 ft3

/s was 

assigned, which indicates a high degree of confidence in 

the observation. The partial-record sites were assigned 

weights corresponding to coefficients of variation (CV) 

from 0.3 to 0.5, representing a low degree of confidence in 

the observations from the partial-record sites; the use of a 

CV causes the observation weight to be a function of the 

magnitude of the observation. For the analysis presented here, 

a uniform weight was used for all water-level observations. 

The uniform weight corresponds to a standard deviation of 

0.01 ft; this value represents a high degree of confidence in the 

water-level observations. 

This analysis is intended to illustrate an application of 

parameter estimation to the evaluation of model calibration 

and uncertainty. The analysis is, to some extent, specific to 

these weights. It should be noted that changing the weighting 

of observations would affect the parameter-estimation solution 

(Hill, 1998). An evaluation of the effect of varying observation 

weights on model calibration and uncertainty is not included 

in this evaluation. 

Use of Inverse Methods to Improve 

Model Calibration

The trial-and-error model was modified to allow for the 

use of parameter estimation to (1) improve model calibration, 

as defined by the statistical fit between observed heads and 

flows and simulated equivalents, and (2) estimate confidence 

intervals for estimated optimal parameters for use in evaluating model-prediction uncertainty. The process involves two 

general steps: (1) a formal analysis of sensitivities of simulated heads and flows to model parameters and (2) the use of 

those sensitivities to estimate optimal parameter values that 

best fit the observed data. These two steps use the SEN and 

PES Processes described in Hill and others (2000). 

Evaluation of Model Sensitivity

The SEN Process uses the sensitivity-equation method 

to calculate the sensitivities of heads and flows to each 

parameter (Hill, 1998) and produces observation sensitivities 

of model-calculated heads and flows at each observation 

specified in the OBS Process with respect to each parameter. 

Model-calculated observation sensitivities, represented as the 

derivative of head or flow with respect to a parameter (δhn

/δbl

or δQn

/δbl), are scaled; dimensionless sensitivities are scaled 

by multiplying the sensitivities by the product of the parameter 

value and observation weight, and 1-percent sensitivities 

are scaled by multiplying the calculated sensitivities by the 

parameter value divided by 100. Dimensionless sensitivities 

are useful for comparing sensitivities of observations with 

different units, such as heads and flows. The SEN Process also 

produces a composite-scaled sensitivity that is a measure of 

the overall sensitivity of each parameter to all observations. 

The composite-scaled sensitivity of a given parameter is the 

root of the mean squared sum of sensitivities at individual 

observations.

One-percent scaled sensitivities of head, as calculated by 

the SEN Process, approximate the change in head resulting 

from a 1-percent change in a given parameter. These sensitivities are independent of observation weights and can be 

computed for all model cells regardless of whether or not a 

cell contains an observation. For this reason, these sensitivities 

are referred to as grid sensitivities because they apply over the 

entire active model grid.

The SEN Process can be used with or without implementation of the PES Process; when used without PES, the 

SEN Process can be used as a diagnostic tool for a variety 

of analyses, including (1) identifying regions of the aquifer 

where additional observations would improve calibration, 

(2) identifying influential observations, or (3) evaluating the 

importance of different parameters to calibration at the locations of individual observations. Prior to the implementation 

of parameter estimation, the SEN Process can be used to 
Use of Inverse Methods to Improve Model Calibration 19

(1) determine which parameters can be adequately estimated 

from available observations and (2) identify observations that 

may not warrant inclusion in the parameter-estimation process. 

The SEN and OBS Processes were used to evaluate grid and observation sensitivities for the parameterized 

model (model 3, table 1). The model consists of one recharge 

parameter (aquifer recharge), 66 hydraulic parameters (33 

horizontal and 33 vertical hydraulic conductivities), and three 

boundary leakance parameters (streams, inland estuaries, and 

open coastal waters) (fig. 6). The parameter COLL, which 

represents sediments within geologic collapse structures, was 

not included in the analysis; these sediments are of limited 

extent and likely would not affect simulated heads and flows. 

Heads and flows were most sensitive to recharge; the composite sensitivity of heads and flows to recharge was about 6.5 

times larger than that of the next most sensitive parameter. The 

hydraulic conductivity parameters are associated with four 

separate vertical groups of parameter zones (fig. 7). Composite-scaled sensitivities of heads and flows at observations in 

the calibration set (fig. 9) to these parameters are illustrated 

in figure 10; the parameters are subdivided by parameter type 

and vertical group. The results show that simulated heads and 

flows generally are insensitive to vertical hydraulic conductivity. Simulated hydraulic conditions are most sensitive to horizontal hydraulic conductivity in the shallowest two parameter 

groups (groups 1 and 2); however, the composite sensitivity 

varies between individual parameters within each group 

(fig. 10A). As an example, composite sensitivities for horizontal hydraulic conductivity in the shallowest part of the aquifer 

(group 1: altitude greater than -10 ft) differ by more than two 

orders of magnitude (fig. 10B). For group 1, simulated heads 

and flows are most sensitive to the hydraulic conductivity of 

the western part of the Sandwich Moraine (SAND_1) and the 

northern part of the Buzzards Bay Moraine (BUZZMN_1); 

least sensitive in group 1 is the hydraulic conductivity of 

Harwich Plain deposits (HARW_1). If the model was to be 

calibrated by using a trial-and-error approach, the formal 

sensitivity analysis shows that the values of vertical hydraulic conductivity are not important for model calibration. The 

results also show that values of horizontal hydraulic conductivity in some areas of the shallow parts of the aquifer 

(groups 1 and 2) as well as boundary leakances are the most 

significant parameters to be considered during model calibration. The sensitivities suggest that, during trial-and-error 

calibration, little effort should be spent adjusting parameters 

with low sensitivities, and more effort should be spent adjusting parameters with higher sensitivities. 

One-percent sensitivities of heads are grid sensitivities 

that can be contoured to identify regions of the model most 

sensitive to a given parameter. Parameters representing the 

hydraulic properties in the Impact Area (fig. 1) include those 

associated with the Sandwich Moraine, the northern part of 

the Buzzards Bay Moraine, the northern part of the Mashpee 

Pitted Plain outwash, the Buzzards Bay Plain outwash, the 

Cape Cod Bay lacustrine deposits, and leakances at coastal 

boundaries (figs. 11A–F). In the central part of the Impact 

Area, the 1-percent sensitivity of head to the hydraulic 

conductivity of the northern part of the Mashpee Pitted 

Plain (MASH1_1), which underlies the Impact Area, is on 

the order of 0.03 ft (fig. 11A). The 1-percent sensitivities of 

heads in the same area to parameters SAND_1, BUZZMN_1, 

and BUZZP_1 are on the order of 0.06, 0.05, and 0.04 ft, 

respectively (figs. 11B, C, and D). The 1-percent sensitivity 

of heads to parameter CCLAK_1 is less than 0.01 (fig. 11E), 

and the sensitivity of heads to coastal leakances (combined 

sensitivity to COAST1 and COAST2) is from 0.01 to 0.02 ft 

(fig. 11F).

Both dimensionless and 1-percent sensitivities of heads to 

each parameter are calculated at each observation in the OBS 

Process. One-percent sensitivities at selected observations 

(fig. 12A) for the eight most sensitive parameters at each 

observation are shown in figures 12B–J. Observation sensitivities are useful in determining the parameters that would be 

most important in trial-and-error calibration at specific locations, as well as observations that likely are influential in the 

parameter-estimation process. The information also provides 

insight into model behavior. Simulated heads at observations 

in the northern part of the model are sensitive to a number 

of different parameters (figs. 12B–E). This result indicates 

that a number of different parameters can be adjusted during 

trial-and-error calibration to match observed heads in the areas 

around those observations. One-percent sensitivities exceed 

0.01 ft for all eight parameters at sites MW58S and MW47M2. 

The hydraulic conductivity parameters SAND, BUZZMN, 

BUZZP, CCLAK, and MASH1 and boundary leakance 

parameters (COAST1 and COAST2) are most significant in 

the northern part of the model. At some observations, a single 

parameter primarily governs simulated heads (figs. 12F and 

H). Simulated heads at observations in the southern part of 

the model generally are not as sensitive to model parameters 

as are heads in the northern part of the model (figs. 12F, G, 

H, and I). As an example, 1-percent sensitivities of simulated 

pond levels at Ashumet Pond do not exceed 0.02 ft for any 

parameter (fig. 12I). This result suggests that the simulated 

water level in the pond likely would not change significantly 

during trial-and-error calibration. Hydraulic conductivity 

parameters representing shallow parts of the aquifer (zones 

in groups 1 and 2) generally are most significant for model 

calibration for most observations. Streamflows at the Quashnet 

River (fig. 12J) are most sensitive to the hydraulic conductivity of the Mashpee outwash deposits (MASH2) and streambed 

leakance (STREAM).

In addition to their use as a diagnostic tool, the observation sensitivities produced by the SEN Process are used in 

the estimation of optimal parameters (PES Process). Prior to 

implementing the PES Process, the sensitivities can be used to 

decide which parameters likely can be estimated and, therefore, should be included in the estimation.
20 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts COMPOSITE-SCALED SENSITIVITY, DIMENSIONLESS

Composite-scaled sensitivity for 

aquifer recharge = 2,044

PARAMETER TYPE

0

100

200

300

400

SAND_1

BUZZMN_1

BUZZP_1

MASH1_1MASH2_1

FALMIC_1

MASH3_1

CCLAK_1

BUZZMS_1

BARN_1

KAME_1

SAND2_1

HARW_1

0

100

200

300

 HORIZONTAL HYDRAULIC CONDUCTIVITY PARAMETER

A All parameters

B Horizontal hydraulic conductivity (Group 1)

 Boundary leakance

Horizontal hydraulic conductivity (group altitude, in feet)

 Vertical hydraulic conductivity

Coastal (COAST1 and COAST2 combined)

Streams

Group 1 (> 60 to -10)

Group 4 (< -240)

Group 3 (-100 to -240)

Group 2 (-10 to -100)

All groups

Horizontal hydraulic conductivity Vertical hydraulic conductivity Boundary 

leakances

Figure 10. Composite-scaled 

sensitivities of heads to (A) all 

parameters (horizontal and vertical 

hydraulic conductivity and boundary 

leakances), and (B) horizontal hydraulic 

conductivity parameters for vertical 

group 1. Parameter values derived from 

trial-and-error calibration. Results are 

for model 3, table 1.
Use of Inverse Methods to Improve Model Calibration 21

41o42'30"

41o38'30"

41o42'30"

41o38'30"

70o37'30" 70o27'30" 70o37'30" 70o27'30" 70o37'30" 70o27'30"

00

5 MILES

5 KILOMETERS

LeakanceHydraulic conductivity

EXPLANATION

Parameter type

1-percent sensitivity

Line of equal 1-percentsensitivity—Contour intervalis 0.005 feet

BUZZP

MASH3

MASH2

SAND

BUZZMS

BUZZMN MASH1

FALMIC

0.01

0.105

D Buzzards Bay outwash (BUZZP_1)

BUZZP

MASH3

MASH2

SAND

BUZZMS

BUZZMN

CCLAK

CCLAK

MASH1

FALMIC

0.01

0.04

0.045

A Mashpee outwash (proximal) (MASH1_1)

BUZZP

MASH3

MASH2

SAND

BUZZMS

BUZZMN

CCLAK

MASH1

FALMIC

0.01

0.65

E Cape Cod Bay lacustrine (CCLAK_1)

BUZZP

MASH3

MASH2

BUZZMS

BUZZMN

CCLAK

MASH1

FALMIC

0.01

0.10

B Sandwich Moraine (SAND_1)

BUZZP

MASH3

MASH2

SAND

BUZZMS

BUZZMN

BUZZMN

CCLAK

MASH1

FALMIC

0.055

0.01

0.01

F Coastal leakances (COAST1 and COAST2)

 Camp Edwards Impact Area

 Camp Edwards Impact Area

BUZZP

MASH3

MASH2

SAND SAND

BUZZMS

CCLAK

MASH1

FALMIC

0.01

0.06

0.065

C Buzzards Bay Moraine (BUZZMN_1)

 Camp Edwards Impact Area

 Camp Edwards Impact Area Camp Edwards Impact Area

 Camp Edwards Impact Area

Coastal (COAST1)Estuary (COAST2)Stream

0.01

KAME KAME KAME

KAME KAME KAME

Figure 11. (A–F) The 1-percent scaled sensitivities of simulated heads to selected parameters (model 3, table 1), western Cape Cod, Massachusetts. Parameter values derived from trial-and-error calibration.
22 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 12. (A) Locations of selected observations on western Cape Cod, Massachusetts, and (B–J) the 1-percent scaled sensitivities for the observation locations. Parameter values derived from trial-and-error calibration. Results are for model 3, table 1.

Observation well and identifier

Streamflow site and identifier

Observation

41o42'30"

41o38'30"

70o37'30" 70o27'30"

BUZZP

MASH3

MASH2

SAND

BUZZMS

BUZZMN

CCLAK

MASH1

FALMIC

A

FSW320

FSW320

ECMP01D

MW65S

MW47M2

MW58S

MIW116

SDW260

AshumetPond

Quashnet River

Quashnet River

00

5 MILES

5 KILOMETERS

FALMIC_2

 BUZZMS_2

 COAST2

 BUZZP_2

 STREAM

FALMIC_1

 FALMIC_3

BUZZMS_1

0

0.070.060.050.040.030.020.01

G FSW320

BUZZP_2

COAST1

BUZZP_1

BUZZMN_1

COAST2

BUZZMN_2

STREAM

SAND_2

0

0.070.060.050.040.030.020.01

H ECMP01D

C MW65S

SAND_2

SAND_1

BUZZP_2

COAST1

BUZZP_1

COAST2

CCLAK_2

BUZZMN_2

0

0.070.060.050.040.030.020.01

MASH2_1

MASH1_2

BUZZMN_1

BUZZMN_2

BUZZP_2

FALMIC_2

MASH1_3

MASH2_2

0

0.070.060.050.040.030.020.01

PARAMETER PARAMETER PARAMETER

I Ashumet Pond

BUZZMN_1

BUZZP_2

BUZZMN_2

SAND_2

SAND_1

BUZZP_1

MASH1_2

MASH2_1

0

0.070.060.050.040.030.020.01

D MW47M2

0

750

3,0002,2501,500

J Quashnet River at Route 28

MASH2_1

 STREAM

MASH1_2

MASH1_3

 COAST2

MASH2_2

MASH2_2

MASH3_1

SAND_1

SAND_2

MASH1_1

MASH1_2

BUZZMN_1

BUZZP_2

BUZZMN_2

BUZZP_1

0

0.070.060.050.040.030.020.01

1-PERCENT SCALED SENSITIVITY,

IN FEET

B MW58S

CCLAK_2

SAND_2

SAND_1

CCLAK_1

COAST2

MORCC_3

BUZZMN_1

MASH2_1

0

0.070.060.050.040.030.020.01

1-PERCENT SCALED SENSITIVITY,

IN FEET

E SDW260

COAST2

STREAM

MASH2_1

MASH3_1

KAME_2

MASH2_3

MASH2_2

KAME_1

0

0.070.060.050.040.030.020.01

1-PERCENT SCALED SENSITIVITY,

IN FEET

F MIW116

Leakance Hydraulic

conductivity

EXPLANATION

Parameter type

CoastalEstuaryStream

1-PERCENT SCALED SENSITIVITY,

IN CUBIC FEET PER DAY

KAME
Use of Inverse Methods to Improve Model Calibration 23

Formulation of the Least-Squares Objective 

Function

The PES Process in MODFLOW-2000 uses nonlinear 

regression to estimate the optimal parameters that best fit the 

observed data. The fit between observed and simulated values 

is quantified by a weighted least-squares objective function. 

The objective function is expressed as

S b y y b P P b i

i

ND

i i p

p

NPR

( ) = − p p  ′( )  

 + −  ′( )  

 = =

∑ ∑   1

2

1

2

, (1)

where

b is a vector of each of the parameters being 

estimated,

ND is the number of observations,

yi is the ith observed value,

 y′

i is the simulated equivalent,

ωi is the weight of the ith observation,

NPR is the number of prior information estimates,

Pp is the pth prior information estimate,

P′

p is the simulated value, and

ωp is the weight of the pth prior information 

estimate (Hill, 1998).

The second term is a penalty term that represents prior 

information about aquifer characteristics; this information 

can be included in the analysis as an estimate of an aquifer 

property and a weight. When prior information is not included, 

the objective function includes only the first term, which is 

the sum of weighted squared residuals. The PES Process uses 

a modified Gauss-Newton optimization method to perform 

the nonlinear regression and estimate the optimal parameter 

values that minimize the objective function (Hill, 1998). This 

is an iterative method in which parameters are repeatedly 

updated and adjusted until a convergence criterion is met, 

indicating that the nonlinear regression has been completed 

and the objection function has been minimized.

Recharge was included as a parameter but was not 

estimated as part of the inverse calibration. Preliminary 

sensitivity analyses showed that simulated heads and flows 

were highly sensitive to recharge, which is the sole source 

of water to the aquifer; the composite-scaled sensitivity 

of recharge was 6.5 times larger than the sensitivity of 

the next most sensitive parameter. Preliminary parameter 

regressions showed that inclusion of recharge as a parameter 

made estimation of reasonable hydraulic conductivity and 

leakance parameters difficult. The recharge parameter value 

was specified as the value resulting from the trial-and-error 

calibration (27.1 in/yr).

Heads and Flows

The decision as to what parameters should be estimated 

is an iterative process in itself. A series of different regressions 

generally is necessary before a reasonable inverse model 

is produced. Common problems that limit the viability of a 

particular regression include (1) highly correlated parameters, 

(2) parameters with low sensitivities, and (3) estimated 

parameters that are unreasonable based on prior knowledge 

of the system. As discussed, recharge was excluded from the 

parameter-estimation regression because of nonconvergence. 

Vertical hydraulic conductivities also were not included 

because the low composite-scaled sensitivities (fig. 10) 

indicated that observation data were insufficient to estimate 

the parameters. All leakance and horizontal hydraulic 

conductivity parameters were initially included in the 

parameter-estimation regression. Initial parameter-estimation 

regressions showed that the only significantly high correlation 

coefficients (greater than 0.98) between parameters were 

between hydraulic conductivity parameters associated with 

vertical groups 1 and 2. A trial-and-error process yielded 

a viable regression in which the primary change from the 

original parameterization scheme (model 3, table 1) was the 

combination of the two vertical groups to prevent these high 

correlations between parameters associated with group 1 and 

the parameters associated with the underlying group 2. The 

high correlations arose because, in a given area of the model, 

a number of different combinations of hydraulic conductivity 

parameters in groups 1 and 2 would result in a similar 

overall transmissivity and a similar fit to observed data. This 

nonuniqueness can result in highly correlated parameters. 

In addition, the two coastal leakance parameters (COAST1 

and COAST2) were combined into a single parameter named 

ALL_COAST. The same spatial distribution of horizontal 

hydraulic conductivity parameter zones, reflecting the 

depositional model, was maintained. The composite-scaled 

sensitivities for the new configuration (model 4, table 1) are 

shown in figure 13A. As before, simulated heads and flows 

were insensitive to vertical hydraulic conductivity. Simulated 

heads generally were most sensitive to horizontal hydraulic 

conductivity, although sensitivities for individual horizontal 

hydraulic conductivity parameters ranged over nearly two 

orders of magnitude (fig. 13A). Of the 23 final parameters 

(21 hydraulic conductivity and 2 leakance parameters), 

15 were included in the regression (fig. 13B). Ten of the 

13 hydraulic conductivity parameters associated with the 

combined vertical groups 1 and 2 (altitude greater than -100 ft) 

were included in the regression; those that were not included 

represent hydraulic properties in the eastern part of the model, 

where observation data were limited and composite-scaled 

sensitivities were low (fig. 13B). Two of the six hydraulic 

conductivity parameters in group 3 and the single parameter 

in group 4 were included in the regression (fig. 13B); the 

remaining hydraulic conductivity parameters were excluded 

owing to low sensitivities that resulted in poor convergence of 

the regression.
24 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 13. Composite-scaled sensitivities of heads (for modified parameter distribution; model 4, 

table 1) to (A) all parameters, and (B) horizontal hydraulic conductivity and boundary leakance 

parameters. Parameters subsequently included in the parameter-estimation regression (estimated) or set 

at specified values (specified) are identified. Parameter values derived from trial-and-error calibration.

 Boundary leakance

 Horizontal hydraulic conductivity

Coastal (estimated)

Streams (estimated)

Groups 1 and 2 (estimated)

Group 4 (estimated)

Group 3 (estimated)

Group 3 (specified)

Groups 1 and 2 (specified)

 Boundary leakance

 Horizontal hydraulic conductivity

 (group altitude, in feet)

 Vertical hydraulic conductivity

Coastal

Streams

Groups 1 and 2 (>60 to -100)

Group 4 (<-240)

Group 3 (-100 to -240)

All groups

A 

B 

ALL_COAST

STREAM

SAND_12

BUZZP_12

BUZZMN_12

MASH2_12MASH1_12

FALMIC_12

CCLAK_12

MASH3_12

MASH1_3

BUZZMS_12

MORCC_3

BARN_12

FALMIC_3

KAME_12

BHP_3

BUZZP_3

MASH2_3

SAND2_12

BUZZMS_3

BASAL_4

COMPOSITE-SCALED SENSITIVITY, DIMENSIONLESS

0

100

200

300

400

500

600

HARW_12

PARAMETER TYPE

PARAMETER

Composite-scaled sensitivity for aquifer recharge = 2,044

0

100

200

300

400

500

600

Horizontal hydraulic conductivity Vertical hydraulic conductivity Boundary

 leakances
Use of Inverse Methods to Improve Model Calibration 25

Improvement in Calibration

Starting parameter values were determined by a spatial 

and vertical averaging of all model cells in each hydraulic 

conductivity parameter zone. The convergence criterion for 

the nonlinear regression was specified as a maximum fractional change in any parameter of 0.01. The objective function 

for the head observations decreased from 7.98 × 106

 to 

4.98 × 106

, indicating an improvement in calibration with 

respect to heads. The absolute mean of the unweighted head 

residuals—a traditional measure of model calibration—

changed from 1.73 ft, based on the initial parameter values 

derived from the trial-and-error model, to 1.41 ft, indicating 

an improvement of 0.32 ft (fig. 14). 

The Gauss-Newton procedure adjusts parameters to minimize the objective function, and the resulting optimal parameter values represent a statistical best fit to the observation 

data. As a result, the change in residuals at individual observations will vary, with residuals decreasing (improving) at some 

observations and increasing at others. The change in residuals represents the difference between unweighted residuals 

from the model that used parameter values derived from the 

trial-and-error calibration and the unweighted residuals from 

the model that used optimal parameter values. The changes in 

unweighted residuals ranged from -4.31 ft (improvement) 

to +2.95 ft; observations where head residuals improved 

outnumbered observations where head residuals worsened 

(fig. 15). The unweighted average residual change of 

-0.32 ft indicates an overall model improvement in head 

calibration. Changes in streamflow residuals ranged from 

-183,557 ft3

/d for the Coonamessett River (indicating 

improvement) to +83,705 ft3

/d for the Marston Mills River. 

The mean change in residuals was -17,660 ft3

/d (or about 

-0.2 ft3

/s), indicating an overall improvement in model fit to 

observed streamflows. The spatial distribution of the changes 

in residuals can identify areas where parameter estimation 

improved model fit over the trial-and-error model. The use 

of parameter estimation improved head calibration in most 

areas of the model (fig. 16). The magnitude of improvement in 

head residuals was largest in the Buzzards Bay and Sandwich 

Moraines and near the Buzzards Bay Plain (outwash) deposits; this is consistent with the composite-scaled sensitivities, 

which were highest for parameters representing the hydraulic 

conductivity of these sediments (fig. 13). Changes in residuals 

generally were smallest in the Mashpee Pitted Plain (outwash) 

deposits (fig. 16). 

Optimal Parameter Estimates

Initial hydraulic conductivity parameter values 

(derived from the trial-and-error model) and final optimal 

parameters estimated from this regression are shown in 

figure 17A. Changes in parameter values over the course of 

the parameter-estimation solution differed among parameters. 

Six of the 13 hydraulic conductivity parameters changed by 

less than 10 ft/d (fig. 17A); of the 7 remaining parameters, 

6 changed by 20 to 60 ft/d. One parameter (BUZZMS_12) 

increased by 444 ft/d to a final value of 638 ft/d (fig. 17A). 

Upper and lower reasonable limits for hydraulic conductivity 

parameters were assumed to be 350 and 10 ft/d, respectively; 

these values are based on aquifer tests, previous modeling 

efforts, and a general understanding of the water-transmitting 

properties of glacial sediments. 

The estimated value (638 ft/d) for parameter 

BUZZMS_12 exceeded the upper reasonable limit for hydraulic conductivity (350 ft/d) (fig. 17A); BUZZMS_12 represents 

the hydraulic conductivity of the southern part of the Buzzards Bay Moraine in the southwestern part of the model. This 

illustrates how model bias can be manifested in unreasonable 

estimated parameter values and how parameter estimation can 

yield diagnostic insight into a model beyond the estimation 

of optimal parameters. In this case, the unreasonable estimate 

indicates that, given observed heads in the area, a larger transmissivity is needed to fit those heads than can be obtained with 

values less than the upper reasonable limit (350 ft/d). This suggests that the simulated bedrock surface may be too shallow to 

allow for a sufficiently large saturated thickness and, therefore, transmissivity to match heads. Depth to bedrock in the Inverse calibration (variant 2)

MODEL (4) VERSION

1.0

0.5

0

1.5

2.0

Trial-and-errorderived parametervalues

Inverse calibration

ABSOLUTE MEAN HEAD RESIDUAL, IN FEET

Figure 14. Absolute mean head residuals 

calculated for the parameterized model (model 4, 

table 1) by using parameter values derived from the 

trial-and-error calibration and optimal parameters 

from the inverse calibration.
26 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 15. Magnitudes of changes calculated in residuals at individual observation locations for (A) heads 

and (B) streamflows for the parameterized model (model 4, table 1) by using parameter values derived from 

the trial-and-error calibration and optimal parameter values from the inverse calibration. CHANGE IN HEAD RESIDUAL, IN FEET Improvement

Average = -0.32 feet

INDIVIDUAL HEAD OBSERVATIONS (NOT IDENTIFIED)

-5

-4

-3

-2

-1

0

1

2

3

4

-2.0

-1.0

0

1.0

Coonamessett River

Phlashes Brook

Childs River

Marston Mills River

Bourne River

Little River

Quashnet River

Mashpee River

Santuit River

Hawes Run

Squeteague Brook

Pocasset River

STREAMFLOW OBSERVATION

CHANGE IN STREAMFLOW RESIDUAL, 

IN CUBIC FEET PER DAY (X 105

)

Backus River

Average = -17,660 cubic feet per day

Improvement

A

B
Use of Inverse Methods to Improve Model Calibration 27

Figure 16. Spatial distribution and magnitudes of changes calculated in absolute head residuals for the parameterized model 

(model 4, table 1) by using parameter values derived from the trial-and-error calibration and optimal parameter values from the 

inverse calibration, western Cape Cod, Massachusetts.

Impact

Area

 Cape Cod Bay

Buzzards Bay

41o

42'

30"

41o

32'

30"

A

B

70o37'30" 70o

27'30"

Base from U.S. Geological Survey digital data,

Universal Transverse Mercator projection,

Zone 19, 1:24,000, 1991

60

50

40

30

20

10

Particle track from

point A to point B

EXPLANATION

Simulated line of equal

water-table altitude—

Contour interval is 5 feet.

Datum is NGVD 29

Model (4) version

Mashpee outwash

Buzzards Bay Moraine

Sandwich Moraine

Buzzards Bay outwash

Cape Cod Bay

lacustrine deposits

Nantucket icecontact deposits

Trial-and-error-derived

parameter values

Inverse calibration

Inverse calibration

(variant 2)

Hydrogeologic unit

Location of head

observation

Change in residual—Blue

indicates improvement

+ / - 1 foot

+ / - 3 feet

+ / - 4 feet

+ / - 2 feet

20

60

50

40

30

20

10

20

20

Massachusetts Military

Reservation boundary B A

Falmouth ice-contact

deposits

Cape Cod Canal
28 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

EXPLANATION

EXPLANATION

0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

0

100

200

300

400

500

600

700

SAND_12

BUZZP_12

BUZZMN_12

MASH1_12MASH2_12

FALMIC_12

MASH3_12

CCLAK_12

MASH1_3

BUZZMS_12

MASH2_3

FALMIC_3

BASAL_4

PARAMETER

Upper reasonable 

limit = 350 ft/d

Specified value = 350 ft/d

(variant 2)

Lower reasonablelimit = 10 ft/d

Horizontal hydraulic 

conductivity

Leakance ALL_COAST STREAM

HORIZONTAL HYDRAULIC CONDUCTIVITY, IN FEET PER DAY

LEAKANCE, IN FEET PER DAY

Specified 

value = 110 ft/d

(variants 1 and 2)

840

0

100

200

300

400

500

600

700

0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

Upper reasonable

limit = 350 ft/d

Lower reasonable

limit = 10 ft/d

Horizontal hydraulic 

conductivity

Leakance

A

B

Parameter values

 Estimated value

 (model 4, table 1)

 Estimated value

 (model 4, variant 1)

 Estimated value

 (model 4, variant 2)

 Specified value

Parameter values

 Values derived from trial-

 and-error calibration 

 Estimated values

+1 Standard deviation

 Estimated value

-1 Standard deviation

Figure 17. (A) Initial and final estimated parameters and precisions of the estimates (as standard deviations) for parameters 

included in the parameter-estimation regression, and (B) estimated and specified parameter values for different parameter 

regressions for model 4 and variants 1 and 2 (table 1).
Use of Inverse Methods to Improve Model Calibration 29

area was determined on the basis of well borings completed 

through 2003; more recent drilling suggests that bedrock in 

the area is deeper than previously thought. Also, drilling logs 

from the area showed that coarse-grained sand and gravel 

extend to near bedrock in the area (Michael Goydas, Environmental Chemical Corporation, Inc., oral commun., 2005). 

The assumption inherent in the model is that sediments show 

a fining-downward sequence in which silty sediments underlie 

coarse-grained sediments. The presence of sandy sediments at 

depth would further increase aquifer transmissivity in the area, 

and the inherent assumption that sediments are finer grained 

at depth would contribute to an unrealistically high estimated 

hydraulic conductivity value for the shallow sediments. This 

suggests that the results of parameter estimation and field data 

both are consistent indicators that changes need to be made to 

the southwestern part of the model. 

Precision of Parameter Estimates

The precision to which the parameters were estimated, 

as indicated by the standard deviations of the estimates, also 

varied by parameter (fig. 17A). Standard deviations for the 

hydraulic conductivity parameters ranged from about 6 ft/d 

for parameter SAND_12 to 342 ft/d for parameter FALMIC_3 

(fig. 17A). The precision and confidence of an estimated 

parameter generally are a function of the goodness of fit to 

observed data and the sensitivity of simulated heads at the 

observations to the parameter, as indicated by the compositescaled sensitivity (fig. 13B). A large composite-scaled sensitivity indicates that the observations provide sufficient information to reliably estimate a parameter; conversely, a small 

sensitivity likely would result in a poorly estimated parameter 

with a large standard deviation. The coefficient of variation, 

which is the standard deviation divided by the optimal parameter value, is a dimensionless measure of precision and can be 

used to compare the precisions of different types of estimated 

parameters, such as hydraulic conductivity and leakance. The 

coefficient of variation was less than 0.1 for three parameters: 

BUZZMN_12, BUZZP_12, and SAND_12. These parameters had the largest composite-scaled sensitivities (fig. 13B), 

suggesting that parameter values could be reliably estimated 

from available observation data. Conversely, coefficients of 

variation exceeded 1 (maximum value of 8.7), indicating a 

standard deviation larger than the estimated value, for five 

parameters: FALMIC_12, MASH1_3, MASH2_3, BASAL_4, 

and FALMIC_3. The latter three parameters had the lowest 

composite-scaled sensitivities of all estimated hydraulic conductivity parameters (fig. 13B).

An important consideration is the precision of individual 

parameter estimates that have a significant effect on simulated 

conditions in areas of the model where predictions will be 

made. As an example, one purpose of the parameter-estimation 

process is to quantify uncertainties associated with particle 

tracks near the Impact Area in the northwestern part of the 

model (figs. 2 and 16). The uncertainty of simulated particle 

positions in the area would be a function of the precision of 

parameters that most affect those simulated positions. The 

precision of estimated parameters for which particle positions 

generally have small sensitivities would not be important 

in the evaluation of prediction uncertainty in those areas. 

Conversely, the precision of parameters that have an important 

control on particle positions in the area would be important in 

the evaluation of prediction uncertainties. Higher precisions 

(smaller standard deviations and coefficients of variation) 

generated by parameter-estimation regressions for parameters 

with larger sensitivities to the predictions of interest are favorable, whereas more imprecise estimates of parameters with 

small prediction sensitivities are of less consequence. 

An example particle track is shown in figure 16; the 

particle track is from the water table in the CIA (point A) to 

the discharge location at the Cape Cod Canal (point B). 

One-percent prediction sensitivities of the final particle position (point B) in the X-Y (horizontal) plane are shown in 

figure 18. The sensitivities in the X-Y plane represent the 

vector sum of the separate X (east-west) and Y (north-south) 

prediction sensitivities calculated by the SEN Process. 

The particle position is most sensitive to the two hydraulic 

conductivity parameters representing moraine sediments 

(SAND_12 and BUZZMN_12) and the parameter representing the hydraulic conductivity of the Buzzards Bay Plain 

deposits (BUZZP_12). The 1-percent sensitivities of the three 

parameters range from 53 to 94 ft, whereas sensitivities of the 

remaining 12 parameters are less than 14 ft (fig. 18). The three 

parameters with the largest sensitivities also are the parameters 

estimated with the highest precision; standard deviations and 

coefficients of variation of the three parameters were less than 

12 ft and 0.1, respectively (fig. 18). The largest measures of 

variability—indicating low precision—are parameters with 

low prediction sensitivities. This result indicates that the low 

precision for some estimated parameters, primarily those 

characterizing the southern part of the model, likely do not 

adversely affect estimates of prediction uncertainty in the area 

around Camp Edwards. 

Modification to the Model Framework

Adjustments were made to the model (model 4, table 1) 

to account for potential bias suggested by the unrealistically 

high estimate for the hydraulic conductivity of parameter 

BUZZMS_12 in the southwestern part of the model. The estimated value of 638 ft/d exceeded the upper reasonable limit of 

350 ft/d. The purpose of this adjustment was to (1) determine 

if a refined conceptual model could improve the estimated 

value of the parameter and (2) evaluate if the model bias in the 

southwestern part of the model would greatly affect estimated 

parameters in the area of interest (Camp Edwards and environs), which is located in the northwestern part of the model. 

Two additional regressions were performed: (1) one in which 

the BUZZMS_12 parameter zone was extended to bedrock 

and estimated along with the other parameters and (2) one in 

which the hydraulic conductivity of parameter BUZZMS_12 

was specified as the upper reasonable limit (350 ft/d) with the 
30 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

remaining parameters estimated as before. The two modified 

regressions are referred to as variant 1 and variant 2, respectively (table 1). In light of more recent drilling, which suggests 

that the aquifer may be more homogenous with depth than was 

previously thought, variant 1 may represent a better conceptual model of the hydrogeologic framework in that area of the 

model. Parameter FALMIC_3, which represents Falmouth 

ice-contact deposits (fig. 6) in group 3 (altitude of -100 to 

-240 ft), was removed from both regressions and specified as 

the value from the trial-and-error calibration (110 ft/d). The 

large standard deviation of about 480 ft/d (fig. 17) indicated 

that the parameter was poorly estimated. 

The estimated optimal parameters were different for the 

two regressions (fig. 17B). For variant 1, the BUZZMS_12 

parameter decreased to 484 ft/d (fig. 17B). This result indicates that extending the parameter zone to bedrock (creating a 

homogenous aquifer in the area) resulted in an estimate closer 

to the upper reasonable limit but still higher by more than 

100 ft/d. The only other parameter affected by the modified parameter configuration was FALMIC_12; the zone for 

FALMIC_12 is adjacent to the zone for BUZZMS_12 

(fig. 6). The parameter value decreased by 30 ft/d. This is 

likely also due, in part, to removal of parameter FALMIC_3 

from the regression and the use of a specified value of 110 ft/d 

for the zone represented by this parameter; this value represents an increase of about 60 ft/d over the original estimated 

value (fig. 17B). The remaining hydraulic conductivity parameters changed by less than 5 ft/d and the leakance parameters 

changed by less than 0.05 ft/d (fig. 17B). 

For variant 2, parameter BUZZMS_12 was specified 

as the upper reasonable limit (350 ft/d). As in variant 1, the 

only parameter affected by the change was FALMIC_12. The 

optimal parameter estimate decreased by 61 ft/d. This change 

is closely correlated with the effective increase of 62 ft/d for 

parameter FALMIC_3 that resulted from the use of a specified 

value of 110 ft/d for the parameter; a correlation coefficient of 

0

10

20

30

40

50

60

70

80

90

100

0

 40

 80

120

160

200

240

280

320

360

400

0

1

2

3

4

5

6

7

8

9

10

Coefficient of variation

X-Y plane

Standard deviation

1-percent sensitivity

Precision of estimate

SAND_12

BUZZP_12

BUZZMN_12

MASH2_12MASH1_12

FALMIC_12

CCLAK_12

MASH3_12

MASH1_3

BUZZMS_12

MASH2_3

FALMIC_3

BASAL_4

PARAMETER

1-PERCENT SENSITIVITY OF X-Y PARTICLE POSITION, IN FEET

ALL_COAST

STREAM

STANDARD DEVIATION, IN FEET

COEFFICIENT OF VARIATION, DIMENSIONLESS

Figure 18. The 1-percent sensitivity of advective-transport predictions in the central part of the Impact Area to estimated 

parameters and precision of the parameter estimates expressed as standard deviations and coefficients of variation. 

Predictions made with variant 2 of model 4 (table 1).
Use of Inverse Methods to Improve Model Calibration 31

0.97 indicates that the two parameters were highly correlated. 

The remaining hydraulic conductivity parameters changed by 

less than 12 ft/d and the leakance parameters changed by less 

than 0.1 ft/d (fig. 17B). The absolute mean head residual for 

variant 2 was about 1.5 ft (fig. 14), indicating a reasonable 

match to the observed heads. It was decided that variant 2 best 

represents the flow system because no hydraulic conductivity 

values exceed the upper reasonable limit.

The parameters of most importance for predictions of 

particle positions near Camp Edwards (parameters SAND_12, 

BUZZMN_12, and BUZZP_12) (fig. 18) generally were 

unaffected by changes to the southwestern part of the model. 

Except for the parameters associated with the southwestern 

part of the model (FALMIC_12, BUZZMS_12, FALMIC_3), 

changes in hydraulic conductivity and leakance parameters 

between the original inversely calibrated model (model 4, 

table 1) and variant 2 were less than 10 and 0.1 ft/d, respectively (fig. 17B). This result suggests that the potential model 

bias evidenced by the unreasonable parameter estimate likely 

does not adversely affect the use of the parameter estimates 

and precisions to evaluate particle-tracking uncertainty in the 

northwestern part of the model.

Advective-Transport Observations

Observed heads and streamflows are good indicators 

of hydrologic conditions in the aquifer at specific times; 

however, care must be taken when these observations are 

used to represent average (steady-state) conditions. Head and 

streamflow observations included in the OBS Process and the 

parameter-estimation regression either were from long-term 

monitoring locations or were measured during periods that 

were estimated to be at or near average hydraulic conditions. 

When available, contaminant plumes can be used as additional 

calibration targets. Plumes are good indicators of long-term 

hydraulic gradients in the aquifer because of the long time 

scale of advective transport in the aquifer as compared to the 

smaller time scale of changes in water levels and streamflows 

(Walter and Masterson, 2003). In the southern part of the 

MMR (fig. 1), contaminant plumes, which generally emanate 

from large, well-defined sources, have previously been used 

as qualitative calibration targets for trial-and-error models 

by tracking particles downgradient from source areas and 

manually adjusting model input parameters until the tracks 

reasonably matched plume center lines (Masterson and others, 

1997b; Masterson and Walter, 2000; Walter and Masterson, 

2003; Walter and Whealan, 2005). 

Observations of advective transport, such as plume paths, 

can be incorporated into parameter estimation by using the 

ADV2 Package (Anderman and Hill, 2001). Contaminant 

plumes are quantitatively useful for inverse model calibration 

if (1) the three-dimensional extent of contamination is well 

characterized by field sampling such that a general plume 

center line can be delineated, (2) the source area is well 

defined, and (3) the source history is well known such that 

traveltimes at one or more locations in the plume can be 

estimated. Subsurface contamination has been delineated 

in a number of areas on Camp Edwards, including the CIA, 

the J-Ranges, and Demo 1 (fig. 2). Of these areas, only 

contamination near Demo 1 fits the criteria for inclusion in 

the parameter-estimation regression. Contamination in the 

remaining areas has been used qualitatively to calibrate the 

trial-and-error model to flow patterns in and around Camp 

Edwards, but was not included in the parameter-estimation 

regression because the subsurface contamination was not well 

delineated owing to the diffuse nature of the source areas, the 

source areas were widely distributed and poorly defined, or 

source histories were not sufficiently well understood. 

Demo 1 was used first as a small-arms training range and 

then as an ordnance training and disposal area between the 

mid-1970s and the late 1990s; these latter activities resulted in 

a plume of RDX in the underlying aquifer that extends about 

2,000 ft downgradient from Demo 1 (AMEC, 2005). The 

inclusion of the Demo 1 plume in the parameter-estimation 

regression provided an opportunity to demonstrate the utility 

of the inverse-modeling method because, unlike plumes 

in other areas of the model, the plume path downgradient 

from Demo 1 was difficult to match during trial-and-error 

calibration. The ADV2 Package allows for observations of 

advective transport to be included in the parameter-estimation 

regression as three-dimensional transport locations at specified 

transport times; an observation is specified as the location, in 

Cartesian space, within a model cell. The regional model is 

oriented along the north-south axis such that the X, Y, and Z 

components defined for the advective-transport observation 

correspond to the east-west, north-south, and vertical 

directions, respectively. Each component of the observation 

location is assigned a weight. For the purposes of this analysis, 

weights of 0.2 ft were assigned for all three components of the 

observation location; this represents a very high weight for the 

advective-transport observation. An estimate of traveltime to 

the observation point also is included in the definition of the 

advective-transport observation. 

Although the ADV2 Package allows for the incorporation 

of multiple observations along a plume path, the leading edge 

of the current (2003) plume (defined by an RDX concentration 

greater than 2 μg/L) was used as the only advective-transport 

observation (fig. 19); there was insufficient information to 

estimate traveltimes at other locations within the plume. It is 

assumed that advective transport is the primary component of 

contaminant transport in the aquifer; this likely is a reasonable assumption given the high ground-water velocities and 

the homogeneity of the sandy aquifer sediments (Walter and 

Whealan, 2005). Velocities are inversely proportional to the 

porosity of aquifer sediments. Davis (1969) reported a range 

of porosity between 25 and 40 percent for unconsolidated sand 

and gravel; Garabedian and others (1991) reported a porosity 

of 0.39 at a location near the southern boundary of the MMR. 

A porosity of 0.35 was used in this analysis. The altitude of the 

observation was -10 ft based on cross sections of the plume 

(AMEC, 2005). Based on the source history, the traveltime 
32 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 19. Spatial extent of the RDX plume downgradient from Demolition Area 1, the advective-transport observation location, 

and particle tracks and advective-transport positions predicted by the parameterized model (model 4, variant 2, table 1) by using 

parameter values derived from the trial-and-error calibration and optimal parameter values from the inverse calibration, western 

Cape Cod, Massachusetts.

0

0 300 METERS

41o

41'19"

41o

41'49"

Demolition 

Area 1 

Advective-transport 

observation

70o 70 32'41" o

34'03"

EXPLANATION

Particle track and traveltime marker

 Model (4) version

 Trial-and-error-derived parameter values

 Inverse calibration to heads and flows (variant 2)

 Advective-transport observation

 Starting location: center of

 Demolition Area 1

 Final location: RDX plume defined as area

 with concentration greater than 2 micrograms 

 per liter. Vertical location is 10 feet below NGVD 29

RDX concentration, in micrograms

per liter

 2–10

 

 Greater than 10

Hydrogeologic unit

 Mashpee outwash

 Buzzards Bay Moraine

Extent of RDX contamination 2003

1,000 FEET

15 years 25 years 20 years

15 years

25 years 20 years

Inverse calibration to heads and flows

Trial-and-error-derived values
Use of Inverse Methods to Improve Model Calibration 33

from the water table to the observation point was estimated 

to be between 15 and 25 years. The simulated particle track 

produced by the parameterized model using parameter values 

derived from trial-and-error calibration did not closely match 

the observed plume. The leading edge of the plume simulated 

with the trial-and-error parameter values for traveltimes of 

15, 20, and 25 years differs by about 1,650, 940, and 1,380 ft, 

respectively, from the leading edge defined by the advectivetransport observation (fig. 19). 

Whereas the original inverse model (model 4, variant 2, 

table 1) (without the advective-transport observation) had a 

better match to heads and flows than the same parameterized 

model using parameter values derived from the trial-and-error 

calibration, with an improvement in the absolute mean head 

residual of 0.32 ft, the inverse model yielded an inferior match 

to the leading edge of the plume defined from observation data 

(fig. 19). The differences in the positions of the leading edge 

of the plume simulated by using the optimal estimated parameters and the leading edge defined by observation data were 

about 1,820, 1,580, and 2,130 ft for transport times of 15, 20, 

and 25 years, respectively (fig. 19). This difference may have 

occurred because there were inconsistencies in the head observation data such that improvement in head calibration resulted 

in a less accurate match to the plume path. The inclusion of 

plume paths in the regression could provide information that 

might partially compensate for these inconsistencies or allow 

a more complex representation of the aquifer and, therefore, 

result in a better match to heads, flows, and plume paths.

Results showed that estimated parameters did not differ 

between representing the northwestern part of the model by 

the original (model 4) and modified regressions (variant 2) 

(fig. 17B). Therefore, the modified regression was used to 

evaluate the incorporation of advective-transport observations 

into the parameter-estimation process; this regression was 

used because all estimated parameter values were within the 

upper and lower reasonable limits. The 1-percent sensitivities 

of the X, Y, and Z particle location to the parameters for 

variant 2 are shown in figures 20A–C. The four most 

important parameters for the X-Y particle location, as 

evidenced by sensitivities, are the hydraulic conductivities of 

the Buzzards Bay and Sandwich Moraines (BUZZMN_12 and 

SAND_12), the Buzzards Bay Plain (BUZZP_12), and the 

central part of the Mashpee Pitted Plain (MASH2_12). The 

east-west (X) particle position (fig. 20A) generally is more 

sensitive to model input parameters than is the north-south 

(Y) position (fig. 20B). The plume path was difficult to match 

during trial-and-error calibration; the trial-and-error model 

yielded a plume center line that was to the north of the plume. 

The locations of the parameter zones associated with the low 

sensitivities of the north-south (Y) particle position (fig. 20B) 

explain why it was difficult to shift particle tracks to the south 

by manually adjusting parameter values. Even though the 

particle path is located within or near the Mashpee Pitted Plain 

(MASH1_12), the northern part of the Buzzards Bay Moraine 

(BUZZMN_12), and the Buzzards Bay Plain (BUZZP_12), 

the north-south (Y) particle position is most sensitive to the 

hydraulic conductivity of the Sandwich Moraine (SAND_12). 

This relation may not be intuitively apparent and illustrates 

how the SEN Process can yield insight into the regional nature 

of the flow system and assist in trial-and-error calibration to 

observed plume paths whether or not parameter estimation 

will be used. It should be noted that the vertical hydraulic 

conductivities of two aquifer zones (BUZZMNV_1 and 

MASH1V_12) were 2 of the 10 most important parameters in 

determining the vertical position of the particle; this indicates 

that inclusion of different types of observations can change the 

importance of different types of parameters. 

When an advective-transport observation is included 

in the parameter-estimation regression, additional terms are 

added to the objective function (eq. 1) that represent the sum 

of the squared weighted residuals for observations in each of 

the X, Y, and Z directions. During parameter estimation, the 

PES Process estimates the parameter values that best fit the 

observations by minimizing the new objective function. For 

an assumed traveltime of 25 years, the changes in estimated 

parameter values varied by parameter (fig. 21A). The largest 

change was for parameter MASH1_12, which represents 

the northern part of the Mashpee Pitted Plain (fig. 6). The 

estimated parameter value increased by 375 to 538 ft/d, well 

above the upper reasonable limit of 350 ft/d (fig. 21A). The 

hydraulic conductivity of the Falmouth ice-contact deposits 

increased by 195 to 345 ft/d—close to the upper reasonable 

limit (fig. 21A). The hydraulic conductivity of deeper 

sediments (altitude of -100 to -240 ft) in the central part of 

the Mashpee Pitted Plain (MASH2_3) decreased by 35 ft/d to 

7 ft/d—close to but below the lower reasonable limit of 

10 ft/d (fig. 21A). Changes in the remaining estimated 

hydraulic conductivity parameters were less than 20 ft/d. 

The particle path produced by the new optimal parameters 

better matched the plume. With an assumed traveltime of 

25 years for the observation, the predicted location of the 

leading edge of the plume was about 160 ft from the observed 

location (fig. 22). The absolute mean head residual increased; 

however, from about 1.50 ft (for variant 2) without inclusion 

of the advective-transport observation to 2.13 ft with the 

advective-transport observation included; these results indicate 

an inferior match to heads (table 1; fig. 23). The additional 

regression did yield favorable results regarding the match to 

the observed plume; however, the regional calibration to heads 

worsened and the regression produced parameter estimates 

outside the reasonable range. This suggests that an adjustment 

to the conceptual model of the geologic framework could 

improve the match to both heads and flows.

A separate regression using a modified aquifer zonation was done to determine if a reasonable match to the 

plume could be obtained while preserving the match to 

observed heads and avoiding parameter estimates outside 

the range of reasonable values. An additional parameter was 

defined for the central part of the Buzzards Bay Moraine 

(BUZZMC_12) (fig. 6). The lithology of the moraines is more 

variable and poorly understood than that of the outwash and 

ice-contact deposits, and simulating the moraines, which are 
34 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 20. Sensitivity of the simulated three-dimensional particle position of the leading edge of a contaminant plume at 

Demolition Area 1 to hydraulic conductivity and leakance parameters for (A–C) original (model 4, variant 2, table 1) and (D–F)

modified (model 5, table 1) aquifer zonation. 1-PERCENT SENSITIVITY OF PARTICLE POSITION, IN FEET

PARAMETER

0

0.01

0.02

0.03

0.04

MASH1_12

MASH1_3

BUZZMN_12

BUZZP_12

SAND_12

BUZZMNV_12

MASH1V_12

MASH2_12

C Z (vertical) position MORCC_3

CCLAK_12

Upward shift

Downward shift

MORCC_3

0

20

40

60

80

SAND_12

BUZZMN_12

MASH2_12

BUZZP_12

ALL_COAST

MASH1_12

CCLAK_12

FALMIC_12

MASH1_3

B Y (north-south) position

Northward shift

Southward shift

0

20

40

60

80

BUZZP_12

SAND_12

MASH2_12MASH1_12

BUZZMN_12

MASH1_3

CCLAK_12

ALL_COAST

A X (east-west) position

MASH3_12

MORCC_3

A–C Original aquifer zonation

Eastward shift

Westward shift

SAND_12

BUZZP_12

BUZZMC_12

MASH1_12MASH2_12

BUZZMN_12

MASH1_3

CCLAK_12

ALL_COAST

MORCC_3

0

20

40

60

80

D–F Modified aquifer zonation

D X (east-west) position

Eastward shift

Westward shift

SAND_12

BUZZMN_12

BUZZMC_12

MASH2_12

BUZZP_12

ALL_COAST

MASH1_12

CCLAK_12

FALMIC_12

MASH1_3

0

20

40

60

80

E Y (north-south) position

Northward shift

Southward shift

MASH1_12

MASH1_3

BUZZMC_12

BUZZMN_12

BUZZP_12

SAND_12

BUZZMNV_12

MASH1V_12

MASH2_12

BUZZP_3

0

0.01

0.02

0.03

0.04

F Z (vertical) position

Upward shift

Downward shift
Use of Inverse Methods to Improve Model Calibration 35

0

100

200

300

400

500

600 EXPLANATION

Specified value

Estimated value

 Head and flow observations only 

 Advective-transport observation

 included. Traveltime = 25 years Upper reasonable 

limit = 350 ft/d

Lower reasonable

limit = 10 ft/d

Leakance

Specified 

value = 110 ft/d

Specified 

value = 350 ft/d

Upper reasonable 

limit = 350 ft/d

Lower reasonable

limit = 10 ft/d

Horizontal hydraulic conductivity

Leakance

Specified 

value = 110 ft/d

Specified 

value = 350 ft/d

[

A Original aquifer zonation

B Modified aquifer zonation

Horizontal hydraulic conductivity

0

100

200

300

400

500

600

SAND_12

BUZZP_12

BUZZMN_12

MASH1_12MASH2_12

FALMIC_12

MASH3_12

CCLAK_12

MASH1_3

BUZZMS_12

MASH2_3

FALMIC_3

BASAL_4

PARAMETER

ALL_COAST

STREAM

BUZZMC_12

BUZZMN_12

MASH1_12MASH2_12MASH3_12

FALMIC_12

BUZZMS_12

SAND_12

BUZZP_12

CCLAK_12

MASH1_3MASH2_3

FALMIC_3

BASAL_4

ALL_COAST

STREAM

0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

LEAKANCE, IN FEET PER DAY

0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

HORIZONTAL HYDRAULIC CONDUCTIVITY, IN FEET PER DAY

EXPLANATION

Traveltime = 15 years

Estimated value 

for original aquifer

zonation and 25-year

traveltime

Traveltime = 20 years

Specified value

Traveltime = 25 years

Estimated values for

modified aquifer zonation

Figure 21. Estimated and specified parameter values for regressions with and without inclusion of the advective-transport 

observation for (A) original aquifer zonation (model 4, variant 2, table 1) and (B) modified aquifer zonation (model 5, table 1) for various 

assumed traveltimes.
36 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 22. (A) Spatial extent of the RDX plume downgradient from Demolition Area 1 and particle tracks and advectivetransport positions predicted by the parameterized models (model 4, variant 2, and model 5; table 1) by using parameter values 

derived from the trial-and-error calibration and optimal parameter values from inverse calibration to heads, flows, and the 

advective-transport observation; and (B) vertical extent of the RDX plume and particle tracks from the inverse models for 

different assumed traveltimes to the advective-transport observation location, western Cape Cod, Massachusetts.

VERTICAL EXAGGERATION X 6

Advective transport observation Advective-transport observation

NGVD 29

 -40

-153

-140

 62

-80

-100

-20

20

40

-60

B

41o41'19"

41o

41'49"

70o 70 32'41" o

34'03"

0

0 300 METERS

1,000 FEET

Approximate vertical extent of RDX contamination

Demolition

 Area 1

 

Advective-transport

observation

EXPLANATION

Particle track and traveltime marker—Original

aquifer zonation (model 4, variant 2, table 1)

 Trial-and-error-derived parameter values

 Inverse calibration to heads and flows

 Inverse calibration to heads, flows,

 and advective-transport observation.

 Traveltime = 25 years

RDX concentration greater than

2 micrograms per liter

Hydrogeologic unit

 Mashpee outwash

 Buzzards Bay Moraine

Particle track and traveltime marker—

Modified aquifer zonation (model 5, table 1)

with inverse model calibrated to heads, flows, and

advective-transport observation

 Traveltime = 25 years

 Traveltime = 20 years

 Traveltime = 15 years

Advective-transport observation

 Starting location: center of Demolition Area 1

 Final location: RDX plume defined as area with 

 concentration greater than 2 micrograms per liter.

 Vertical location is -10 feet. Datum is NGVD 29

15 years

25 years

(original zonation)

25 years

15 years

25 years

20 years

20 years

25 years

20 years

15 years

A

Extent of RDX contamination

0

0 300 METERS

1,000 FEET

ALTITUDE, IN FEET

Inverse calibration to heads and flows

Trial-and-error-derived values

Inverse calibration to heads, flows, and advective-transport observation
Use of Inverse Methods to Improve Model Calibration 37

Figure 23. Absolute mean head residuals for versions of the 

inversely calibrated regional model with and without inclusion 

of the advective-transport observation and with different 

assumed traveltimes to the observation location. ABSOLUTE MEAN HEAD RESIDUAL, IN FEET

PARAMETER-ESTIMATION VERSION

1.0

0.5

0

1.5

2.0

2.5

 Advective-transport

observation included

(25-year traveltime) 

Head and flow 

observations

25-year traveltime

Advective-transport observation

included 20-year traveltime

15-year traveltime

Modified aquifer

zonation (model 5, table 1)

Original aquifer zonation

(model 4, variant 2, table 1)

(2) the ability of a model to produce simulated equivalents that 

closely match observed hydrologic conditions. The inverse 

methods discussed here address primarily the latter component 

of uncertainty. As was illustrated by the addition of a new 

hydraulic conductivity zone to the model, changing the definition of the aquifer zonation, or conceptualization, changes 

parameter sensitivities which, in turn, change estimates of 

optimal parameter values and precisions. Therefore, estimates 

of optimal parameters and precisions are specific to a given 

model design and will vary if the hydrogeologic framework 

or boundaries are changed. This consideration illustrates the 

importance of understanding the hydrogeology of the system being investigated in the proper application of inverse 

methods to model calibration. It should be remembered that 

components of uncertainty arise from inevitable discrepancies 

between the real aquifer system and the numerical model that 

represents the aquifer.

hydraulically important as shown by the parameter sensitivities (figs. 10 and 11), as single units may not be valid. Introducing the additional parameter changes both the magnitude 

of parameter sensitivities as well as the order of important 

parameters (figs. 20D–F). The sensitivities of the east-west 

(X) particle position are most affected by the modified aquifer 

zonation; the magnitude of sensitivities of the north-south (Y) 

position generally did not change. The additional parameter 

(BUZZMC_12) had the second and third highest sensitivities 

for the X and Y particle positions, respectively. This illustrates 

the importance of the conceptualization of the aquifer system 

and how it can affect the parameter-estimation regression, and 

underscores the need to understand regional geology when 

inverse calibration methods are used.

Parameter-estimation regressions that incorporated the 

modified aquifer zonation were run with assumed traveltimes 

of 15, 20, and 25 years to the advective-transport observation 

location. The resulting parameter estimates for the three 

regressions varied for each regression; however, all the 

optimized parameters were within the reasonable range of 

values (fig. 21B). The largest change in a parameter value 

was for parameter MASH1_12, which decreased from 

538 ft/d—above the upper reasonable limit—to between 202 

and 246 ft/d (fig. 21B). Particle tracks produced by using 

optimal parameters from the three regressions more closely 

matched the plume center line than did the particle track from 

the regression in which the advective-transport observation 

was included and the original aquifer zonation was used 

(fig. 22A). The model-predicted leading edge of the plume 

differed from the observed location by 160, 110, and 80 ft for 

assumed traveltimes of 15, 20, and 25 years, respectively 

(fig. 22A). The absolute mean residuals for heads were 1.81, 

1.75, and 1.58 ft for the regressions based on traveltimes of 

15, 20, and 25 years, respectively (table 1; fig. 23); both the 

match to the observed plume and the match to observed heads 

were better with the modified aquifer zonation. 

The vertical particle paths were affected by the estimated 

traveltime to the advective-transport observation (fig. 22B). 

Based on plume cross sections, the estimated altitude of 

the advective-transport observation was about -10 ft. The 

model-predicted altitudes of the leading edge of the plume for 

traveltimes of 15, 20, and 25 years were -12, -27, and -34 ft, 

respectively. A traveltime of 15 years yielded the best match 

to the vertical location of the leading edge of the plume; 

however, the absolute mean head residual was higher—

indicating an inferior match to heads—than mean head 

residuals for assumed traveltimes of 20 and 25 years (fig. 23).

Limitations of Analysis

There are a number of considerations when inverse methods are applied to model calibration and analyses of model 

uncertainty. Model uncertainties arise from a number of 

different sources, including (1) inaccuracies in the representation of the real aquifer system by the numerical model and 
38 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Observation weights also are important in estimating 

optimal parameter values and precisions. An observation’s 

weight, which is the inverse of the estimated precision of the 

observation, determines the importance of the observation 

residual in the objective function and the effect of the 

observation on the regression. The precision of a steady-state 

observation has several components, including measurement 

error, temporal variability, and the geographic accuracy 

of the observation location. For this analysis, heads at all 

observations were weighted equally. With the exception of 

streamflow at the long-term monitoring site, which was given 

a higher weight, flows from partial-record sites also were 

weighted equally. If heads and streamflows were differentially 

weighted based on the availability of additional information, 

then the parameter-estimation regression would yield different 

estimated values and precisions. The importance of the 

weighting of different types of observations is illustrated 

by inclusion of the advective-transport observation in the 

analysis. The results showed that improvement in calibration 

to the observed plume near Demo 1 resulted in a less accurate 

calibration to heads. Depending on how these different types 

of observations are weighted, different results can be obtained 

regarding the simulated match to these different types of 

observations. In this analysis, it was decided that the plume 

observation was a better indicator of long-term hydrologic 

conditions than were heads and flows; the advective-transport 

observation was given a higher weight, which resulted in a 

poorer match to the heads. 

The valid use of parameter estimation is based on 

the assumption that there is a reasonable match between 

observations and simulated equivalents and that observation 

errors are random and normally distributed (Hill, 1998). The 

absolute mean unweighted hydraulic-head residual for the 

optimized model (model 4, table 1) inversely calibrated to 

head and flow observations was about 1.4 ft, or about 

2 percent of the total head gradient for the model. This result 

indicates that, according to traditional measures of model 

calibration, the model is well calibrated. Weighted hydraulic 

head observations and simulated equivalents plot close to 

the 1:1 line (R2

=0.99), which also indicates a reasonable 

calibration (fig. 24A). A plot of weighted hydraulic head 

residuals and simulated values (fig. 24B) shows that weighted 

residuals are randomly distributed around a value of zero 

and that there are no discernable trends (R2

=0.00) suggesting 

model bias. The normality of the residuals can be evaluated 

qualitatively by a plot of weighted residuals and normal 

probability; the summary statistic RN

2

( )is calculated by the 

PES Process to assist in the evaluation of normality (Hill, 

1998, eq. 25). Hydraulic head residuals generally plotted along 

a straight line (indicating a normal distribution) (fig. 25); 

however, outliers in the upper and lower ranges plotted off the 

line (fig. 25). The value of the normal statistic RN

2

( ) was 0.94 

when all 128 observations were included in the regression. 

This value was less than the critical value of 0.98, indicating 

a greater than 5-percent chance that the head residuals are 

not normally distributed. To look at the effect of outlying 

observations on this measure of normality, a second regression 

was run in which eight head observations with large residuals 

(about 6 percent of the total) were removed. The decision of 

what constituted an outlying observation was determined by 

looking at major breaks in slope on figure 25. As before, the 

residuals plotted generally along a straight line and the normal 

statistic RN

2

( )increased to 0.99. Given the sample size, a value 

greater than 0.98 indicates a normal distribution at the 0.05 

significance level (Hill, 1998, table D1).

The application of parameter statistics to the evaluation 

of linear confidence intervals on the estimated parameters 

is most useful when models are linear. Most ground-water 

problems are inherently nonlinear, and linear confidence 

intervals for these problems are considered approximations of 

the actual confidence intervals (Hill, 1998). There is no direct 

measure of the linearity of confidence intervals; however, 

the Beale’s statistic (Hill, 1994) is used as a surrogate measure by quantifying the linearity of the confidence region of 

the parameters. The statistic can be computed by the utility 

program BEALE-2000 (Hill and others, 2000). A Beale’s 

statistic could not be computed for the estimated parameters 

because the utility program, which uses the variance-covariance matrix produced by the PES Process, yielded unrealistic hydraulic conductivity fields that could not be solved in 

forward runs; this result was due to the low precisions of some 

of the parameters, such as BUZZMS_12, FALMIC_12, and 

MASH1_3 (fig. 17A). These parameters were removed from 

the regression and specified as the optimal values, and the 

remaining parameters were log-transformed so that the utility 

program could be run. The Beale’s statistic for this run was 

246.1; therefore, the statistic of the original regression likely 

was greater than 246.1. This magnitude of the Beale’s statistic 

indicates that the model is highly nonlinear. This is a common 

situation when applying a linear approximation, such as the 

Gauss-Newton optimization method, to the solution of a nonlinear regression. Although the model as a whole is nonlinear, 

observations in the region of interest (Camp Edwards, in the 

northwestern part of the model) are sensitive to parameters 

that are well estimated in the regression as shown by small 

standard deviations, and the nonlinearity may affect prediction uncertainties less in this area than elsewhere (figs. 17 and 

18). Nevertheless, the nonlinearity of the model indicates that 

confidence intervals on each estimated parameter should be 

considered approximate.
Use of Inverse Methods to Improve Model Calibration 39

2,000 3,000 4,000 7,000

-800

-600

-400

-200

0

200

400

600

800

1,000 5,000 6,000

WEIGHTED OBSERVATION, DIMENSIONLESS

WEIGHTED SIMULATED EQUIVALENT, DIMENSIONLESS

0

0

1,000 2,000 3,000 4,000 5,000 6,000 7,000

0

1,000

2,000

3,000

4,000

5,000

6,000

7,000

WEIGHTED RESIDUAL, DIMENSIONLESS

WEIGHTED SIMULATED EQUIVALENT, DIMENSIONLESS

A

B

R2=0.99

R2

=0.00

Figure 24. Comparisons of (A) weighted observations and weighted simulated equivalents, and 

(B) weighted simulated equivalents and weighted residuals for hydraulic head observations. 

Results are for model 4, table 1.
40 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

particle tracking from known areas of contamination to 

identify downgradient drilling locations, and reverse tracking 

from subsurface contaminant detections to identify possible 

source areas (Walter and Masterson, 2003). Particle tracking 

also has been used to delineate recharge areas to production 

wells. Hypothetical forward and reverse particle tracks are 

shown in figure 26. A forward particle starts at a hypothetical 

source area at the water table (time=t0

) and reaches a downgradient point at time t6 (fig. 26A). A reverse particle track 

starts at a well (time=t0

) and terminates at a recharge location 

at the water table (time=tfinal) (fig. 26B). A confidence interval 

can be computed (by the utility program YCINT-2000) for a 

point in each of the X, Y, and Z directions from the parameter 

statistics and the sensitivities of the particle location to model 

parameters (fig. 26). Confidence intervals can be computed for 

any point along the forward particle path (t1

−t6) (fig. 26A) and 

at the endpoint of the reverse track (tfinal) (fig. 26B). A linear 

confidence interval is a function of the standard deviation 

computed for the model prediction. For a normal distribution 

with an estimated standard deviation, the regions within ±1, 2, 

and 3 standard deviations of the mean correspond to about 

68-, 95-, and 99-percent confidence intervals, respectively. 

Evaluation of Model-Prediction 

Uncertainty

Prediction uncertainties can be calculated from 

measures of precision for the estimated parameters. Linear 

confidence intervals on the parameters and the sensitivities 

of the predictions to those parameters are used to estimate 

confidence intervals for the predictions by using the utility 

program YCINT-2000 (Hill and others, 2000). Because the 

model is nonlinear, the resulting confidence intervals on 

model predictions are considered approximate. Confidence 

intervals can be evaluated for any type of model prediction, 

including heads, flows, and advective-transport positions. The 

model predictions do not have to be located in the same area 

as the observations; however, estimates of uncertainties likely 

are more reliable when predictions are made in regions with 

reliable observations. 

The model predictions of most interest at Camp Edwards 

are predictions of advective transport from particle-tracking 

analyses. Two different types of particle tracking have been 

used at the site to support data-collection efforts: forward 

Weighted residual 

(outlier)

Weighted residual

NORMAL ORDER STATISTIC

WEIGHTED RESIDUAL, DIMENSIONLESS

-3 -2 -1 0 1 2 3

-800

-600

-400

-200

0

200

400

600

800

Figure 25. Normal-order statistic and weighted residuals for hydraulic head observations. Results are for 

model 4, table 1.
Evaluation of Model-Prediction Uncertainty 41

Figure 26. Uncertainties associated with (A) forward and (B) reverse particle tracks represented as 

confidence intervals on particle positions.

Source

x

z

y

Flow line

t

1

t

2 t

3 t

4 t

5

t

0

t

6

X-Y plane

X-Z plane

A Forward particle track

Well Flow line

t

0

x

y

X-Y plane

B Reverse particle track

t

final

Recharge

location

X

Z

Direction of

Y

ground-water flow

Direction of

ground-water flow

Particle track (map view)

Starting position of particle

Final position of particle

Particle position at

 intermediate time

tfinal

t 0

t 4

Distance corresponding to

 plus or minus 2 standard

 deviations (σ) in the X, Y, 

 or Z direction

Ellipse of uncertainty

 corresponding to a 95-percent

 confidence interval around

 a particle position in the X-Y,

 X-Z, or Y-Z plane 

EXPLANATION

x
42 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Forward Particle Tracking

Forward particle tracking refers to the movement of 

particles through the model domain in the direction of groundwater flow and forward in time toward discharge locations; 

this methodology can provide predictions of advective 

transport downgradient from known source areas and is useful 

in determining possible drilling and subsurface sampling locations (Walter and Masterson, 2003). An example particle track 

from an area of known soil contamination within the CIA 

toward a discharge location near the Cape Cod Canal is shown 

in figure 27. The GIS post-processing software MODTOOLS 

(Orzol, 1997) was used to convert particle information from 

the particle-tracking program MODPATH for several points 

along the track into X and Y locations in geographic space 

and Z locations in altitude. The X and Y directions correspond to east-west and north-south axes, respectively, because 

the model grid is aligned with geographic north. The utility 

program YCINT-2000 was used to calculate standard deviations for each particle position along the track; these standard 

deviations, in turn, were used to estimate approximate 68-, 

95-, and 99-percent confidence intervals in the X and Y directions and 95-percent confidence intervals in the Z direction. 

These confidence intervals can be displayed graphically in the 

X-Y (geographic) and the X-Z or Y-Z (cross-section) planes as 

ellipses of uncertainty at each point (fig. 26) along the particle 

track. Uncertainty ellipses around multiple points along a particle track can be coalesced into swaths of uncertainty along 

the track in both the horizontal and vertical planes (figs. 27A 

and B). This approach is useful when source histories, and 

therefore traveltimes, are not well known.

The areas of uncertainty, defined by ellipses in the 

X-Y plane around each point along the particle track form a 

cone-shaped area of uncertainty along the particle track for 

each confidence interval (fig. 27A). The confidence intervals 

increase in width with increasing transport distance from the 

particle source. The 95-percent confidence intervals, which 

are common measures of uncertainty, in the X (east-west) 

and Y (north-south) directions were 46.2 and 44.7 ft wide, 

respectively, around the advective-transport point nearest 

the source and 2,083 and 1,970 ft wide, respectively, around 

the predicted discharge point (fig. 27A). The definition 

of uncertainties enhances the ability of decisionmakers to 

effectively use available resources when locating drilling sites 

to characterize subsurface contamination downgradient from 

known contaminant sources.

This increase in uncertainty along the particle track 

is, in part, a function of increasing parameter sensitivities 

with increasing transport distance. One-percent sensitivities, 

expressed as the vector sum of the sensitivities in the 

X (east-west) and Y (north-south) directions at three particle 

positions along the particle track (positions 1, 2, and 3 in 

fig. 27A) are shown in figure 28. The particle positions at the 

three locations are sensitive to the same parameters and in the 

same ascending order; however, the magnitudes of the sensitivities increase with increasing transport distance from the 

source (fig. 28). At the location closest to the source—about 

360 ft downgradient—the maximum 1-percent sensitivity 

(expressed as a vector sum of the sensitivities in the X and Y 

directions) was less than 5 ft. The maximum 1-percent sensitivities at distances of 10,860 and 19,240 ft downgradient from 

the source were 65 and 86 ft, respectively (fig. 28). Simulated 

particle positions near the source would not be expected to 

vary significantly with changes in model parameters because 

transport distances from the source are small. As transport 

distances increase, simulated particle positions (relative to the 

source) would vary more with changes in model parameters 

because the effects of the parameter values on the particle path 

are cumulative.

Uncertainties on vertical particle positions can be 

expressed as a band of uncertainty along the cross-sectional 

view of the particle path (fig. 27B). Like uncertainties in the 

X-Y plane, parameter sensitivities and uncertainties in vertical particle positions also increase with increasing transport 

distance. The 95-percent confidence intervals defined as 

distances above and below the simulated vertical particle positions, ranged from about 0.2 ft close to the source to 6.6 ft at 

the discharge location (fig. 27B). 
Evaluation of Model-Prediction Uncertainty 43

Figure 27. (A) Forward particle track and advective-transport positions near the Central Impact Area and 68-, 95-, and 

99-percent confidence intervals associated with the predictions; and (B) vertical forward particle track and vertical 95-percent 

confidence intervals associated with the prediction, western Cape Cod, Massachusetts. Results are for model 4, variant 2, table 1.

30

40

50

60

10

20

Camp Edwards 

Central Impact Area

 

20 years

40 years

30 years

60 years

50 years

10 years

NGVD 29

 -40

 62

-100

-123

A

B

Position 2

Position 1

Position 3

ALTITUDE, IN FEET 

1,000 METERS

0

0

5,000 FEET 

1,000 METERS

0

0

5,000 FEET 

70O36'42" 70O

31'54"

41O

42'

22"

41O

45'

01"

EXPLANATION

EXPLANATION

 Area of 68-percent confidence

 Area of 95-percent confidence

 Area of 99-percent confidence

 

Massachusetts Military Reservation

Impact Area

Area of 95-percent

confidence band

Simulated vertical

particle track

Upper limit of

vertical 95-percent 

confidence interval

Lower limit of

vertical 95-percent 

confidence interval

20

Cross section

((X,Y)-Z plane)

Simulated line of equal water-table altitude—

Contour interval is 5 feet. Datum is NGVD 29

Simulated particle track

Simulated advective-transport position

Advective-transport position shown in 

figure 28

Map view, X-Y plane

VERTICAL SCALE GREATLY EXAGGERATED

Areas of confidence
44 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 28. (A–C) The 1-percent sensitivity of simulated particle positions at three locations along the forward particle track 

to selected model parameters. Results are for model 4, variant 2, table 1. Positions shown in figure 27. 1-PERCENT SENSITIVITY, IN FEET

PARAMETER

BUZZP_12

SAND_12

MASH2_12MASH1_12

BUZZMN_12

MASH1_3

CCLAK_12

ALL_COAST

MASH3_12

FALMIC_12

0

20

40

60

80

100

BUZZP_12

SAND_12

MASH2_12MASH1_12

BUZZMN_12

MASH1_3

CCLAK_12

ALL_COAST

MASH3_12

FALMIC_12

BUZZP_12

SAND_12

MASH2_12MASH1_12

BUZZMN_12

MASH1_3

CCLAK_12

ALL_COAST

MASH3_12

FALMIC_12

A Advective-transport position 1 B Advective-transport position 2 C Advective-transport position 3

Reverse Particle Tracking

Reverse particle tracking refers to the movement of 

particles in the direction opposite of ground-water flow, 

or backward in time, towards recharge locations; reverse 

particle tracking can be used to identify recharge locations 

for subsurface contaminant detections and can assist in 

determining areas for soil sampling (Walter and Masterson, 

2003). Reverse particle tracks from several deep wells at 

Camp Edwards where RDX has been detected and the particletrack endpoints at the water table are shown in figure 29. Each 

discrete endpoint represents the simulated recharge location 

at the water table of an RDX-contaminated observation 

well. The uncertainty of the simulated recharge location can 

be represented as an ellipse, the axes of which represent 

confidence intervals in the X (east-west) and Y (north-south) 

directions (fig. 29). There is no meaningful uncertainty in the 

Z (vertical) direction because the point terminates at the water 

table. Instead of discrete points, the display of uncertainties 

allows for the definition of areas of confidence around the 

points (fig. 29). Approximate areas of 68-, 95-, and 99-percent 

confidence totaled 63.7 acres (0.10 mi2

), 218.6 acres 

(0.34 mi2

), and 396.6 acres (0.62 mi2), respectively 

(fig. 29). The use of these measures of uncertainty in planning 

prospective source-characterization activities allows for 

more confidence that locations of sampling sites adequately 

encompass potential areas of contamination.

Reverse particle tracking also can be used to identify 

areas at the water table that contribute water to production 

wells; this information can be used by decisionmakers to 

identify potential contaminant sources that could affect water 

quality in a well. A 100-year recharge area to production 

well #3 (fig. 2), located to the northeast of the Impact Area 

and pumping at a rate of about 0.56 Mgal/d, is shown in 

figure 30. The recharge area was produced by uniformly 

seeding particles at the water table and tracking the particles 

forward to discharge locations. The final locations of particles 

terminating in the production well were converted to starting 

locations, and the particles were tracked in the reverse 

direction to recharge locations at the water table. The 100-year 

recharge area, representing the area at the water table that 

contributes water to the well within 100 years of advective 

transport from this area, is about 262.4 acres (0.41 mi2

) 

(fig. 30). Uncertainties in the simulated recharge area can be 

represented as areas of confidence by merging uncertainty 

ellipses for the individual particles, similar to those shown 

in figure 29, into a single area. The resulting approximate 

areas of 68-, 95-, and 99-percent confidence totaled 

425.9 acres (0.66 mi2

), 657.4 acres (1.03 mi2), and 952.6 acres 

(1.49 mi2

), respectively (fig. 30). These areas of confidence 

represent increases of 61, 151, and 263 percent over the area 

encompassed by the simulated 100-year recharge area. Given 

a specified pumping and recharge rate, the recharge area does 

not change in size. As a result, the confidence areas can be 

thought to represent areas from which ground water could 

recharge the well with the specified confidence.
Evaluation of Model-Prediction Uncertainty 45

Figure 29. Simulated recharge locations for eight subsurface contaminant detections in and around the 

Central Impact Area and approximate areas of 68-, 95-, and 99-percent confidence associated with the 

predicted recharge locations, western Cape Cod, Massachusetts. Results are for model 4, variant 2, table 1.

EXPLANATION

Area of 68-percent confidence

Area of 95-percent confidence

Area of 99-percent confidence

Simulated reverse particle track

Simulated line of equal water-table altitude—

Contour interval is 5 feet. Datum is NGVD 29

Well with detection of Royal

Dutch Explosive (RDX) and identifier

Simulated recharge location

60

MW23

MW43

MW37

MW50

MW39

MW39

MW2

MW1

MW38

60

65

55

Impact Area 

Central Impact

Area

 

Sandwich

Bourne

70O 33'36" 70O

31'51"

41O 42'19"

41O 43'51"

1,000 2,000 FEET

300

0

0 600 METERS
46 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

Figure 30. Simulated recharge area to a production well near the Impact Area and approximate areas 

of 68-, 95-, and 99-percent confidence associated with the predicted recharge area, western Cape Cod, 

Massachusetts. Results are for model 4, variant 2, table 1. Location of production well #3 shown in figure 2.

Snake

Pond

Camp Edwards Impact

Area

Production well #3 

65

0

0 1,000 METERS

5,000 FEET

Massachusetts Military Reservation boundary

EXPLANATION

100-year recharge area

Areas of confidence

 Area of 68-percent confidence

 Area of 95-percent confidence

 Area of 99-percent confidence

Simulated line of equal watertable altitude—Contour interval 

is 5 feet. Datum is NGVD 29

65

70o31'11" 70o

29'04"

41o

41'31"

41o

43'33"
Summary and Conclusions 47

Summary and Conclusions

Camp Edwards, which is located in the northern part of 

the Massachusetts Military Reservation, has been a military 

training facility for more than 90 years, and historical training, weapons testing, and disposal activities have resulted in 

the contamination of soils and ground water in and around the 

facility. The underlying glacial aquifer, which consists primarily of sand and gravel, is the sole source of potable water for 

local communities. Given the large recharge rates and permeable sediments, contaminants can move by advective transport 

for large distances and potentially threaten water supplies. 

Contaminant sources, which include unexploded ordnance 

(UXO) and scattered munitions-disposal sites, often are difficult to fully characterize, and patterns of soil and ground-water 

contamination frequently resemble patterns consistent with 

nonpoint sources. Ground-water models, in conjunction with 

particle tracking, have been used by the Impact Area Groundwater Study Program (IAGWSP) at the site to determine 

(1) drilling locations and the placement of observation wells 

downgradient from known sources of soil and ground-water 

contamination, and (2) areas for soil sampling and source 

characterization by identifying potential source areas for subsurface contamination detected from reconnaissance drilling. 

The sporadic nature of some sources and the limited number 

of contaminant plumes as indicators of advective-transport 

patterns increase the reliance on these numerical tools to 

support field activities. 

Numerical models developed by different groups and calibrated by using trial-and-error methods often yield different 

predictions of advective transport even though the models are 

reasonably calibrated to the same hydrologic data. In addition, 

models calibrated by trial-and-error methods lack measurements of model uncertainty and do not provide quantitative 

information regarding the confidence in predictions of advective transport. The result is that decisionmakers may need to 

reconcile differing model predictions and plan field activities 

without quantitative information regarding prediction uncertainty that could allow for a more efficient use of resources. 

In 2003, the U.S. Geological Survey (USGS), in cooperation with the IAGWSP of the U.S. Army Environmental 

Command (USAEC), began an evaluation of model calibration and model-prediction uncertainty at Camp Edwards. 

The objective of the investigation was to evaluate the utility 

of inverse methods to improve model calibration and assess 

model-prediction uncertainty. An existing regional model of 

western Cape Cod was modified so that model inputs, including hydraulic conductivities and boundary leakances, were 

represented as parameters. The parameter-estimation capabilities integrated into MODFLOW-2000 were used in the 

analysis. The OBS Process was used to incorporate a network 

of 128 water-level and 19 streamflow observations, and the 

ADV2 Package was used to incorporate advective-transport 

observations into the model. The SEN Process was used to 

formally evaluate the sensitivity of simulated heads, flows, 

and advective transport to model input parameters and to 

provide sensitivities required for parameter estimation. The 

PES Process was used to estimate optimal parameter values 

that yielded the best statistical fit to the observed data; results 

from the PES Process were used in conjunction with the utility 

program YCINT-2000 to evaluate uncertainties associated 

with predictions of advective transport. The conclusions of the 

investigation are summarized below.

Formal sensitivity analyses using the SEN Process 

enhanced understanding of model sensitivities; the analyses 

are useful for both trial-and-error and inverse calibration methods. The formal sensitivity analyses showed that simulated 

heads and flows were most sensitive to recharge and horizontal hydraulic conductivity. It should be noted that recharge 

was not included as an estimated parameter because inclusion 

resulted in nonconvergence of initial parameter-estimation 

regressions. The parameters of most importance to model 

calibration were the hydraulic conductivities of the Buzzards 

Bay outwash plain and the Sandwich and Buzzards Moraines. 

Conversely, simulated heads and flows were not sensitive to 

vertical hydraulic conductivity and, therefore, those parameters were not important in model calibration. Simulated heads 

and flows were moderately sensitive to boundary leakances, 

particularly in areas near the boundaries. 

Parameter estimation (inverse calibration) using the PES 

Process, which uses nonlinear regression to minimize a leastsquares objective function, was useful in improving calibration 

to observed heads and streamflows. The absolute mean of the 

head residuals calculated by using optimal parameter values 

decreased 0.32 ft (from 1.73 to 1.41 ft) relative to the head 

residual from the parameterized model that used parameter 

values derived from the trial-and-error calibration; individual 

changes in heads ranged from about -4.3 (a value less than 0 

indicates improvement) to 3.0 ft. The absolute mean residual 

for streamflow improved by about 0.2 ft3

/s; individual changes 

in the residuals ranged from about -2.1 to 1.0 ft3

/s. The precision of the parameter estimates, represented as coefficients of 

variation, ranged from 8.7 to less than 0.1. Advective-transport 

predictions for Camp Edwards were most sensitive to the 

parameters with the highest precision (lowest coefficients of 

variation), which indicated that the model is well suited to 

address model uncertainties in the area.

The parameter-estimation process produced one 

estimated parameter value—the hydraulic conductivity of the 

southern Buzzards Bay Moraine deposits—that exceeded the 

upper reasonable limit of 350 ft/d, suggesting a model bias 

in the southwestern part of the model that may be related to 

an underestimate of bedrock depths in that area. This result 

illustrates that parameter estimation not only allows for 

estimates of parameters that give the best statistical fit to 

observation data, but also can identify potential biases in a 

model. The bias was not considered relevant to the analyses 

of ground-water flow around Camp Edwards, which is located 

in the northwestern part of the model, because (1) predictions 

of advective transport in and around Camp Edwards were 

insensitive to the parameter in question, and (2) alternate 
48 Inverse Modeling Methods for Calibration and Evaluation of Uncertainty, Cape Cod, Massachusetts

geologic representations in which the unreasonable estimate 

was reduced or eliminated did not affect parameter estimates 

around Camp Edwards. 

The incorporation of an advective-transport observation 

using the ADV2 Package allowed for an improved calibration 

to an observed plume path. The path of the contaminant 

plume downgradient from Demolition Area 1, in the southern 

part of Camp Edwards, had been difficult to match with 

the trial-and-error model. Including the leading edge of the 

plume as an observation produced a calibrated model that 

reasonably simulated the plume path; however, the absolute 

mean residual of heads increased to 2.13 ft, larger than that in 

the parameterized model using parameter values derived from 

the trial-and-error calibration. Modifying the representation 

of hydraulic conductivities in the model by including an 

additional hydrogeologic unit within the Buzzards Bay 

Moraine produced a model that matched the plume better than 

the previous inverse model and had an absolute mean residual 

of 1.58 ft, representing an improvement over the model using 

the trial-and-error-derived parameter values. 

Information from the PES Process was used to quantify 

uncertainties, expressed as about 68-, 95-, and 99-percent 

confidence intervals, on predictions of advective transport 

in and around the Impact Area. Confidence intervals on 

three-dimensional simulated particle positions are produced 

by the program YCINT-2000; separate confidence intervals 

are produced for the X (east-west), Y (north-south), and 

Z (vertical) directions. The confidence intervals can be 

graphically represented as ellipses around individual particle 

positions in the X-Y (geographic) plane and in the X-Z or 

Y-Z (vertical) planes. The merging of individual uncertainty 

ellipses allows forward particle tracks, which can be used 

to locate observation wells downgradient from known 

contaminant sources, to be displayed in map or crosssectional view as a cone of uncertainty around a particle 

path. Reverse particle tracks, which can be used to identify 

recharge locations for subsurface contaminant detections and 

contributing areas to production wells, can be geographically 

displayed as areas at the water table around the discrete 

particle endpoints.

Acknowledgments

The authors wish to thank the staff of the Impact Area 

Groundwater Study Program of the U.S. Army Environmental 

Command and the National Guard Bureau for their support 

and assistance during this investigation. Thanks are also 

extended to the staff of the Installation Restoration Program 

of the Air Force Center for Engineering and the Environment 

for providing data to assist in updating the regional model of 

western Cape Cod. Also, the authors thank Michael Goydas 

of Environmental Chemical Corporation, Inc., and the staff at 

AMEC Earth & Environmental, Inc., for facilitating 

the transmittal of lithologic and hydrologic data used in 

the investigation.

References Cited

AMEC Earth & Environmental, Inc. (AMEC), 2001, Central 

Impact Area groundwater report: Prepared for the National 

Guard Bureau Impact Area Groundwater Study Program, 

Technical Memorandum 01-6.

AMEC Earth & Environmental, Inc. (AMEC), 2002, J-1, J-3 

and L Ranges additional delineation report No. 1: Prepared 

for the National Guard Bureau Impact Area Groundwater 

Study Program, Report MMR-5834.

AMEC Earth & Environmental, Inc. (AMEC), 2005, 2005 

System performance and ecological impact monitoring 

(SPEIM) report, rapid response action systems Demo 1 

groundwater operable unit: Prepared for the National Guard 

Bureau Impact Area Groundwater Study Program, Report 

MMR-9842.

AMEC Earth & Environmental, Inc. (AMEC), 2006, Northwest Corner study area remedial investigation report: 

Prepared for the National Guard Bureau Impact Area 

Groundwater Study Program, Report MMR-9964.

Anderman, E.R., and Hill, M.C., 2001, MODFLOW-2000, 

The U.S. Geological Survey modular ground-water 

model—Documentation of the advective-transport observation (ADV2) package, version 2: U.S. Geological Survey 

Open-File Report 01–54, 67 p.

Barlow, P.M., and Dickerman, D.C., 2001, Numerical 

simulation and conjuctive-management models of the 

Hunt-Annaquatucket-Pettaquamscutt stream-aquifer system, 

Rhode Island: U.S. Geological Survey Professional Paper 

1636, 88 p.

Davis, S.N., 1969, Porosity and permeability of natural materials, in DeWiest, R.J.M., ed., Flow through porous media: 

New York, Academic Press, p. 54–89.

DeSimone, L.A., Walter, D.A., Eggleston, J.R., and Nimiroski, 

M.T., 2002, Simulation of ground-water flow and evaluation of water-management alternatives in the Upper Charles 

River Basin, eastern Massachusetts: U.S. Geological Survey Water-Resources Investigations Report 02–4234, 94 p.

Garabedian, S.P., LeBlanc, D.R., Gelhar, L.W., and Celia, 

M.A., 1991, Large-scale natural-gradient tracer test in 

sand and gravel, Cape Cod, Massachusetts—2. Analysis of 

spatial moments for a nonreactive tracer: Water Resources 

Research, v. 27, no. 5, p. 911–924.

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.
References Cited 49

Hill, M.C., 1994, Five computer programs for testing weighted 

residuals and calculating linear confidence and prediction intervals on results from the ground-water parameter 

estimation program MODFLOWP: U.S. Geological Survey 

Open-File Report 93–481, 81 p.

Hill, M.C., 1998, Methods and guidelines for effective model 

calibration: U.S. Geological Survey Water-Resources 

Investigations Report 98–4005, 91 p.

Hill, M.C., Banta, E.R., Harbaugh, A.W., and Anderman, 

E.R., 2000, MODFLOW-2000, The U.S. Geological Survey 

modular ground-water model—User guide to the observation, sensitivity, and parameter-estimation processes and 

three post-processing programs: U.S. Geological Survey 

Open-File Report 00–184, 209 p.

LeBlanc, D.R., 1984, Sewage plume in a sand and gravel 

aquifer, Cape Cod, Massachusetts: U.S. Geological Survey 

Water-Supply Paper 2218, 28 p.

LeBlanc, D.R., Garabedian, S.P., Hess, K.M., Gelhar, L.W., 

Quadri, R.D., Stollenwerk, K.G., and Wood, W.W., 1991, 

Large-scale natural-gradient tracer test in sand and gravel, 

Cape Cod, Massachusetts—l. Experimental design and 

tracer movement: Water Resources Research, v. 27, no. 5, 

p. 895–910.

LeBlanc, D.R., Guswa, J.H., Frimpter, M.H., and Londquist, 

C.J., 1986, Ground-water resources of Cape Cod, Massachusetts: U.S. Geological Survey Hydrologic-Investigations 

Atlas HA-692, 4 pl.

Masterson, J.P., and Barlow, P.M., 1997, Effects of simulated 

ground-water pumping and recharge on ground-water flow 

in Cape Cod, Martha’s Vineyard, and Nantucket Island 

Basins, Massachusetts: U.S. Geological Survey WaterSupply Paper 2447, 78 p.

Masterson, J.P., and Walter, D.A., 2000, Delineation of 

ground-water recharge areas, western Cape Cod, Massachusetts: U.S. Geological Survey Water-Resources Investigations Report 00–4000, 1 pl.

Masterson, J.P., Stone, B.D., Walter, D.A., and Savoie, 

Jennifer, 1997a, Hydrogeologic framework of western Cape 

Cod, Massachusetts: U.S. Geological Survey HydrologicInvestigations Atlas HA-741, 1 pl.

Masterson, J.P., Walter, D.A., and Savoie, Jennifer, 1997b, 

Use of particle tracking to improve numerical model calibration and to analyze ground-water flow and contaminant 

migration, Massachusetts Military Reservation, western 

Cape Cod, Massachusetts: U.S. Geological Survey WaterSupply Paper 2482, 50 p.

McDonald, M.G., and Harbaugh, A.W., 1988, A modular 

three-dimensional finite-difference ground-water-flow 

model: U.S. Geological Survey Techniques of WaterResources Investigations, book 6, chap. A1, 586 p.

Oldale, R.N., 1992, Cape Cod and the Islands—The geologic 

story: East Orleans, Massachusetts, Parnassus Imprints, 

205 p.

Oldale, R.N., and Barlow, R.A., 1986, Geologic map of 

Cape Cod and the Islands, Massachusetts: U.S. Geological 

Survey Miscellaneous Investigations Series Map I-1763, 

scale 1:100,000.

Oldale, R.N., and O’Hara, C.J., 1984, Glaciotectonic origin 

of the Massachusetts coastal end moraines and a fluctuating late Wisconsinan ice margin: Geological Society of 

America Bulletin, v. 95, p. 61–74.

Orzol, L.L., 1997, User’s guide for MODTOOLS—Computer 

programs for translating data of MODFLOW and 

MODPATH into geographic information systems files: 

U.S. Geological Survey Open-File Report 94–464, 89 p.

Pollock, D.W., 1994, User’s guide for MODPATH/

MODPATH-PLOT, version 3—A particle tracking postprocessing package for MODFLOW, the U.S. Geological 

Survey finite-difference ground-water flow model: 

U.S. Geological Survey Open-File Report 94–464, 234 p.

Prudic, D.E., 1989, Documentation of a computer program to 

simulate stream-aquifer relations using a modular, finitedifference, ground-water flow model: U.S. Geological 

Survey Open-File Report 88–729, 113 p.

Uchupi, Elazar, Giese, G.S., and Aubrey, D.G., 1996, The late 

Quaternary construction of Cape Cod, Massachusetts—

A reconsideration of the W.M. Davis model: Geological 

Society of America Special Paper 309, 69 p.

U.S. Army Environmental Command (USAEC), 2007a, 

Impact Area Groundwater Study Program, Overview and 

update—Transitioning to cleanup: Camp Edwards, 

Massachusetts Military Reservation, U.S. Army Environmental Command, Fact Sheet 2007–05, 8 p.

U.S. Army Environmental Command (USAEC), 2007b, 

Impact Area Groundwater Study Program, investigation and 

cleanup update: Camp Edwards, Massachusetts Military 

Reservation, U.S. Army Environmental Command, Fact 

Sheet 2007–05, 32 p.

Walter, D.A., and Masterson, J.P., 2003, Simulation of 

advective flow under steady-state and transient recharge 

conditions, Camp Edwards, Massachusetts Military Reservation, Cape Cod, Massachusetts: U.S. Geological Survey 

Water-Resources Investigations Report 03–4053, 51 p.

Walter, D.A., and Whealan, A.T., 2005, Simulated water 

sources and effects of pumping on surface and ground 

water, Sagamore and Monomoy flow lenses, Cape Cod, 

Massachusetts: U.S. Geological Survey Scientific 

Investigations Report 2004–5181, 85 p.

Walter, D.A., Masterson, J.P., and Hess, K.M., 2004, 

Ground-water recharge areas and travel times to pumped 

wells, ponds, streams, and coastal water bodies, Cape Cod, 

Massachusetts: U.S. Geological Survey Scientific Investigations Map I-2857, 1 pl.