Numerical Simulation of Ground-Water 

Flow Through . - 4^. Glacial ,, ,- . Deposits _£- and

Crystalline Bedrock in the Mteor Lake 

Area, ., Oraflon ...... , Count?, . • - ^/- -f Mew • • .... Mamt>shire -j|^. .....
NUMERICAL SIMULATION OF GROUND-WATER 

FLOW THROUGH GLACIAL DEPOSITS AND 

CRYSTALLINE BEDROCK IN THE MIRROR LAKE 

AREA, GRAFTON COUNTY, NEW HAMPSHIRE

By CLAIRE R. TIEDEMAN, DANIEL J. GOODE, and PAUL A. HSIEH

U.S. GEOLOGICAL SURVEY PROFESSIONAL PAPER 1572

1997
U.S. DEPARTMENT OF THE INTERIOR 

BRUCE BABBITT, Secretary

U.S. GEOLOGICAL SURVEY 

Gordon P. Eaton, Director

The use of firm, trade, and brand names in this report is for 

identification purposes only and does not constitute 

endorsement by the U.S. Government.

Library of Congress Cataloging-in-Publications Data

Tiedeman, Claire R.

Numerical simulation of ground-water flow through glacial deposits and crystalline bedrock in the Mirror 

Lake area, Grafton County, New Hampshire / by Claire R. Tiedeman, Daniel J. Goode, and Paul A. Hsieh.

p. cm -- (U.S. Geological Survey professional paper; 1572) 

Includes bibliographical references. 

Supt. of Docs, no.: I 19.16: 1572

1. Groundwater flow-New Hampshire-Mirror Lake (Grafton County) 

I. Goode, Daniel J. II. Hsieh, P.A. (Paul A.) III. Title. IV. Series. 

GB1197.7.T54 1997

551.49'09742'3-dc21 97-873

CIP

ISBN 0-607-86639-X

For sale by U.S. Geological Survey

Box 25286 

Denver, CO 80225-0286
CONTENTS

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

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

Background............................................................................... 2

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

Study Area................................................................................ 3

Physical Setting........................................................................ 3

Previous and Related Studies ................................................... 7

Acknowledgments.................................................................... 7

Ground-Water Hydrology............................................................... 7

Hydrogeologic Setting.............................................................. 7

Glacial Deposits............................................................... 8

Bedrock............................................................................ 10

Lake Sediments................................................................ 11

Ground-Water Flow.................................................................. 11

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

Discharge......................................................................... 12

Hydraulic Head......................................................................... 12

Data Collection................................................................ 12

Temporal Variation of Hydraulic Head............................ 16

Long-Term Average Hydraulic Head............................... 17

Model of Steady-State Ground-Water Flow................................... 17

Assumptions............................................................................. 17

Governing Equation.................................................................. 20

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

Boundary Conditions................................................................ 21

Parameter Structure.................................................................. 25

Model Calibration........................................................................... 28

Calibration Data........................................................................ 28

Nonlinear Least-Squares Regression Method.......................... 28

Optimal Estimates of Recharge from Precipitation

and Hydraulic Conductivity................................................ 29

Parameter Uncertainty.............................................................. 30

Goodness of Model Fit............................................................. 32

Simulation of Ground-Water Flow................................................. 35

Water Table............................................................................... 35

Mirror Lake Ground-Water Basin ............................................ 35

Ground-Water Budget............................................................... 41

Flow Paths Through the Ground-Water Basin.......................... 46

Summary........................................................................................ 47

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

III
IV CONTENTS

FIGURES

1-3. Maps showing:

1. Location, extent, and topography of study area............................................................................................... 4

2. Mirror Lake surface-water basin; sub-basins W, NW, and E; area of lacustrine mud sediments; and

area of sand and gravel deposits along southern shore of Mirror Lake........................................................... 5

3. Locations of flumes, piezometers, and bedrock wells as of 1994.................................................................... 6

4. Hydrogeologic section A-B-C-D ............................................................................................................................ 8

5. Map showing thickness of glacial deposits, area of lacustrine mud sediments, and areas of bedrock outcrop

androadcutin the vicinity of Mirror Lake.............................................................................................................. 9

6-7. Graphs showing:

6. Histogram of hydraulic conductivities of bedrock determined by single-well hydraulic tests ....................... 10

7. Hydraulic head from 1991 through 1993 in piezometer FS1-17 (A) and in packer-isolated intervals of

bedrock well FS3 (B)...................................................................................................................................... 16

8. Map showing long-term average altitude of water table in the vicinity of Mirror Lake......................................... 18

9. Section A-B showing long-term average hydraulic heads...................................................................................... 19

10. Map showing model grid and boundary conditions................................................................................................ 22

11. Section along row 25 of the model grid showing the five model layers................................................................. 24

12. Section along row 50 of the model grid showing glacial deposits, Mirror Lake,

lacustrine mud sediments, and bedrock in layers 1, 2, and 3.................................................................................. 24

13. Map showing delineation of the bedrock zone beneath upper hillsides and hilltops used in

parameter structure D.............................................................................................................................................. 27

14-15. Graphs showing:

14. Approximate, individual, 95-percent confidence intervals for (A) recharge from precipitation and

(B) hydraulic conductivities............................................................................................................................. 31

15. Normal probability plots of weighted residuals............................................................................................... 33

16. Maps showing distribution of hydraulic head residuals in model layers 1 through 5 in calibration A................... 36

17. Section A-B showing observed and simulated hydraulic head............................................................................... 38

18-20. Maps showing:

18. Water table for simulation A, and areas where water table lies in glacial deposits......................................... 39

19. Water table for simulation D, and areas where water table lies in glacial deposits......................................... 40

20. Horizontal extent of Mirror Lake ground-water basin in simulation A........................................................... 42

21. Diagram showing three-dimensional extent of Mirror Lake ground-water basin in simulation A......................... 43

22. Map showing horizontal extent of Mirror Lake ground-water basin in simulation D............................................ 44

23-24. Diagrams showing:

23. Three-dimensional extent of Mirror Lake ground-water basin in simulation D.............................................. 45

24. Components of ground-water budget for the Mirror Lake ground-water basin.............................................. 46

25. Schematic vertical section showing flow paths through the Mirror Lake ground-water basin............................... 47
CONTENTS

TABLES

1. Ranges of hydraulic conductivities of glacial deposits as determined by slug tests.................................................. 8

2. Long-term average streamflo w and baseflo w in Streams W, NW, and E at their inlets to Mirror Lake.................... 12

3. Altitude of screen midpoint, long-term average hydraulic head, and averaging method for piezometers................ 13

4. Altitude of top, bottom, and most transmissive fracture; long-term average hydraulic head; and

	averaging method for well intervals .......................................................................................................................... 14

5. Hydraulic conductivities used in parameter structures A, B, C, and D..................................................................... 25

6. Optimal estimates of recharge from precipitation and hydraulic conductivities....................................................... 29

7. Modified Beale's measure, RN2, and critical value of RN2 ........................................................................................ 32

8. Residual statistics used to measure overall goodness of model fit............................................................................ 34

9. Observed and simulated baseflow in Streams W, NW, and E.................................................................................... 34

10. Ground-water budget of Mirror Lake ground-water basin in simulations A and D.................................................. 45

CONVERSION FACTORS AND VERTICAL DATUM

Multiply

centimeter (cm)

centimeter per year (cm/year)

cubic meter per year (m3/year)

kilometer (km)

meter (m)

meter per second (m/s)

square kilometer (km2)

square meter (m2)

By

0.394

0.394

35.320

0.622

3.281

3.281

0.387

10.765

To obtain

inch

inch per year

cubic foot per year

mile

foot

foot per second

square mile

square foot

Temperature, given in degrees Celsius (°C), can be converted to degrees Fahrenheit (°F) by the following equation:

°F = 9/5 (°C) + 32

Sea level: In this report, "sea level" refers to the National Geodetic Vertical Datum of 1929—a geodetic datum 

derived from a general adjustment of the first-order level nets of both the United States and Canada, formerly called Sea 

Level Datum of 1929.

NUMERICAL SIMULATION OF GROUND-WATER FLOW

THROUGH GLACIAL DEPOSITS AND CRYSTALLINE

BEDROCK IN THE MIRROR LAKE AREA,

GRAFTON COUNTY, NEW HAMPSHIRE

By Claire R. Tiedeman, Daniel J. Goode, and Paul A. Hsieh

ABSTRACT

This report documents the development 

of a computer model to simulate steady-state 

(long-term average) flow of ground water in the 

vicinity of Mirror Lake, which lies at the eastern 

end of the Hubbard Brook valley in central New 

Hampshire. The 10-km2 study area includes 

Mirror Lake, the three streams that flow into 

Mirror Lake, Leeman's Brook, Paradise Brook, 

and parts of Hubbard Brook and the Pemigewasset River. The topography of the area is 

characterized by steep hillsides and relatively flat 

valleys. Major hydrogeologic units include 

glacial deposits, composed of till containing 

pockets of sand and gravel, and fractured crystalline bedrock, composed of schist intruded by 

granite, pegmatite, and lamprophyre. Ground 

water occurs in both the glacial deposits and 

bedrock. Precipitation and snowmelt infiltrate to 

the water table on the hillsides, flow downslope 

through the saturated glacial deposits and 

fractured bedrock, and discharge to streams and to 

Mirror Lake.

The model domain includes the glacial 

deposits, the uppermost 150 m of bedrock, Mirror 

Lake, the layer of organic sediments on the lake 

bottom, and streams and rivers within the study 

area. A streamflow routing package was included 

in the model to simulate baseflow in streams and 

interaction between streams and ground water. 

Recharge from precipitation is assumed to be 

areally uniform, and riparian evapotranspiration 

along stream banks is assumed negligible. The

spatial distribution of hydraulic conductivity is 

represented by dividing the model domain into 

several zones, each having uniform hydraulic 

properties. Local variations in recharge and 

hydraulic conductivities are ignored; therefore, 

the simulation results characterize the general 

ground-water system, not local details of groundwater movement.

The model was calibrated using a nonlinear 

regression method to match hydraulic heads 

measured in piezometers and wells, and baseflow 

in the three inlet streams to Mirror Lake. Model 

calibration indicates that recharge from precipitation to the water table is 26 to 28 cm/year. 

Hydraulic conductivities are 1.7 x 10~6 to 2.7 x 

10~6 m/s for glacial deposits, about 3 x 10~7 m/s 

for bedrock beneath lower hillsides and valleys, 

and about 6 x 10~8 m/s for bedrock beneath upper 

hillsides and hilltops. Analysis of parameter 

uncertainty indicates that the above values are 

well constrained, at least within the context of 

regression analysis. In the regression, several 

attributes of the ground-water flow model are 

assumed perfectly known. The hydraulic conductivity for bedrock beneath upper hillsides and 

hilltops was determined from few data, and 

additional data are needed to further confirm this 

result. Model fit was not improved by introducing 

a 10-to-l ratio of horizontal-to-vertical anisotropy 

in the hydraulic conductivity of the glacial 

deposits, or by varying hydraulic conductivity 

with depth in the modeled part (uppermost 150 m) 

of the bedrock.
NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

The calibrated model was used to delineate 

the Mirror Lake ground-water basin, defined as 

the volume of subsurface through which ground 

water flows from the water table to Mirror Lake 

or its inlet streams. Results indicate that Mirror 

Lake and its inlet streams drain an area of groundwater recharge that is about 1.5 times the area of 

the surface-water basin. The ground-water basin 

extends far up the hillside on the northwestern 

part of the study area. Ground water from this 

area flows at depth under Norris Brook to 

discharge into Mirror Lake or its inlet streams. 

As a result, the Mirror Lake ground-water basin 

extends beneath the adjacent ground-water basin 

that drains into Norris Brook.

Model simulation indicates that approxio

mately 300,000 m /year of precipitation recharges 

the Mirror Lake ground-water basin. About half 

the recharge enters the basin in areas where the 

simulated water table lies in glacial deposits; the 

other half enters the basin in areas where the 

simulated water table lies in bedrock. Flow from 

glacial deposits to bedrock occurs primarily 

beneath upper hillsides. Flow from bedrock to 

glacial deposits occurs primarily beneath lower 

hillsides and valleys. About 40 percent of the 

recharge (120,000 m3/year) moves through the 

basin along flow paths that remain in glacial 

deposits; about 60 percent (180,000 m3/year) 

moves along flow paths that pass through 

bedrock. Under steady state, total recharge to the 

basin equals total discharge from the basin. Of 

the total discharge, slightly more than half 

(170,000 m3/year) flows into the three Mirror 

Lake inlet streams; slightly less than half 

(130,000 m3/year) flows directly into Mirror 

Lake.

INTRODUCTION

BACKGROUND

The area in the vicinity of Mirror Lake, which 

lies at the eastern end of the Hubbard Brook valley in

central New Hampshire, has long been the subject of 

ecologic and hydrologic investigations. In 1955, the 

Forest Service of the U.S. Department of Agriculture 

established the Hubbard Brook Experimental Forest, 

which includes nearly the entire Hubbard Brook 

valley, as a center for watershed-management research 

in the northeastern United States. During the 1960's, 

the Forest Service together with several academic 

institutions developed a long-term, multidisciplinary 

study at the Experimental Forest to investigate the 

structure and function of northern hardwood forest 

ecosystems. Since the late 1970's, the U.S. Geological 

Survey has been conducting ground-water studies in 

the Mirror Lake area. This report presents findings 

from a part of the U.S. Geological Survey studies.

The initial study by the U.S. Geological Survey 

intended to characterize the hydrologic system 

associated with Mirror Lake, with emphasis on the 

interaction between the lake and ground water. 

Beginning in 1979, Winter (1984) led a long-term 

effort to monitor lake evaporation, surface-water 

flow into and out of Mirror Lake, and ground-water 

conditions. The monitoring effort involved installation of a weather station, three stream gages, and a 

network of piezometers and bedrock wells. In 

addition, geophysical surveys and aquifer tests were 

carried out to determine the thickness and hydraulic 

conductivity of the glacial deposits. Results of this 

study established a basic understanding of the 

geohydrologic setting and the water budget of Mirror 

Lake.

In 1989, the U.S. Geological Survey began an 

additional study to investigate ground-water flow 

through fractured crystalline rocks that underlie the 

Mirror Lake area. This study seeks to develop 

hydrologic, geologic, geophysical, and geochemical 

methods to characterize ground-water flow and 

chemical transport in fractured rocks (Shapiro and 

Hsieh, 1991). Although the Mirror Lake area is not 

contaminated, the techniques and understanding 

developed in the study are transferable to other 

fractured rock sites that have been contaminated with 

chemical wastes. One component of the fracturedrock hydrology study is to develop a computer model 

to simulate the three-dimensional, steady-state flow of 

ground water through glacial deposits and crystalline 

bedrock in the Mirror Lake area. Results of this 

modeling study are presented in this report.
INTRODUCTION

PURPOSE AND SCOPE

The purpose of the modeling study is to 

estimate, for the Mirror Lake area, (1) the rate of 

recharge from precipitation to the water table, (2) the 

hydraulic conductivity of glacial deposits and 

bedrock, (3) the three-dimensional extent of the 

Mirror Lake ground-water basin, and (4) the groundwater budget. The recharge rate is an average value 

over the entire study area. The hydraulic conductivities are effective values that characterize entire 

hydrogeologic units or large zones within units. The 

Mirror Lake ground-water basin refers to the volume 

of subsurface through which ground water flows from 

the water table to Mirror Lake, or to one of the streams 

that flow into Mirror Lake. The ground-water budget 

gives an accounting of recharge to the ground-water 

basin, discharge from the basin, and flow between 

hydrogeologic units in the basin.

The scope of the modeling study includes:

(1) conceptualizing the ground-water system,

(2) developing a computer model of the ground-water 

system, (3) calibrating the model to determine optimal 

parameters, and (4) simulating ground-water flow to 

determine the ground-water basin and ground-water 

budget. Conceptualization of the ground-water 

system is based on previous studies and hydrologic 

data that include: hydraulic heads measured in 

piezometers and wells, streamflow measured at stream 

gages, hydraulic properties determined from aquifer 

tests, and lithologic information obtained from field 

observations, well drilling, and geophysical measurements. Ground-water flow is simulated by the 

U.S. Geological Survey Modular Three-Dimensional 

Finite-Difference Ground-Water Flow Model 

(McDonald and Harbaugh, 1988), commonly known 

as MODFLOW. The model is calibrated using the 

computer program MODFLOWP (Hill, 1992), which 

employs nonlinear least-squares regression to 

determine optimal parameters. Included in the 

calibration process is an assessment of parameter 

uncertainty and goodness of model fit. The threedimensional extent of the Mirror Lake ground-water 

basin is delineated by tracking ground-water flow 

paths with the computer program MODPATH 

(Pollock, 1994), a particle tracking post-processor 

package for MODFLOW. The ground-water budget is 

calculated using the computer program 

ZONEBUDGET (Harbaugh, 1990), which calculates

subregional water budgets using simulation results 

from MODFLOW.

STUDY AREA

The study area is in central New Hampshire 

about 15 km north of the town of Plymouth (fig. 1). It

r\

occupies about 10 km, and includes Mirror Lake. 

The study area is bounded by the Pemigewasset River 

to the east, and Leeman's Brook to the north. The 

western boundary is formed by a small stream locally 

known as Paradise Brook, and a north-northeastward 

flowing segment of Hubbard Brook. The southern 

boundary is formed by the ridge south of Hubbard 

Brook. For the purpose of this report, the study area is 

referred to as the "Mirror Lake area."

Within the Mirror Lake area, the U.S. Geological Survey studies have focused on the Mirror Lake 

ground-water basin. At the outset, however, the 

extent of the ground-water basin was not well known. 

Therefore, the surface-water basin, an area of

SJ

0.85 km outlined by the topographic divide 

surrounding Mirror Lake (fig. 2), was initially chosen 

as the area of focused investigation. As of 1994, most 

piezometers and bedrock wells in the Mirror Lake area 

were drilled within the surface-water basin (fig. 3). 

However, because it was expected that the surfacewater basin might not entirely coincide with the 

ground-water basin, several piezometers and bedrock 

wells were drilled outside the surface-water basin.

PHYSICAL SETTING

The topography of the Mirror Lake area is 

characterized by steep hillsides and relatively flat 

valleys. The higher altitudes of the study area are 

dominated by ridges and bedrock knobs, which 

typically have a thin mantle of glacial deposits. 

Relatively flat kame terraces occur at several levels on 

the steep hillsides. At the lower altitudes, thicker 

glacial and alluvial deposits cover the slopes and 

valleys. Land-surface altitude ranges from 180 m 

above sea level at the Pemigewasset River to 700 m 

above sea level at the ridge top near the northwestern 

corner of the study area.

Within the Mirror Lake surface-water basin, 

land-surface altitude ranges from 213m above sea 

level at the lake shore to 470 m above sea level at the
BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

7F431 7F421

43°58'

43°57'

43°56'

Base from U.S. Geological Survey

1:24,000 Woodstock, New Hampshire quadrangle

Land surface elevations converted from feet to meters

72'

0 1 KILOMETER 

J_________

1/2 1 MILE 

CONTOUR INTERVAL 40 METERS

EXPLANATION

FS7 o Bedrock well

~^\ Stream, perennial or intermittent

Figure 1. Location, extent, and topography of study area.

43
INTRODUCTION

71°42'

43°57

Boundary of 

X Mirro/L 

. surf^6e-X«ater

34b/.

SJ

Base from U.S. Geological Survey Q

1:24,000 Woodstock, New Hampshire quadrangle L

Land surface elevations converted from feet to meters Q

EXPLANATION

0.25 0.5 KILOMETER

1/4 

CONTOUR INTERVAL 40 METERS

1/2 MILE

AREA OF LACUSTRINE MUD SEDIMENTS

AREA OF SAND AND GRAVEL DEPOSITS ALONG SOUTHERN SHORE OF MIRROR LAKE 

Road 

^ Stream, perennial or intermittent

Figure 2. Mirror Lake surface-water basin; sub-basins W, NW, and E; area of lacustrine mud 

sediments; and area of sand and gravel deposits along southern shore of Mirror Lake.
NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

43°57' -

12o 70 4583

100 90 *508 

6«

130

Base from U.S. Geological Survey

1:24,000 Woodstock, New Hampshire quadrangle „

Land surface elevations converted from feet to meters

CONTOUR INTERVAL 40 METERS

EXPLANATION

INSTRUMENTATION SITE AND IDENTIFIER 

EA Flume 

W2« Piezometer

T1 • Bedrock well and Piezometer nest 

H1O Bedrock well 

A—— D LINE OF GEOLOGIC SECTION (See fig. 4)

1/2 MILE

Figure 3. Locations of flumes, piezometers, and bedrock wells as of 1994.
GROUND-WATER HYDROLOGY

highest point of the drainage divide. The surfacewater basin contains three sub-basins, which have 

been designated as W, NW, and E (fig. 2). Each 

sub-basin contains a perennial stream that flows into 

Mirror Lake. The streams are designated by the same 

names as their respective sub-basins. Streams W and 

NW have several tributaries, some of which are 

intermittent. Physiographic characteristics of the three 

sub-basins are described by Winter (1984).

Mirror Lake is a kettle lake occupying about 

150,000 m2 (Winter, 1985). The lake has a maximum 

depth of 11 m and an average depth of about 6 m 

(Winter, 1984). A small dam at the lake's outlet on the 

southern shore maintains a near-constant lake level. 

At the dam, the lake discharges into an outlet stream, 

which flows into Hubbard Brook.

The climate of the Mirror Lake area is classified 

as humid continental characterized by short, cool 

summers and long, cold winters (Likens and others, 

1977). Mean air temperatures are 19°C in July and 

-9°C in January (Federer and others, 1990). At the 

Forest Service station about 500 m west of Mirror 

Lake (fig. 2), the average precipitation from 1978 to 

1986 was 124 cm/year (Federer and others, 1990). 

Monthly precipitation is fairly uniform throughout the 

year. Typically, 8 to 12 cm of rain or snow falls each 

month. A continuous snowpack that is as much as 

1.5 m deep accumulates each winter.

PREVIOUS AND RELATED STUDIES

As of 1996, findings from U.S. Geological 

Survey studies in the Mirror Lake area have been 

presented in reports, journal papers, and other publications. As noted above, Winter's (1984) report 

detailing studies from 1979 to 1983 established the 

basic understanding of the geohydrologic setting of 

Mirror Lake. Winter and others (1989) analyzed 

streamflow records of Streams W and NW. 

Rosenberry and Winter (1993) estimated the water 

budget of Mirror Lake. Mau and Winter (1997) 

estimated recharge and baseflow from streamflow 

hydrographs. Shattuck (1991) evaluated the effect of 

evapotranspiration on the interaction of ground water 

with Stream W. Results of borehole geophysical 

logging in the Mirror Lake area are reported by Paillet 

(1985), Hardin and others (1987), Paillet and others 

(1987), and Paillet and Kapucu (1989). Results of 

applying surface geophysical methods to detect

subsurface fractures in the Mirror Lake area are 

reported by Haeni and others (1993) and Lane and 

others (1995). Harte (1992) simulated ground-water 

flow in a vertical section of a hillside in the study area. 

Hsieh and others (1993) prepared a field guide 

describing methods of investigations employed in 

the U.S. Geological Survey studies. A collection of 

15 short papers describing preliminary results of 

fractured rock studies at the Mirror Lake area can be 

found in the conference proceedings edited by 

Morganwalp and Aronson (1996).

Research findings at the Hubbard Brook Experimental Forest are reported in numerous papers, 

reports, and books. A publication prepared by the 

U.S. Department of Agriculture Forest Service (1991) 

gives a selected listing. A complete list of titles up to 

1991 has been compiled by Likens (1991). Some of 

the research results have been synthesized in three 

books, the latest of which is edited by Likens (1985).

ACKNOWLEDGMENTS

The authors are grateful to individuals of the 

U.S. Department of Agriculture Forest Service, 

especially Wayne Martin, Chris Eager, and the late 

Robert Pierce, for use of facilities at the Hubbard 

Brook Experimental Forest. We thank the Science 

Advisory Committee of the Hubbard Brook 

Ecosystem Study for permission to carry out drilling 

and field investigations. Many individuals of the 

U.S. Geological Survey have contributed to this study. 

They include Alien Shapiro, Tom Winter, Don 

Rosenberry, Fred Paillet, Chris Barton, Phil Harte, 

Carole Johnson, and Rick Perkins. Special thanks are 

due to property owners in the Mirror Lake area who 

graciously allow access to their property. Without 

their cooperation, data collection would have been 

significantly limited.

GROUND-WATER HYDROLOGY

HYDROGEOLOGIC SETTING

The two major hydrogeologic units in the 

Mirror Lake area are glacial deposits and crystalline 

bedrock (fig. 4). In addition, a layer of organic mud 

sediments occurs on much of the bottom of Mirror
BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

Vertical exaggeration x 3.2

Figure 4. Hydrogeologic 

section A-B-C-D. (See 

figure 3 for location of 

line of section.)

500 1,000 FEET

Lake. Alluvium, which occurs along the Pemigewasset River, is not included in the present study.

GLACIAL DEPOSITS

Glacial deposits cover much of the land surface 

in the Mirror Lake area. Thickness of the glacial 

deposits ranges from 0 m where bedrock crops out to 

about 50 m at a moraine that forms the ridge north of 

Mirror Lake. Glacial deposits are believed to be 

absent from the central part of the bottom of Mirror 

Lake. In this area, lacustrine mud sediments are 

presumed to extend from lake bottom to bedrock. 

Thickness of glacial deposits, compiled from data 

provided by Winter (1984), PT. Harte (written 

communications, 1994), and C.C. Barton (written 

communications, 1994), is shown in figure 5.

The glacial deposits consist primarily of till, but 

local deposits of sand and gravel are present within 

and overlying the till. The till is highly heterogeneous, 

containing clay, silt, sand, gravel, and boulders. Fineto coarse-grained sands form kame terraces at several 

levels on the hillsides. A deposit of sand and gravel 

lies directly south of Mirror Lake, in the area between 

Mirror Lake and Hubbard Brook (fig 2). This deposit 

is believed to be a glacial delta deposit that formed 

when Hubbard Brook adjusted its grade to the Pemigewasset River during deglaciation (Winter, 1985). The 

deposit plays an important role in the hydrology of 

Mirror Lake, as water seeps from the lake into the 

highly permeable sand and gravel and flows toward 

Hubbard Brook.

T.C. Winter and D.O. Rosenberry (written 

communication, 1994) and Wilson (1991) performed 

slug tests (see Freeze and Cherry, 1979, p. 339) in 

piezometers to determine the hydraulic conductivity of 

the glacial deposits. Winter and Rosenberry tested 

16 piezometers and analyzed their test data using the 

method of Bouwer and Rice (1976). Wilson tested 

33 piezometers and analyzed the data using the 

method of Cooper and others (1967). The ranges of 

hydraulic conductivities determined from 15 piezometers in which slug tests were carried out in both studies 

are shown in table 1. Both sets of results indicate that 

the sand and gravel are about ten to one thousand 

times more permeable than the till. However, the 

ranges determined by Wilson are approximately three 

to six times the ranges determined by Winter and 

Rosenberry. The disagreement is probably due to 

different methods of analysis. In particular, the 

assumption of purely horizontal flow in Wilson's 

analysis (method of Cooper and others, 1967) might 

not accurately represent field conditions during the 

slug tests. To re-analyze some of Wilson's data,

Table 1. Ranges of hydraulic conductivities of glacial deposits as determined by slug tests. 

[Hydraulic conductivity in meters per second]

Type of glacial deposit

Range of hydraulic conductivity determined by

Winter and Rosenberry 

(written communication, 1994)

Wilson 

(1991)

Till

Sand and gravel

2.1xlO-7 tol.3xl(r6 

4.6 x 1(T6 to 2.9 x lO'4

5.2 x 10'7 to 7.3 x 1(T6 

2.7 x 1(T5 to 1.3 x 10'3
GROUND-WATER HYDROLOGY

43°57

71°42' 

I

Base from U.S. Geological Survey

1:24,000 Woodstock, New Hampshire quadrangle I 

Land surface elevations converted from feet to meters A

0.25 0.5 KILOMETER

1/4 1/2 MILE

CONTOUR INTERVAL 40 METERS

EXPLANATION

-15— LINE OF EQUAL THICKNESS OF GLACIAL DEPOSITS-Interval is 5 meters. Dashed where inferred. 

B AREA OF LACUSTRINE MUD SEDIMENTS

^ Bedrock outcrop or roadcut

• Piezometer

• Bedrock well and Piezometer nest 

O Bedrock well

Figure 5. Thickness of glacial deposits, area of lacustrine mud sediments, and areas of bedrock outcrop 

and roadcut in the vicinity of Mirror Lake.
10 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

Demir (1992) used a numerical model that considers 

both horizontal and vertical flow around the piezometer. He obtained hydraulic conductivity values that 

are 0.4 to 0.5 times those of Wilson. The lower values 

are in closer agreement with the results of Winter and 

Rosenberry, and are probably more representative than 

are Wilson's results of the glacial deposits in the 

Mirror Lake area.

BEDROCK

The bedrock underlying the Mirror Lake area is 

composed of four major rock types: schist, granite, 

pegmatite, and lamprophyre. The schist is part of 

the Rangeley formation of early Silurian age (438 to 

428 million years old) (Lyons and others, 1986). It is 

intruded extensively by the Concord granite of late 

Devonian age (370 to 365 million years old). Both the 

schist and granite are intruded by pegmatite dikes, 

which are possibly a residual differentiate of the 

Concord granite. All three rocks are intruded by lesser 

amounts of lamprophyre, a fine-grained volcanic dike 

rock. The age of the lamprophyre is not specifically 

known but is likely between middle Jurassic and early 

Cretaceous (190 to 95 million years old), based on 

ages of similar dikes in the area (McHone, 1984).

The spatial distribution of rock types is 

highly complex. Exposed rocks on roadcuts along 

Interstate 93 (fig. 5) show that the bedrock composition is highly variable over distances of tens of meters. 

The granitic intrusions have the form of dikes, 

irregular pods, and anastomosing fingers, ranging 

from centimeters to tens of meters thick. The 

pegmatite and lamprophyre dikes are centimeters to 

meters thick. Drill cuttings and video image logs from 

bedrock wells show that there is generally poor 

correlation of rock types between adjacent wells 

drilled tens of meters apart. Because of these 

complexities, the distribution of rock type was not 

determined in this study. Instead, the entire bedrock 

complex is considered as one hydrogeologic unit.

Fractures are the major conduits through which 

ground water flows through the bedrock, as the intact 

rock has extremely low hydraulic conductivity. 

Barton (see Hsieh and others, 1993) mapped and 

analyzed fractures on road cuts along Interstate 93 and 

found a dominant strike in the N 30° E direction. 

However, fracture dips are highly variable. The 

granite appears more fractured than the schist. More 

than half of the exposed fractures have blind terminations, ending in the rock matrix without intersecting 

other fractures. This feature suggests that the fracture 

network is poorly connected.

Hydraulic conductivity of the fractured bedrock 

has been determined by single-well and multiple-well 

hydraulic tests. In a single-well hydraulic test (Hsieh 

and others, 1993), two packers are used to isolate a 

4- to 5-m interval in a well. Hydraulic conductivity of 

the interval is determined by injecting or withdrawing 

water for approximately 20 minutes while measuring 

the flow rate and change in hydraulic head. Because 

of the relatively short test duration, these tests characterize a small volume of the bedrock within several 

meters of the well. Results of 208 single-well tests in 

the Mirror Lake area show that hydraulic conductivity 

of the fractured bedrock varies over at least 6 orders of 

magnitude, from less than 1 x 10" 10 m/s (lower limit of 

measurement of test equipment) to 5 x 10 m/s. A 

histogram of the test results is shown in figure 6. The 

large range of values indicates that, at a scale of

50

Below -10-9 -8 -7 

Detection 

Limit

-6 -5

LOG HYDRAULIC CONDUCTIVITY,

IN METERS PER SECOND 

Figure 6. Histogram of hydraulic conductivities of 

bedrock determined by single-well hydraulic tests.
GROUND-WATER HYDROLOGY 11

several meters, the hydraulic conductivity of the 

fractured bedrock is extremely heterogeneous.

In a multiple-well hydraulic test, water is 

withdrawn from a packer-isolated interval in one well, 

while changes in hydraulic head are measured in 

packer-isolated intervals in adjacent wells. Multiplewell hydraulic tests have been conducted at the 

FSE well field (fig. 3), where 13 wells were drilled in a 

100-m by 100-m area. Test results indicate that the 

bedrock underlying the well field contains 4 highly 

conductive fracture clusters, each cluster occupying 

a near-horizontal, tabular volume approximately 

1.5 m thick by 20 to 50 m in horizontal extent, and 

having a hydraulic conductivity of 6 x 10"5 m/s 

(Hsieh and Shapiro, 1994). The highly conductive 

fracture clusters are hydraulically connected to each 

other via a less conductive fracture network haying 

horizontal and vertical hydraulic conductivities of 3 x 

10" 8 and 2 x 10"7 m/s, respectively. These findings 

suggest that, over distances of tens of meters, ground 

water could travel through the highly conductive 

fracture clusters. However, over distances of 100 m or 

more, ground water must flow through the less 

conductive fracture network. Consequently, on a scale 

of 100 m, the effective hydraulic conductivity of the 

bedrock underlying the FSE well field is strongly 

controlled by the less conductive fracture network. 

Hsieh and Shapiro (1994) estimated this effective 

hydraulic conductivity to be 2 x 10"7 m/s.

LAKE SEDIMENTS

Lacustrine mud sediments lie over much of the 

bottom of Mirror Lake. The extent of these sediments 

was mapped by Moeller (1978) and is shown in 

figure 2. The maximum thickness of the lacustrine 

mud is approximately 12m (Davis and Ford, 1982). 

In the central part of the lake bottom, the lacustrine 

mud is believed to lie directly over bedrock 

(Rosenberry and Winter, 1993); away from the central 

region, the lacustrine mud overlies glacial deposits 

(fig. 5). Along the near-shore lake bottom (littoral 

zone), the lacustrine mud is absent. Hydraulic 

conductivity of the lacustrine mud is significantly less 

than that of the glacial deposits. On the basis of 

hydraulic properties of similar lake sediments reported 

in the literature (Boelter, 1972; Frape and Patterson, 

1981), Rosenberry and Winter (1993) estimated the 

hydraulic conductivity of the Mirror Lake lacustrine 

mud to be between 10" 10 to 10"9 m/s.

GROUND-WATER FLOW

Ground-water flow in the Mirror Lake area is 

typical of mountain-and-valley terrain of New 

England uplands. On hillsides, water from precipitation and snowmelt infiltrates to the water table, flows 

downslope through the saturated glacial deposits and 

fractured bedrock, and discharges to streams and to 

Mirror Lake. Along a part of the south side of Mirror 

Lake, lake water seeps into the subsurface and flows 

toward Hubbard Brook. Part of this subsurface flow 

discharges into a wetland that drains into the outlet 

stream of Mirror Lake.

RECHARGE

Recharge is the flow of water into the groundwater system, that is, the saturated zone below the 

water table or surface-water bodies. The three sources 

of recharge in the Mirror Lake area are: (1) infiltration 

from precipitation (including snowmelt) to the water 

table, (2) subsurface seepage from Mirror Lake, and 

(3) infiltration of stream water along losing reaches of 

streams. Infiltration from precipitation is the largest 

source of recharge. By contrast, infiltration of stream 

water is relatively minor and occurs only along short, 

local reaches of streams.

Although monthly precipitation in the Mirror 

Lake area is fairly uniform throughout the year, 

recharge from precipitation (infiltration that reaches 

the water table) varies seasonally. In winter, recharge 

is lower because infiltration from snowmelt at the base 

of the snowpack is minimal. In spring, recharge is 

greater due to rapid snowmelt and minimal 

evapotranspiration. In summer, recharge is lower 

because much of the precipitation is lost to the 

atmosphere by evapotranspiration from the soil zone. 

In late autumn, decreased evapotranspiration leads 

again to increased recharge. The amount of recharge 

from precipitation and snowmelt in the Mirror Lake 

area was estimated by Mau and Winter (1997) from 

streamflow recession analysis. Their results indicate 

that recharge from precipitation and snowmelt varies 

from 20 to 56 cm/year, or about 16 to 45 percent of 

precipitation.

Water seeps from Mirror Lake into the subsurface along a part of the south side of the lake. In this 

area, a deposit of highly permeable sand and gravel 

allows rapid infiltration of lake water into the groundwater system (fig. 2). On the basis of hydraulic heads
12 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

observed in piezometers and wells, hydraulic 

conductivities determined from slug tests, and glacialdeposit thickness inferred from seismic refraction 

surveys, Rosenberry and Winter (1993) estimated 

the seepage from Mirror Lake to be approximately 

260,000 m3/year. However, they noted that this 

value is subject to large uncertainty, approaching 

±100 percent.

Although streams in the Mirror Lake area 

generally act as ground-water drains (that is, they 

receive ground-water discharge), stream water 

infiltrates into the subsurface along short, local 

reaches. An example of such a location is where 

Stream W flows near piezometer site W3 (fig. 3), at 

the edge of a terrace. Shattuck (1991) examined the 

ground-water flow pattern in this area and found that 

where the stream flows on till, ground water enters the 

stream along both its sides and its bottom. However, 

at a bend where the stream flows on sand, ground 

water enters the stream along the uphill bank, but 

stream water infiltrates into the ground along the 

downhill bank. Although the amount of stream-water 

infiltration is minor compared with other sources of 

recharge, its presence illustrates the complexity that 

can exist in a ground-water flow pattern on a hillside.

DISCHARGE

Ground water in the Mirror Lake area 

discharges to streams and to Mirror Lake. Along a 

stream bank, a portion of the discharge is consumed by 

riparian vegetation, and the remainder enters the 

stream to become baseflow. For Streams W, NW, and 

E, baseflow at the inlet to Mirror Lake represents 

ground-water discharge to the stream minus the sum 

of riparian evapotranspiration and any infiltration of 

stream water along losing reaches.

Baseflow in Streams W and NW at the inlet to 

Mirror Lake has been estimated by Mau and Winter 

(1997) using methods developed by the Institute of 

Hydrology (1980a, b) and Rutledge and Daniel 

(1994). They concluded that baseflow accounts for 

about 25 to 57 percent of streamflow. In the present 

study, baseflow is assumed to be 40 percent of 

streamflow over the long term. Long-term average 

streamflow (Rosenberry and Winter, 1993) and 

estimated long-term average baseflow in Streams W, 

NW, and E at their inlets to Mirror Lake are shown in 

table 2.

The amount of riparian evapotranspiration in the

Mirror Lake area is not well known. Shattuck (1991) 

showed that in the summer months, riparian 

evapotranspiration causes diurnal variations in water 

levels in shallow water table wells near Stream W, but 

did not quantify the evapotranspiration rate. Rutledge 

(1993) analyzed a ten-year record of streamflow in 

166 surface-water basins in the Appalachian 

Mountains and Piedmont area of the eastern United 

States, and found that for each basin, approximately 

90 percent of ground-water recharge becomes 

baseflow in streams, while the remaining 10 percent of 

recharge is consumed by riparian (near-stream) 

evapotranspiration.

In addition to streams, Mirror Lake is also a sink 

for ground-water discharge. Most ground water seeps 

into Mirror Lake from glacial deposits along the west, 

north, east, and part of the south sides of the lake 

bottom, where the lacustrine mud sediments are absent 

(fig. 5). Minute amounts of ground water enter Mirror 

Lake by upward seepage through the lacustrine mud, 

or directly from bedrock where it crops out on the east 

side of the lake (fig. 5). Rosenberry and Winter (1993, 

fig. 3) estimated that ground-water seepage to Mirror 

Lake is approximately 80,000 m3/year. This value, 

like the estimate of subsurface seepage from Mirror 

Lake, is also subject to large uncertainty, approaching 

±100 percent.

Table 2. Long-term average streamflow and baseflow in Streams 

W, NW, and E at their inlets to Mirror Lake.

[Streamflow and baseflow in cubic meters per year]

Stream

W

NW

E

Total

Long-term average 

streamflow

157,000

249,000

12,000

418,000

Long-term 

average baseflow1

63,000

100,000

5,000

168,000

Estimated as 40 percent of streamflow.

HYDRAULIC HEAD

DATA COLLECTION

Since 1979, hydraulic heads have been 

measured weekly in piezometers and bedrock wells in 

the Mirror Lake area. Piezometer and well construction are described by Winter (1984). From 1979 

through 1993, about 70 piezometers having either 0.6 

or 1.0 m long screens were installed in the glacial
GROUND-WATER HYDROLOGY 13

deposits. The locations of 48 piezometers used for 

model calibration are shown in figure 3. Approximately half of these piezometers (those with names 

starting with W) are screened at or just below the 

water table. The remainder are grouped into nests that 

are installed at various depths adjacent to bedrock 

wells. The names of these piezometers consist of the 

bedrock well name and a suffix equal to the depth of 

the screen midpoint, in feet below land surface. 

Hydraulic head in a piezometer is assumed to 

represent hydraulic head at the midpoint of the 

piezometer screen. Altitudes of the piezometer screen 

midpoints are given in table 3. The depths of the 

screen midpoints range from 1 m to 40 m below land 

surface.

From 1979 through 1993, 30 bedrock wells 

were installed in the Mirror Lake area. Depths of 

wells range from 46 m to 305 m below land surface. 

Of the 30 wells, 13 are clustered in a 100 m by 

100 m area known as the FSE well field, 4 are 

clustered in another 100 m by 100 m area known as 

the CO well field, and the remainder are distributed 

over a 1 km by 1 km area in the vicinity of Mirror 

Lake (fig. 3). Prior to 1989, hydraulic heads were 

measured in open bedrock wells. Since 1989, inflatable packers have been installed in many of the 

bedrock wells to allow multi-level monitoring (Hsieh 

and others, 1996). Packers are positioned in a well to 

hydraulically isolate high transmissivity fractures 

from each other. Hydraulic head in a packer-isolated 

interval is assumed to represent hydraulic head in the

Table 3. Altitude of screen midpoint, long-term average hydraulic head, and averaging method for piezometers. 

[Altitude and hydraulic head in meters above sea level]

Piezometer

CO1-11

CO1-18

FS1-17

FS1-25

FS1-35

FS3-11

FS3-22

FS3-29

FSE2-23

FSE2-43

FSE5-36

FSE6-41

FSE6-59

IS1-18

IS1-38

Kl-8

Kl-39

K2-21

K2-41

K3-22

K3-61

Rl-36

Rl-55

TR1-63

Altitude 

of screen 

midpoint

214.3

211.8

255.7

253.9

250.9

271.0

267.6

265.5

233.5

226.8

231.3

229.8

224.3

220.0

213.1

211.8

202.3

208.4

202.3

200.5

188.6

246.0

239.9

229.8

Long-term 

average 

hydraulic 

head

215.71

215.52

259.62

259.90

260.81

271.53

271.13

270.48

234.40

234.14

237.30

237.57

237.36

221.70

221.81

212.78

211.82

210.40

210.48

207.08

205.78

249.15

250.71

233.90

Averaging 

method

1991 adjusted

1991 adjusted

Direct

Direct

Direct

Direct

Direct

1991 adjusted

Direct

Direct

Dec 93 adjusted

1991 adjusted

1991 adjusted

Dec 93 adjusted

Dec 93 adjusted

Direct

Direct

Direct

Direct

Direct

Direct

1993 adjusted

1993 adjusted

Direct

Piezometer

TR1-132

TR2-23

W2

W3

W3A

W5

W6

W7

W8

W9

W10

Wll

W13

W14

W15

W16

W18

W19

W20

W21

W24

W25

W26

W27

Altitude 

of screen 

midpoint

208.2

225.9

251.9

251.5

256.0

209.6

228.4

212.3

212.3

211.8

209.6

232.4

209.5

207.1

261.0

257.9

227.1

214.0

213.6

214.2

208.5

254.5

251.5

262.4

Long-term 

average 

hydraulic 

head

233.02

232.09

252.08

257.84

257.10

210.69

231.70

212.94

212.77

212.45

211.22

234.11

210.09

207.25

261.86

262.33

227.78

214.19

213.77

213.95

208.52

256.70

252.72

263.23

Averaging 

method

1991 adjusted

1993 adjusted

Direct

Direct

Direct

Direct

Direct

Direct

Direct

Direct

Direct

Direct

Direct

Direct

Direct

Direct

Direct

Direct

Direct

Direct

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted
14 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

most transmissive fracture of the interval. This 

assumption is reasonable because the most transmissive fracture in an interval is typically at least one 

order of magnitude more transmissive then the other 

fractures in the interval. Consequently, hydraulic head 

in the most transmissive fracture dominates hydraulic 

head in the interval. Locations of packer-isolated 

intervals in each well, and the positions of the most 

transmissive fracture in each interval are shown in 

table 4. In each well, packer-isolated intervals are 

named in alphabetical order, beginning with the 

uppermost interval. For example, the three intervals 

of well FS1 are denoted, from top to bottom, as FSl-a,

FS 1-b, and FS 1-c. An interval name without a trailing 

letter denotes a well without packers. In this case, the 

interval is the entire open-hole section of the well.

In September 1995, an additional bedrock well, 

FS7, was drilled on the ridge near the northwest corner 

of the study area (fig. 1). This well provides a 

hydraulic head measurement near the highest altitudes 

of the Mirror Lake area, where hydraulic head data 

were heretofore unavailable due to limited drilling 

access. Because the well was drilled near the end of 

this modeling study, hydraulic head measured in the 

well was not used in calibration runs completed prior

Table 4. Altitude of top, bottom, and most transmissive fracture; long-term average hydraulic head; and averaging method for well 

intervals .

[Altitude and hydraulic head in meters above sea level]

Well interval

COl-a

CO 1-b

CO2-a

CO2-b

CO3-a

CO3-b

FSEl-a

FSEl-b

FSEl-c

FSE2-a

FSE2-b

FSE2-C

FSE3-a

FSE3-b

FSE3-C

FSE4-a

FSE4-b

FSE4-C

FSE4-d

FSE5

FSE6-a and b1

FSE7

FSE8-a

FSE8-b

FSE9-a

FSE9-b

FSElO-a

Altitude at top 

of interval

208

165

209

183

208

189

223

207

173

224

208

174

223

208

174

223

208

174

134

225

220

224

221

190

223

197

222

Altitude at 

bottom of 

interval

166

41

184

155

189

175

208

174

131

208

174

132

208

174

132

209

175

135

11

181

166

166

190

166

198

166

212

Altitude of 

most 

transmissive 

fracture

184

88

209

166

205

185

222

194

140

222

202

172

222

192

143

215

197

140

107

207

210

222

211

167

200

177

216

Long-term 

average 

hydraulic 

head

216.29

215.07

215.38

216.66

214.16

215.75

233.97

233.83

235.43

233.99

233.81

235.61

233.95

233.79

236.99

233.83

233.83

237.13

240.42

233.84

235.52

235.00

234.87

234.97

234.99

234.98

235.65

Averaging method

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

Dec 93 adjusted

Dec 93 adjusted

Dec 93 adjusted

Dec 93 adjusted

Dec 93 adjusted

Dec 93 adjusted
GROUND-WATER HYDROLOGY 15

Table 4. Altitude of top, bottom, and most transmissive fracture; long-term average hydraulic head; and averaging method for well 

intervals—Continued.

[Altitude and hydraulic head in meters above sea level]

Well interval

FSElO-b

FSEll-a

FSEll-b

FSE12-a

FSE12-b

FSE13-a 

FSE13-b

FSl-a 

FSl-b

FSl-c

FS2-a

FS2-b

FS2-c

FS3-a

FS3-b

FS3-c

FS4

FS7

Hl-a

Hl-b

Hl-c

ISl-a

ISl-b

ISl-c

Kl

K2

K3

Rl-a

Rl-b

Rl-c

Tl-a

Tl-b

Tl-c

TRl-a

TRl-b

TRl-c

TR2-a

TR2-c

Altitude at top 

of interval

211

218

183

223

178

216 

198

245 

229

195

246

223

189

262

247

214

343

660 (approx.)

223

192

159

208

202

149

200

198

186

237

225

192

223

196

169

197

187

157

217

165

Altitude at 

bottom of 

interval

168

184

158

179

156

199 

164

229 

196

124

223

190

103

248

215

48

45

610 (approx.)

193

159

140

203

150

72

168

166

154

226

192

63

197

169

76

188

157

58

199

1

Altitude of 

most 

transmissive 

fracture

208

197

163

187

158

203 

193

235

223

166

226

212

135

253

244

162

not known

not known

198

176

155

207

151

101

not known

not known

not known

234

221

109

208

179

161

194

163

101

201

148

Long-term 

average 

hydraulic 

head

235.32

234.98

235.00

235.34

236.40

234.9 1 1

256.151

256.50

247.32

247.59

240.92

269.13

268.43

265.05

350.54

650 (approx.)

220.86

218.50

220.24

221.60

219.28

217.63

212.65

211.04

208.34

250.06

250.06

251.39

224.95

224.79

224.90

231.10

228.96

225.72

233.02

233.50

Averaging method

Dec 93 adjusted

Dec 93 adjusted

Dec 93 adjusted

Dec 93 adjusted

Dec 93 adjusted

Dec 93 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

Dec 93 adjusted

single measurement, 

Sept 95

1993 adjusted

1993 adjusted

1993 adjusted

1993 adjusted

1993 adjusted

1993 adjusted

Direct

Direct

Direct

1993 adjusted

1993 adjusted

1993 adjusted

1993 adjusted

1993 adjusted

1993 adjusted

1991 adjusted

1991 adjusted

1991 adjusted

1993 adjusted

1993 adjusted

^wo intervals combined into one due to failure of packer to hydraulically isolate the two intervals.
16 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

to September 1995. Instead, an additional calibration 

run was carried out to include the new data.

TEMPORAL VARIATION OF HYDRAULIC HEAD

Hydraulic head in glacial deposits and bedrock 

of the Mirror Lake area exhibits both seasonal and

shorter-period variations. To illustrate these 

variations, hydrographs from 1991 through 1993 in 

piezometer FS1-17 and in packer-isolated intervals of 

bedrock well FS3 are shown in figure 7. Seasonal 

variation in these hydrographs is caused by the 

seasonal variation in recharge from precipitation. In 

spring, hydraulic head typically rises to an annual

261.0

258.5

1991 1992 1993

B

>3-b and FS3-C 

to July 1993

264.0

1991 1992 1993

Figure 7. Hydraulic head from 1991 through 1993 in piezometer FS1-17 (A) and in 

packer-isolated intervals of bedrock well FS3 (B).
MODEL OF STEADY-STATE GROUND-WATER FLOW 17

maximum due to high recharge from snowmelt and 

rainfall. In summer, hydraulic head declines because 

of lower recharge caused by increased evapotranspiration. In autumn, hydraulic head usually rises again, 

as evapotranspiration decreases, and typically reaches 

a peak that is smaller than the spring peak. In winter, 

hydraulic head generally declines due to lower 

recharge, but might also rise if mild temperatures in 

midwinter melt some of the snowpack. Superimposed on the seasonal variation in hydraulic head are 

shorter-period variations that are due to individual 

periods of precipitation and snowmelt.

LONG-TERM AVERAGE HYDRAULIC HEAD

The steady-state flow model developed in this 

study is designed to simulate long-term average flow 

conditions. Therefore, the model should be calibrated 

to long-term average hydraulic heads. Because data 

have been collected in piezometers and wells over 

varying durations, two different methods are used to 

determine long-term averages. The results are given 

in tables 3 and 4.

For piezometers and wells that have five or 

more years of hydraulic-head measurements, a 

"direct" method is used to determine long-term 

averages. This method involves first computing the 

annual average hydraulic head for all calendar years 

in which there are at least 26 weekly measurements. 

Then the long-term average hydraulic head is 

computed as the average of the annual averages. This 

method is denoted as "Direct" in tables 3 and 4.

For piezometers and wells that have less than 

5 years of hydraulic-head measurements, an "adjustment" method is used to estimate long-term averages. 

In many cases, data are sufficient to determine the 

annual average hydraulic head only for 1991. 

Because this annual average might not be representative of the long-term average, an adjustment is 

applied. The adjustment is based on a best-fit, linear 

relation between the 1991 annual average hydraulic 

head and the long-term average hydraulic head for 

wells and piezometers for which both averages can be 

calculated. This linear relation is then used to 

estimate the long-term average hydraulic heads in 

cases where data are sufficient to determine only the 

1991 annual average hydraulic head. This method is 

denoted by "1991 adjusted" in tables 3 and 4. A 

similar method is applied to a few cases for which 

data are sufficient to determine only an annual

average hydraulic head for 1993, or a monthly average 

hydraulic head for December 1993. The latter 

methods are denoted respectively by "1993 adjusted" 

and "Dec 93 adjusted" in tables 3 and 4.

The long-term average altitude of the water 

table as determined from the long-term average 

hydraulic head data is shown in figure 8. The water 

table is generally within 20 m of land surface, and in 

many locations, less than 10m below land surface. 

Hydraulic head measured at the newly drilled well 

FS7 indicates that the water table may remain shallow 

all the way to the highest altitude at the northwest 

corner of the study area. Near this ridge top, hydraulic 

head in well FS7 is approximately 15m below land 

surface.

The distribution of long-term average hydraulic 

heads in the vertical section A-B through four wells on 

a hillside west of Mirror Lake (see fig. 8) is shown in 

figure 9. Although the hydraulic gradient indicates a 

predominance of lateral flow in the downhill direction, 

vertical components of flow are also apparent. 

Ground water flows downward from glacial deposits 

to bedrock along most of the section from well FS3 to 

well FS2. Near FSE4, however, upward flow is significant in the bedrock.

MODEL OF STEADY-STATE 

GROUND-WATER FLOW

Ground-water flow in the glacial deposits and 

crystalline bedrock of the Mirror Lake area is 

simulated by the U.S. Geological Survey Modular 

Three-Dimensional Finite-Difference Ground-Water 

Flow Model (McDonald and Harbaugh, 1988), 

commonly known as MODFLOW. The steady-state 

model in this study simulates long-term average flow 

conditions. Seasonal and shorter-period variations are 

not considered. Because MODFLOW is a continuum 

model, its application to fractured rock presupposes 

that the rock can be treated as an equivalent porous 

medium.

ASSUMPTIONS

The model developed in this study is based on 

four simplifying assumptions. These assumptions are: 

(1) recharge from precipitation is areally uniform,
18 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

43°57'

71°42' 

I

Base from U.S. Geological Survey 0 

1:24,000 Woodstock, New Hampshire quadrangle [•

0.25 

I

0.5 KILOMETER 

I__________ * ~ i i i

Land surface elevations converted from feet to meters 0 1/4 1/2 MILE

TOPOGRAPHIC CONTOUR INTERVAL 40 METERS

EXPLANATION

211-2 • PIEZOMETER

237.4 • BEDROCK WELL AND PIEZOMETER NEST

350.5 o BEDROCK WELL FS4

Number adjacent to symbol indicates long-term average altitude of water table, in meters above sea level.

—230— WATER TABLE CONTOUR-Shows long-term average altitude of water table.

Contour interval 10 meters. Datum is sea level

A——B LINE OF HYDROLOGIC SECTION (See fig. 9) 

Figure 8. Long-term average altitude of water table in the vicinity of Mirror Lake.
MODEL OF STEADY-STATE GROUND-WATER FLOW 19

FSE4

0 300 FEET

No Vertical Exaggeration

EXPLANATION

——240— HYDRAULIC HEAD CONTOUR-Shows long-term 

average hydraulic head. Contour interval 5 

meters. Datum is sea level.

•a HYDRAULIC HEAD MEASUREMENT LOCATION-Denotes well

interval and location of the highest-transmissivity fracture in the 

well interval.

Figure 9. Long-term average hydraulic heads in section A-B. (See figure 8 for location of line of 

section.)

(2) riparian evapotranspiration is a small component 

of ground-water discharge to streams and can be 

neglected, (3) flow in the bedrock occurs primarily in 

the uppermost 150 m, and (4) the distribution of 

hydraulic conductivity can be represented by dividing 

the model domain into several zones each having 

uniform properties.

The assumption of spatially uniform recharge 

from precipitation is adopted because of the lack of 

more detailed information. Although variations in 

precipitation, slope, vegetation, and soil could lead to 

variations in infiltration, studies of such spatial 

variability have not been carried out in the Mirror 

Lake area. To the extent that recharge from precipitation is assumed uniform in the model, the simulation 

results might not accurately represent local details in 

the ground-water system. However, if the recharge 

variability is not great, then the simulation results 

should provide a useful characterization of the overall 

ground-water system and ground-water budget.

The assumption of negligible riparian 

evapotranspiration is based on the work of Rutledge 

(1993), who estimated that for several surface-water 

basins in the Appalachians and Piedmont, approximately 10 percent of recharge to ground water is 

consumed by riparian evapotranspiration, and the 

remaining 90 percent of recharge becomes baseflow in 

streams. In the model of ground-water flow in the 

Mirror Lake area, riparian evapotranspiration is 

neglected. Therefore, at the inlet to Mirror Lake, it is 

assumed that baseflow in the W, NW, or E Stream 

equals ground-water discharge to the stream less any 

infiltration of stream water along losing reaches.

The assumption that flow in the bedrock occurs 

primarily in the uppermost 150 m is based on results 

of hydraulic tests using packers to determine the 

hydraulic conductivity distribution in wells to depths 

of about 250 m. The test results indicate that most of 

the permeable fractures are within 150 m of the 

bedrock surface. Although permeable fractures are
20 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

present at greater depths, they are relatively 

uncommon. By limiting ground-water flow in the 

bedrock to the uppermost 150 m, the model in this 

study ignores the possibility of deeper ground-water 

flow at larger regional scales that extend beyond the 

study area.

The distribution of hydraulic conductivity in the 

model is represented by a "zonation" approach. In this 

approach, the model domain is divided into several 

zones. Within each zone, hydraulic conductivity is 

uniform but can be isotropic or anisotropic. The zones 

and the associated hydraulic conductivities define the 

"parameter structure." The parameter structure 

characterizes the large-scale variation of hydraulic 

conductivity in the model domain. Variations on 

scales smaller than the zone dimension are neglected. 

Consequently, the simulation results characterize the 

overall ground-water system, but not local details.

GOVERNING EQUATION

In the absence of internal sources and sinks, the 

equation of three-dimensional, steady-state groundwater flow is (Freeze and Cherry, 1979, p. 64):

5~ AK /-/5~ + + 5~ KK H^~ + + 5~ AK ~ - - nU ' 

ax\ ndxJ dy\ tldyJ dz

where

h = hydraulic head (L), 

Kfj = horizontal hydraulic conductivity (L/T), 

Ky = vertical hydraulic conductivity (L/T), 

x, y = cartesian coordinates in the horizontal

directions (L),

z = cartesian coordinate in the vertical direction (L). 

Implicit in equation 1 is the assumption that two 

principal directions of hydraulic conductivity are 

horizontal and the third is vertical. Hydraulic conductivity is isotropic in the horizontal plane, but could 

exhibit horizontal-to-vertical anisotropy. In the 

model, hydraulic conductivities of the different zones 

are designated by additional subscripts.

Equation 1 is solved using MODFLOW with the 

preconditioned conjugate gradient method (PCG2) 

(Hill, 1990). The MODFLOW computer code is 

modified slightly so that the drying of finite-difference 

cells is delayed until the iterative solution loop is 

exited. During the solution loop, if the saturated

thickness of a cell becomes zero, the horizontal 

hydraulic conductivity is set to zero, but the vertical 

hydraulic conductivity is not changed. This modification allows cells to become resaturated from below. 

When the solution loop is exited, cells with saturated 

thickness equal to zero are removed. If the PCG2 

solver has not converged, additional passes through 

the iterative solution loop, with delayed drying of 

cells, are performed until convergence is achieved.

SPATIAL DISCRETIZATION

The three-dimensional model domain extends 

horizontally over the Mirror Lake area and vertically 

from the water table to 150 m below the bedrock 

surface. In the horizontal dimensions, the model 

domain is discretized into a grid of 89 rows by 

85 columns of rectangular cells (fig. 10). The central 

part of the grid has square cells 17 m on each side. 

Outside the central part, cell dimension increases 

towards the perimeter in two stages, from 17 m to 

50 m, and then from 50 m to 150 m. The finer, central 

part of the grid allows more accurate model simulation 

in the area of focused investigation near Mirror Lake. 

The coarser, outer part of the grid allows coverage of 

the entire study area using a reasonable number of 

cells.

In the vertical dimension, the model domain is 

discretized into five model layers. The upper two 

model layers (layers 1 and 2) represent glacial 

deposits, Mirror Lake, the lacustrine mud sediments, 

and the Pemigewasset River. The lower three model 

layers (layers 3,4 and 5) represent the bedrock. 

Conceptually, the five model layers slope with the land 

surface as illustrated by the vertical section along 

row 25 of the model grid (fig. 11). In the actual 

MODFLOW simulation, the layer slope is ignored, 

and thus an error is introduced. However, Harte 

(1994) investigated the magnitude of this error and 

found it to be relatively small.

In the central part of the grid where the 

thickness of glacial deposits is known (fig. 5), layers 1 

and 2 are each assigned half the thickness of the 

glacial deposits. In the outer part of the grid where the 

thickness of glacial deposits is unknown, layers 1 and 

2 are each assigned a uniform thickness of 4.5 m. 

Where bedrock crops out, layers 1 and 2 are absent. 

Mirror Lake is represented by cells in layer 1, and the 

lacustrine mud sediments are represented by cells in
MODEL OF STEADY-STATE GROUND-WATER FLOW 21

layer 2. As shown in figure 12, the lacustrine mud is 

assumed to lie directly on bedrock. This approximation ignores the fringe region where these sediments 

overlie glacial deposits. The Pemigewasset River is 

represented by the eastern-most cell of each model 

row in layer 1 (figure 10). The bedrock is represented 

by the lower three model layers (layers 3,4, and 5). 

The top of layer 3 is the bedrock surface, and the 

bottom of layer 5 is 150 m below the bedrock surface. 

Each of layers 3, 4, and 5 has a uniform thickness that 

is respectively 30 m, 60 m, and 60 m.

BOUNDARY CONDITIONS

The vertical surfaces that bound the sides of the 

model domain are specified as no-flow boundaries 

(fig. 10). On the east, north, and west sides, these 

vertical boundary surfaces approximately lie beneath 

the Pemigewasset River, Leeman's Brook, Paradise 

Brook, and a reach of Hubbard Brook. The no-flow 

boundary condition beneath these water bodies is a 

result of the assumption that ground water does not 

flow across the river and streams, but discharges into 

them. On the south side and near the northwest corner 

of the model domain, the vertical boundary surfaces 

lie beneath major topographic divides. The no-flow 

boundary condition in these areas is a result of the 

assumption that the topographic divides are also 

ground-water divides. The bottom surface of the 

model domain (150 m below bedrock surface) is also 

specified as a no-flow boundary. The use of this 

boundary condition ignores ground water flow in the 

bedrock at depths greater than 150 m below the 

bedrock surface.

Mirror Lake and the Pemigewasset River are 

both represented in layer 1 by cells that are assigned 

known hydraulic heads. Cells representing Mirror 

Lake are assigned a hydraulic head of 213 m above sea 

level, which is the average lake surface altitude from 

1982 through 1990. Cells representing the Pemigewasset River are assigned a hydraulic head of 177 m 

above sea level, which is the river stage obtained from 

the U.S. Geological Survey Woodstock topographic 

map (7-1/2 minute series). The small gradient in the 

river stage over the study area is neglected. The use of 

constant-head cells to represent Mirror Lake and the 

Pemigewasset River precludes determination of a 

water balance for these surface water bodies.

The position of the water table is not specified 

but is computed in the model simulation. At each 

model cell in the horizontal dimension, the computed 

water table may lie in any of the five model layers. 

During the simulation, model cells that lie above the 

water table become "inactive" and are eliminated from 

the simulation. Model cells that contain or lie below 

the water table remain "active" in the simulation.

Recharge from precipitation is represented as 

an areally uniform flux applied to the uppermost active 

cell throughout the model domain, except for cells 

representing Mirror Lake. At these lake cells, precipitation enters the lake, and there is no direct recharge to 

the ground-water system. The rate of recharge from 

precipitation is estimated during model calibration.

Streams in the model domain include Hubbard 

Brook, Norris Brook, Paradise Brook, Leeman's 

Brook, Streams W, NW and E, and their tributaries. In 

the context of the present ground-water flow model, 

streams are defined as locations where the topography 

forms a stream channel, regardless of whether or not 

flow is perennial in the channel. Therefore, certain 

stream segments in the ground-water model might not 

qualify as streams on a topographic map. The amount 

of ground-water discharge into a stream is determined 

by model simulation.

Streams are simulated in the model by the 

Streamflow Routing Package of MODFLOW, which 

simulates the interaction between streams and ground 

water, and accounts for the amount of flow in streams 

(Prudic, 1989). In this approach, streams are considered to lie over model cells in layer 1. Each stream is 

divided into many reaches, each reach corresponding 

to an individual model cell over which the stream 

traverses (fig. 10). For a reach, flow from the 

underlying model cell to the stream is computed by 

(Prudic, 1989, p. 7):

(2)

where

/j

Cs = streambed conductance (L /T), 

Hs = average stream stage in the reach (L), and 

h = hydraulic head in the cell (L).

The average stream stage is determined from the 

topographic map of the Mirror Lake area. If the 

average stage of a reach is lower than the hydraulic
22

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

71°43'

43°58' —

1 KILOMETER 

I_________

0

I

1/2 1MILE

EXPLANATION

• Location of observed hydraulic head

-— No-flow boundary 

^^ Stream channel

43°57'

43°56' —

Base from U.S. Geological Survey

1:24,000 Woodstock, New Hampshire quadrangle

Figure 10. Model grid and boundary conditions.
MODEL OF STEADY-STATE GROUND-WATER FLOW 23

71°42'
24 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

EXPLANATION 

I No flow boundary

1 KILOMETER 

J_________

1/2

Vertical Exaggeration 5x

Figure 11. Section along row 25 of the model grid showing the five model layers.

1MILE

40

72

0 100

I T I

0 250 500 FEET

VERTICAL EXAGGERATION 5X

200 METERS

Figure 12. Section along row 50 of the model grid showing glacial deposits, Mirror Lake, lacustrine mud 

sediments, and bedrock in layers 1, 2, and 3.
MODEL OF STEADY-STATE GROUND-WATER FLOW 25

head in the underlying model cell, ground water 

discharges into the stream (Qs is negative). 

Conversely, if the average stage of a reach is higher 

than the hydraulic head in the underlying model cell, 

stream water infiltrates into the ground-water system 

(Qs is positive), providing there is streamflow in the 

reach. After Qs is computed for a reach, it is added to 

the streamflow from the adjacent upstream reach, and 

the sum is routed to the adjacent downstream reach. 

Because the streamflow simulated in this manner is 

derived from ground-water discharge, it represents 

baseflow in the stream.

In theory, the streambed conductance can be 

determined by (Prudic, 1989, p. 7)

(3) B

where

Ks = hydraulic conductivity of the stream bed (L/T), 

Bs = thickness of the streambed (L), 

Ws = width of the stream reach (L), and 

Ls = length of the stream reach (L).

In the present study, Ks and Bs are not known, so 

equation 3 cannot be directly applied. However, work 

of Shattuck (1991) indicates that there is little 

resistance to flow between the glacial deposits and the 

streams in the Mirror Lake watershed; that is, 

streambed conductance is high. Therefore, the 

streambed conductance is set to a large value for all 

stream reaches. This approach causes hydraulic heads 

in the cells below the streams to be nearly the same as 

the stream stage, in areas of the model where the 

computed Qs is not zero.

PARAMETER STRUCTURE

Distribution of hydraulic conductivity in the 

model is represented by dividing the model domain 

into zones, each having homogeneous hydraulic 

conductivity. The zones and hydraulic conductivities 

together define the parameter structure. Four different 

parameter structures are considered for the groundwater flow model of the Mirror Lake area. Model 

calibration is carried out for each parameter structure.

The simplest parameter structure is referred to 

as parameter structure A (table 5). The zones for this 

parameter structure correspond to the three hydrogeoTable 5. Hydraulic conductivities used in parameter structures A, B, C, and D.

[Hydraulic conductivity in meters per second; estimated, value estimated during model calibration; —, not used in parameter structure]

Hydraulic conductivity 

parameter

Value in parameter structure

B D

LGD

^GOtf

^GD,V

Glacial Deposits (excluding sand and gravel south of Mirror Lake)

estimated ~ estimated 

estimated

estimated

KSG

Sand and Gravel south of Mirror Lake

5 x 10'5 5 x 10'5 5 x ID'5

5 x 10'5

^LM

Lacustrine Mud Sediments

1 x 10'9 1 x 1(T9

1 x 10'9 1 x ID'9

Bedrock

*BR

KBR,3

KBR,4

KBR,5

KBR,U

KBR.L

estimatea esumatea

estimated

estimated

estimated

-

--

—

-

--

-

estimated

estimated
26 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

logic units in the Mirror Lake area (glacial deposits, 

lacustrine mud, and bedrock), with one exception. 

Because the sand and gravel south of Mirror Lake 

(fig. 2) exert an important control on the subsurface 

seepage from the lake, this deposit is considered a 

separate zone from the rest of the glacial deposits. 

Therefore, the zonation scheme for parameter 

structure A divides the model domain into four zones: 

glacial deposits, sand and gravel, lacustrine mud, and 

bedrock. The hydraulic conductivity of each zone is 

assumed to be isotropic. The hydraulic conductivities 

in parameter structure A are: 

KGD = hydraulic conductivity of glacial deposits,

excluding sand and gravel south of Mirror

Lake (L/T), 

KSG = hydraulic conductivity of sand and gravel south

of Mirror Lake (L/T), 

KLM ~ hydraulic conductivity of lacustrine mud

sediments (L/T), and 

KBR = hydraulic conductivity of bedrock (L/T).

Because KSG and KLM characterize relatively small 

parts of the model domain, it is unlikely that the 

hydraulic-head and flow data provide enough information to support estimation of these model parameters. 

Thus, their values are assumed to be known. The 

value of KSG is fixed at 5 x 10~5 m/s, which is the value 

used by Rosenberry and Winter (1993) to calculate 

subsurface seepage from Mirror Lake. The value of 

KLM is fixed at 1 x 10~9 m/s, which is within the range 

estimated by Rosenberry and Winter (1993). During 

model calibration, the values of KGD and KBR are 

estimated along with recharge from precipitation.

Although parameter structure A offers the 

appeal of simplicity, it might not adequately represent 

the hydraulic conductivity in the model domain. To 

investigate parameter structures of greater complexity, 

three additional cases are considered. They are 

referred to as parameter structures B, C, and D 

(table 5).

Parameter structure B is identical to A except 

the glacial deposits (excluding sand and gravel south 

of Mirror Lake) are assumed anisotropic. The 

components of hydraulic conductivity of glacial 

deposits in parameter structure B are: 

KGDH = horizontal hydraulic conductivity of glacial 

deposits, excluding sand and gravel south 

of Mirror Lake (L/T), and

KGDV = vertical hydraulic conductivity of glacial

deposits, excluding sand and gravel south

of Mirror Lake (L/T).

Fine-grained layers in the glacial deposits might 

impede vertical flow of water, therefore, the vertical 

hydraulic conductivity is assumed to be less than the 

horizontal hydraulic conductivity. On the basis of 

estimates by Rosenberry and Winter (1993), KGD y is 

assumed to be one-tenth of KGD>H. During model 

calibration, the value ofKGDH is estimated along with 

KBR and recharge from precipitation.

Parameter structure C is designed to investigate 

vertical variability of bedrock hydraulic conductivity. 

The modeled part of bedrock (uppermost 150 m) is 

vertically subdivided into three zones corresponding 

to model layers 3,4 and 5. The bedrock hydraulic 

conductivities in parameter structure C are: 

KBR 3 = hydraulic conductivity of bedrock in model

layer 3 (L/T), 

KBR4 = hydraulic conductivity of bedrock in model

layer 4 (L/T), and 

KBR 5 = hydraulic conductivity of bedrock in model

layer 5 (L/T).

The hydraulic conductivity zonation in model layers 1 

and 2 is identical to that of parameter structure A. 

Isotropy of hydraulic conductivity is assumed in all 

zones. During model calibration, the values ofKBRj, 

KBR4 and KBR5 are estimated along with KGD and 

recharge from precipitation.

Parameter structure D is designed to investigate 

the lateral variability of bedrock hydraulic conductivity. This parameter structure was prompted by data 

from the newly drilled well FS7 on the ridge at the 

northwestern corner of the study area (fig. 1). At this 

well site, the water table was found to be about 15 m 

below land surface. The shallow water table on the 

ridge indicates that the water table lies close to land 

surface over the entire study area. As shown by 

topographic contours in figure 1, upper hillsides (at 

altitudes greater than 360 m above sea level) are 

significantly steeper than lower hillsides (at altitudes 

less than 360 m above sea level) in the study area. 

Therefore, the water-table in upper hillsides is 

expected to be significantly steeper than the watertable in lower hillsides. Using uniform recharge from 

precipitation and homogeneous hydraulic conductivity 

in the bedrock, it is not possible to simulate a water
MODEL OF STEADY-STATE GROUND-WATER FLOW 27

table that has an abrupt change in slope. Therefore, in 

parameter structure D, the bedrock is laterally 

subdivided into two zones (fig. 13): (1) a zone beneath 

upper hillsides and hilltops (stippled), and (2) a zone 

beneath lower hillsides and valleys (not stippled). In 

plan view, the boundary separating the two zones 

approximately coincides with the 360-m contour line. 

This boundary extends vertically throughout the 

bedrock (model layers 3, 4 and 5). The bedrock 

hydraulic conductivities in parameter structure D are:

KBRU = hydraulic conductivity of bedrock beneath

upper hillsides and hilltops(L/T), and 

KBR,L = hydraulic conductivity of bedrock beneath

lower hillsides and valleys (L/T). 

The hydraulic conductivity zonation in model layers 1 

and 2 is identical to that of parameter structure A. 

Isotropy of hydraulic conductivity is assumed in all 

zones. During model calibration, the values ofKBRU 

and KBRL are estimated along with KGD and recharge 

from precipitation.

71°43'

43°58' —

43°57' —,

43°56'

Base from U.S. Geological Survey

1:24,000 Woodstock, New Hampshire quadrangle

Land surface elevations converted from feet to meters

EXPLANATION

[""I ZONE BENEATH UPPER HILLSIDES AND HILLTOPS

1 KILOMETER 

I________

1/2 1 MILE

CONTOUR INTERVAL 40 METERS

Figure 13. Delineation of the bedrock zone beneath upper hillsides and hilltops used in parameter 

structure D.
28 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

MODEL CALIBRATION

Model calibration (or parameter estimation) is 

the process of determining unknown model parameters so that simulated quantities match observed 

quantities (calibration data) to an acceptable level. 

The parameters determined in this study are (1) rate of 

recharge from precipitation, (2) hydraulic conductivity 

of glacial deposits, and (3) hydraulic conductivity of 

bedrock. Calibration is accomplished by using the 

computer program MODFLOWP (Hill, 1992), which 

performs nonlinear least-squares regression. Included 

in the calibration process is an assessment of 

parameter uncertainty and goodness of model fit. 

Calibration is carried out for the four parameter 

structures A, B, C, and D. These calibrations are 

referred to as calibrations A, B, C, and D.

CALIBRATION DATA

The two data sets used for model calibration in 

this study are (1) long-term average baseflow in 

Streams W, NW, and E at their inlets to Mirror Lake 

(table 2), and (2) long-term average hydraulic heads 

measured in piezometers and wells (tables 3 and 4). 

For brevity, these two data sets are referred to as 

(1) observed baseflow and (2) observed hydraulic 

head. Because well FS7 was drilled near the end of 

this study, when calibrations A, B, and C were already 

completed, hydraulic head measured in the new well 

was included only in calibration D.

As noted previously, the measurement location 

for hydraulic head is taken to be the midpoint of a 

piezometer screen, or the position of the most 

transmissive fracture in a well interval. During model 

calibration, observed hydraulic heads are compared 

with simulated hydraulic heads at model cells that 

contain the measurement locations. If two or more 

measurement locations lie in the same model cell, the 

observed hydraulic heads are first averaged and then 

compared with the simulated hydraulic head. If a 

measurement location lies in a model cell that has 

become unsaturated during model calibration, the 

observed hydraulic head is compared with the 

simulated hydraulic head in the model cell one layer 

below the unsaturated cell.

NONLINEAR LEAST-SQUARES 

REGRESSION METHOD

The nonlinear least-squares regression method 

used for model calibration is based on the formulation 

of Cooley and Naff (1990). This formulation assumes 

that a true model exists, and errors (or disturbances) 

are random, unbiased, and normally distributed 

(Cooley and Naff, 1990, p. 164-165). In principle, the 

variance-covariance matrix for the errors can be used 

to weight the observations. In the present study, this 

matrix is unknown, and weighting is assigned in a 

subjective manner on the basis of knowledge of the 

field data.

Parameters to be estimated (recharge from 

precipitation and hydraulic conductivities) are denoted 

by the following notation: 

bi = /* parameter to be estimated, 

b = vector of parameters to be estimated, 

b- = optimal estimate of the Ith parameter, and 

b = vector of optimal parameters. 

In addition, 

p = number of parameters to be estimated

(dimensionless), 

n = total number of observations of hydraulic head

and baseflow (dimensionless). 

In the actual regression, the common log of hydraulic 

conductivity is estimated; the result is then exponentiated to yield the estimate for hydraulic conductivity. 

For example, in calibration A, the parameter vector b 

consists of \og(KGD\ \og(KBR), and recharge from 

precipitation. In this case, p equals 3.

The optimal parameters are those that minimize 

the sum of squared weighted differences between 

observed and simulated quantities:

S(b) = (4)

/= 1 i = 1

where

Hi = observed hydraulic head at location i (L), 

hj(b) = simulated hydraulic head at location / computed

with parameter vector b, (L), 

nh = number of hydraulic head observations

(dimensionless), 

Qi = observed baseflow in stream i (L3/T),
MODEL CALIBRATION 29

qf(b) = simulated baseflow in stream i computed using

parameter vector b (L3/T), 

nq = number of streams (dimensionless), and

w;1/2 = weight on difference between observed and

simulated baseflow in stream / (T/L2). 

In calibrations A, B and C, nh is 90. In calibration D, 

nh is 91. The value of nq is 3 in all cases. The total 

number of observations, «, is the sum of n^ and nq.

Because hydraulic head and baseflow are 

different quantities having different units, weighting is 

necessary to combine both into a single optimization 

criterion (equation 4). For Streams W and NW, the 

weights are chosen such that a 1 percent difference 

between observed and simulated baseflow is equivalent to a 1 m difference between observed and 

simulated hydraulic head. For Stream E, the weight is 

such that a 5 percent difference between observed and 

simulated baseflow is equivalent to a 1 m difference 

between observed and simulated hydraulic head. The 

unequal weighting is appropriate because Stream E 

drains a much smaller sub-basin than Streams W or 

NW, and therefore exerts a much smaller influence on 

the overall ground-water system. The unequal 

weighting allows a larger tolerance when matching 

observed and simulated baseflow to the E Stream than 

to Streams W and NW.

The minimization of S(b) is accomplished with 

the computer program MODFLOWP (Hill, 1992), 

which uses a modified Gauss-Newton method to 

determine the optimal MODFLOW parameters. As 

noted previously, the model in this study requires that 

at each cell in the horizontal dimension, the computed 

water table be allowed to lie in any of the five model 

layers. Consequently, the saturated thickness of each

model layer is determined during the model simulation. However, the version of MODFLOWP used in 

this study does not allow the water-table position to 

vary during a model simulation, and thus requires that 

the saturated thickness of each layer be specified. To 

overcome this problem, an iterative scheme is used. 

Starting with a set of assumed parameters, 

MODFLOW is executed, the water table position is 

computed during the simulation, and saturated 

thicknesses are computed and stored. Next, a 

MODFLOWP calibration is run with the saturated 

thicknesses set equal to those obtained from the prior 

MODFLOW simulation. If the parameters estimated 

by MODFLOWP differ significantly from those used 

in MODFLOW, the procedure is repeated. In the next 

iteration, parameters estimated by MODFLOWP are 

used in MODFLOW to obtain new saturated 

thicknesses. The new saturated thicknesses are then 

used in MODFLOWP to obtain new parameter 

estimates. This iterative process continues until the 

parameter estimates no longer change from one 

iteration to the next.

OPTIMAL ESTIMATES OF RECHARGE

FROM PRECIPITATION AND HYDRAULIC

CONDUCTIVITY

Examination of the results of model calibration 

(table 6) shows that, for many model parameters, the 

four calibrations A, B, C, and D yield optimal 

estimates that are similar to each other. The four 

optimal estimates of recharge from precipitation all 

fall within a narrow range from 26 to 28 cm/year.

Table 6. Optimal estimates of recharge from precipitation and hydraulic conductivities.

[Recharge in centimeters per year; hydraulic conductivity in meters per second; —, not used in parameter structure]

Optimal estimate in calibration

Parameter A

Recharge from precipitation 27

KGD 2.67 x 10'6

KGD>H

r 1 KGD,V

KBR 2.93 x 10'7

KBR, 3

KBR, 4

KBR, 5

KBR, u

KBRL

B

26

--

2.57 x 10'6

2.57 x 10'7

3.05 x 10'7

--

--

--

-

--

C

27

2.54 x 10'6

—

--

--

3.10 xlO'7

2.89 x 10'7

2.78 x 10'7

--

--

D

28

1.74xlO'6

—

--

--

-

--

--

6.26 x 10'8

3.21 x 10'7

In parameter structure B, KGDV is preset to equal one-tenth of KGDtH.
30 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

These values represent approximately 20 percent of 

the long-term average precipitation measured at the 

Forest Service station and are near the lower limit of 

the range (16 to 45 percent) estimated by Mau and 

Winter (1997).

In the parameter structures where glacial 

deposits are assumed isotropic (calibrations A, C and 

D), the optimal estimates of KGD are also similar to 

each other, ranging from 1.7 x 10~6 to 2.7 x 10~6 m/s. 

These values lie between the slug-test estimates for till 

and for sand and gravel as determined by Winter and 

Rosenberry (table 1). If glacial deposits are assumed 

anisotropic (calibration B), the optimal estimate of 

KGDH is similar to those of KGD. The value of KGD v 

was preset to equal one-tenth of KGD H.

Optimal estimates of bedrock hydraulic conductivity also show similarities. In the parameter 

structures where the bedrock is assumed uniform 

(calibrations A and B), the optimal estimates of >~l

KBR are both about 3 x 10 m/s. This value is 

consistent with the effective hydraulic conductivity of

2 x 10"7 m/s determined by Hsieh and Shapiro (1994) 

for the bedrock underlying the FSE well field. If the 

bedrock is vertically subdivided by model layers 

(calibration C), the optimal estimates of KBR 3 , KBR 4 

and KBR5 are all about 3 x 10"7 m/s. The absence of a 

significant trend in these values indicates that the 

calibration data in the present study do not support the 

hypothesis that hydraulic conductivity varies with 

depth in the modeled part (uppermost 150 m) of 

bedrock. If the bedrock is laterally subdivided 

according to land-surface altitude (calibration D), 

the optimal estimate of KBR^ is similar to those of 

KBR, but the optimal estimate of KBRU is about 

one-fifth of KBR^. These results indicate that

3 x 10"7 m/s is a reasonable hydraulic conductivity for 

bedrock beneath lower hillsides and valleys, but a 

lesser value (6 x 10"8 m/s) might be more appropriate 

for bedrock beneath upper hillsides and hilltops. 

However, because the estimation of KBR u is based 

largely on the single hydraulic head measured at the 

newly drilled well FS7, additional data are required to 

further confirm this result.

PARAMETER UNCERTAINTY

Within the context of regression analysis, 

parameter uncertainty can be assessed by calculating a 

confidence interval for each parameter. These

confidence intervals should be interpreted with 

caution, because they are computed under the assumption that errors can be approximated by residuals 

(differences between observed and simulated quantities). In this approach, model attributes such as 

geometry, zonation, boundary conditions, and values 

of non-estimated parameters (table 5) are assumed 

perfectly known. Errors or uncertainties in these 

model attributes could greatly increase the uncertainty 

of estimated parameters.

The confidence intervals computed in this study 

are known as approximate, individual confidence 

intervals. The intervals are approximate because the 

computation assumes that the model is effectively 

linear in the vicinity of the optimal parameters. The 

intervals are individual in that they are computed for 

each parameter while holding all other parameters at 

their optimal values. In addition, the computation 

assumes that errors are independent and normally 

distributed. Because the above assumptions might not 

be satisfied in calibrations of the ground-water model 

for this study, their validity is examined.

To calculate confidence intervals, it is necessary 

to first calculate the vaiiance-covariance matrix for the 

parameters, which is (Cooley and Naff, 1990, p. 167, 

eq. 5.4-13)

V(b) = (5)

where f is defined as (Cooley and Naff, 1990, p. 166, 

eq. 5.4-5)

,2 _ S (b)

3 *~ "n-p (6)

where S (b) is computed using equation 4 and the 

optimal parameters b, co is the n by n diagonal matrix 

of weights (which equals one for hydraulic heads and 

w^2 for baseflow in stream /), and X is the n by p 

sensitivity matrix containing derivatives of simulated 

hydraulic head or baseflow at the measurement 

locations with respect to model parameters b, 

evaluated at the optimal parameters b. The i 

element on the main diagonal of V (b) is the 

variance, V#, for the ith parameter. After the variancecovariance matrix is computed, the approximate, 

individual, 100(1-a) percent confidence interval for 

the ith parameter is given by (Hill, 1992, page 58)
MODEL CALIBRATION 31

(7)

where t(n-p,l-a/2) is the 100(1-0/2) percentage point 

of the Student-t distribution with n-p degrees of 

freedom.

The approximate, individual, 95-percent 

confidence intervals for recharge from precipitation 

and hydraulic conductivities are shown in figure 14. 

In all four calibrations, the confidence intervals for 

recharge span approximately the same width — about 

15 percent of the optimal estimate — indicating that the 

four recharge estimates have similar levels of 

uncertainty. For hydraulic conductivities, confidence 

intervals in calibrations A, B and D are much narrower 

than those in calibration C. This comparison indicates 

that hydraulic conductivity estimates are much better 

constrained in calibrations A, B and D, but they are 

subject to much greater uncertainty in calibration C. 

The relatively high uncertainty in calibration C is an 

additional indication that the calibration data do not 

support the hypothesis that hydraulic conductivity 

varies with depth in the modeled portion (uppermost 

150 m) of bedrock.

Because the confidence intervals in figure 14 are 

computed under the assumptions that (1) the model is 

effectively linear about the optimal parameters and (2) 

errors are independent and normally distributed, the 

validity of these assumptions should be examined. In 

this study, the degree of model nonlinearity is 

examined by the "modified Beale's measure" (Cooley 

and Naff, 1990, p. 187). Independence and normality 

of errors are inferred from analysis of weighted 

residuals, which are weighted differences between 

observed and simulated quantities (hydraulic heads 

and baseflows).

Modified Beale's measures calculated by the 

computer program BEALEP (Hill, 1994, p. 45-54) are 

shown in table 7. On the basis of criteria suggested by 

Cooley and Naff (1990, p. 189), these modified 

Beale's measures are interpreted as follows: the model 

is effectively linear if the modified Beale's measure is 

less than 0.03; the model is highly nonlinear if the 

modified Beale's measure is greater than 0.37. In 

calibrations A, B, and D, the modified Beale's 

measure lies between 0.03 and 0.37, indicating that 

although the model is nonlinear, it is not highly 

nonlinear. In calibration C, the modified Beale's

O

fe 

cc.

99

< O

D

0 10 20 30 40 50 

RECHARGE, IN CENTIMETERS PER YEAR

B

o 

fe

GC

99 —i 

< O

D

KBR~

KBR~

Tf i

BR,3

BR,4

kBR,5

K *-** 

^BR,L

K w BR,U

MKGD

" KGD

' ——— " KGD

"

"KGD

1CT8 10~7

HYDRAULIC CONDUCTIVITY, 

IN METERS PER SECOND

1(T5

Figure 14. Approximate, individual, 95-percent 

confidence intervals for (A) recharge from precipitation 

and (B) hydraulic conductivities.
32 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

measure is significantly larger than 0.37, indicating a 

high degree of nonlinearity.

Table 7. Modified Beale's measure, RN2, and critical value ofRN2 

[All statistics dimensionless]

Statistic

Modified Beale's

measure 

R 2 KN

Critical value of RN2

Value in calibration

A

0.14

0.925 

0.973

B

0.22

0.919 

0.973

C

21.7

0.943 

0.973

D

0.26

0.899 

0.974

In the analysis of weighted residuals, the 

residual for hydraulic head is assigned a weight of 

unity, and the residual for baseflow in stream / is 

assigned the weight w,- 1 /2. To examine normality, 

weighted residuals for each parameter structure are 

plotted on normal probability plots (fig. 15). For 

parameter structures A and C, the weighted residuals 

form a nearly straight line on the normal probability 

plots. The plots for parameter structures B and D each 

contain one weighted baseflow residual that is an 

outlier, but the remainder of the weighted residuals for 

these two parameter structures plot nearly on a straight 

line. Therefore, the normal probability plots indicate 

that the distribution of weighted residuals for all the 

parameter structures is close to normal. To examine 

both normality and independence, a hypothesis test 

suggested by Hill (1992, p. 63) is used. This test is 

based on the correlation coefficient (denoted by RN ) 

between the weighted residuals ordered from smallest 

to largest and the order statistics from the normal

^

probability distribution function. If the RN falls 

below a critical value (given by Hill, 1992, p. 64) the 

hypothesis that the weighted residuals are independent 

and normally distributed is rejected. Table 7 shows 

that for all four parameter structures, RN lies close to, 

but nonetheless falls below, the critical value. 

Therefore, strictly speaking, the hypothesis of 

independence and normality is rejected. However, if 

the normal probability plots (fig. 15) are considered 

together with the hypothesis test, the overall results 

indicate that the weighted residuals (and by inference, 

the errors) might be close to independent and normally 

distributed.

The above analyses of the modified Beale's 

measure and weighted residuals indicate that the 

assumptions of linearity, independence, and normality

are not strictly valid. Nonetheless, in calibrations A, 

B, and D, these assumptions appear nearly satisfied. 

Thus the computed confidence intervals could be 

considered rough indicators of parameter uncertainty. 

In calibration C, the linearity assumption appears 

seriously incorrect. Thus, the computed confidence 

intervals could be subject to large inaccuracies.

GOODNESS OF MODEL FIT

The overall goodness of model fit can be 

assessed by examining statistics of residuals (observed 

values minus simulated values). In general, a lower 

residual statistic indicates a better overall model fit. 

The sum of squared weighted residuals, S (b),

/v

defined by equation 4 with b being the vector of 

optimal parameters, gives an overall measure of model 

fit to all observations. In addition, the sum of squared 

residuals of hydraulic head, defined as

SH= I [",-¥*>] '

/= 1

gives an overall measure of model fit to hydraulic 

head. The sum of squared weighted residuals of 

baseflow, defined as

sq =

<*>

(9)

/ = i

gives an overall measure of model fit to baseflow. 

Finally, the standard error of hydraulic head, defined 

as

/2

sh = (10)

gives the average (root-mean-square) discrepancy 

between observed and simulated hydraulic heads. The 

standard error of baseflow cannot be calculated 

because of insufficient data for statistical analysis.

Computed values of residual statistics for 

calibrations A through D, and the maximum and 

minimum residual of hydraulic head in each calibra-
MODEL CALIBRATION 33

Q 

co

LLJ 

CC

Q 

LU

LU

20

18

16

14

12

10

8

6

4

2

0 

_2

-4

-6

-8 

20

18

16

14

12

10

8

6

4

2

0

-2

-4

-6

-8

\^ \ n rn i 11 n i r

PARAMETER STRUCTURE A

WEIGHTED 

BASEFLOW 

RESIDUAL

HYDRAULIC _ 

HEAD • O 

RESIDUAL

q 

c>

PARAMETER STRUCTURE C

i \ \ \ \\\ n i i r

PARAMETER STRUCTURE B

i i i i i i

i i i i n i ni r

PARAMETER STRUCTURE D

to

ir>o o

CO

o o o in

O5

O5

oi

O5

o o o 

co in r--

O IT) 

O5 O5

O5 

O5

O5

O5

O) 

O5

O5 

O5

O) 

O5

EXPECTED CUMULATIVE PROBABILITY, IN PERCENT

Figure 15. Normal probability plots of weighted residuals.

tion are given in table 8. Because calibration D 

includes an extra hydraulic head measured at the 

newly drilled well FS7, it is considered separately 

from the other three calibrations. Comparison of 

residual statistics for calibrations A, B, and C shows 

that S (b) is highest in B, and is the same in A and C. 

The poorer fit in calibration B is due to a poorer fit to 

baseflow. The fit to hydraulic head is about the same 

in all three calibrations. The standard error of 

hydraulic head in all three cases is about 2.8 m, which 

is relatively small compared with the 145 m of 

elevation difference between the highest and lowest

observed hydraulic head. Observed and simulated 

baseflow differ by less than 10 percent in Streams W 

and NW (table 9), which together account for more 

than 90 percent of total baseflow. These results lead to 

two conclusions. First, calibrations A, B, and C all 

yield relatively good overall model fits. Second, the 

model fit in calibration A is not improved by 

introducing a 10-to-l ratio of horizontal-to-vertical 

anisotropy in hydraulic conductivity of glacial 

deposits (calibration B), or by varying hydraulic 

conductivity with depth in the modeled part 

(uppermost 150 m) of the bedrock (calibration C).
34 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

Table 8. Residual statistics used to measure overall goodness of model fit.

[Sum of squared weighted residuals and sum of squared residuals in square meters; standard error, maximum residual and minimum residual of hydraulic 

head in meters]

Statistic

Sum of squared weighted residuals, S(b)

Sum of squared residuals of hydraulic head, S^

Sum of squared weighted residuals of baseflow, Sq

Standard error of hydraulic head, s^

Minimum residual of hydraulic head3

Maximum residual of hydraulic head3

Value in calibration

A1

799

665

134

2.76

-4.0

10.0

B1

1017

677

340

2.79

-4.9

9.8

C1

799

641

158

2.75

-3.9

9.7

D2

1048

636

412

2.70

-4.7

9.5

ldoes not include hydraulic-head residual at well FS7. 

Includes hydraulic-head residual at well FS7.

3Positive residual indicates that observed value is higher than simulated value. Negative residual indicates that observed value is lower than 

simulated value.

Table 9. Observed and simulated baseflow in Streams W, NW, and E. 

[Baseflow in cubic meters per year]

Baseflow simulation with parameter structure

W 63,000

NW 100,000

E 5,000

Total 168,000

A

60,200

97,100

7,300

164,600

B

63,500

92,300

8,800

164,600

C

59,500

96,700

7,300

163,500

D

62,800

92,300

9,200

164,300

If taken at face value, the residual statistics for 

calibration D show that the overall model fit in calibration D is no better (in fact, a bit worse) than the overall 

model fit in calibration A. Such an assessment is 

correct as long as it is confined to lower hillsides and 

valleys, that is, the region in which hydraulic heads are 

measured for calibration A. However, calibration D 

includes an extra hydraulic-head measured at the 

newly drilled well FS7, located near the highest point 

of the study area. The simulated hydraulic head at 

FS7 is within 1 m of the observed hydraulic head. In 

this regard, calibration D yields a more realistic 

simulation than calibration A at higher altitudes. This 

conclusion, however, should be considered preliminary and requires further confirmation using additional 

observations of hydraulic heads at higher altitudes.

In addition to an overall assessment of model fit, 

a more detailed assessment can be made by examining

the spatial distribution of hydraulic-head residuals. 

Figure 16 shows this distribution in the five model 

layers over an area in which most of the piezometer 

and wells are located. Residuals are shown only for 

calibration A, because similar results are obtained for 

the other three calibrations. Although most of the 

residuals in figure 16 are relatively small (many fall 

within the range from -4 m to 4 m), their spatial distribution in model layers 1 and 2 (glacial deposits) shows 

an obvious pattern. In these two layers, residuals in 

the area northwest of Mirror Lake are generally 

positive (observed values higher than simulated 

values), and residuals in the area south of Mirror Lake 

are generally negative (observed values lower than 

simulated values). This pattern suggests that residuals 

in layers 1 and 2 are not random, but that they exhibit 

trends due to limitations imposed by model assumptions of uniform recharge from precipitation and
SIMULATION OF GROUND-WATER FLOW 35

homogeneous hydraulic conductivity in each 

hydrogeologic unit. Each of model layers 3,4, and 5 

contains too few residuals to draw inferences about 

their spatial distribution.

The lack of fit between simulated and observed 

hydraulic heads at the local scale is further illustrated 

in figure 17, which shows hydraulic-head contours in 

the vertical section AB (see fig. 8) on a hillside west of 

Mirror Lake. In general, simulated contours match 

observed contours in conveying the sense of lateral 

ground-water flow in the downslope direction. 

However, detailed examination of the two sets of 

contours shows distinct differences. The observed 

contours have an irregular pattern, suggesting that 

flow is locally affected by nonuniform recharge and/or 

heterogeneous hydraulic conductivity. The simulated 

contours have a more regular pattern, because of the 

model assumptions of uniform recharge and homogeneous hydraulic conductivity. This comparison further 

shows that the model developed in this study characterizes the generalized ground-water flow system, not 

local details.

SIMULATION OF 

GROUND-WATER FLOW

The major objective of simulating ground-water 

flow in the Mirror Lake area is to determine the threedimensional extent of the Mirror Lake ground-water 

basin and the ground-water budget for the basin. 

Because four parameter structures are used in model 

calibration, four sets of simulation results are 

obtained. In this study, it is found that simulations 

using parameter structures A, B, and C yield similar 

results. Therefore, only simulations using parameter 

structures A and D are presented. These two simulations are referred to as simulation A and simulation D. 

Although simulation D is expected to be more realistic 

than simulation A at higher altitudes of the Mirror 

Lake area, both simulations are presented for comparison.

WATER TABLE

The water-table configuration for simulation A 

is shown in figure 18, and that for simulation D is 

shown in figure 19. The two water-table configurations are similar at lower altitudes, but differ significantly at upper altitudes. At the northwestern corner 

of the model domain, the water table in simulation A 

rises to 560 m above sea level (about 100 m below 

land surface), whereas the water table in simulation D 

rises to 640 m above sea level (about 20 m below land 

surface). The shallower water table in simulation D 

results from using a lower hydraulic conductivity for 

bedrock beneath upper hillsides in order to match the 

newly observed hydraulic head at well FS7. This new 

observation was not used in the calibration process for 

simulation A.

Areas where the simulated water table lies in the 

glacial deposits are also shown in figures 18 and 19. 

In these areas, recharge from precipitation enters the 

glacial deposits. Outside these areas, the simulated 

water table lies in the bedrock, and recharge from 

precipitation enters the bedrock, either by infiltration 

through the overlying unsaturated glacial deposits, or 

directly from land surface where bedrock crops out. 

Areas where the glacial deposits are absent or unsaturated occur mostly on upper hillsides and away from 

streams. Because the water table is deeper in simulation A than in simulation D, unsaturated glacial 

deposits cover a larger area in simulation A. 

Consequently, more of the water table lies in bedrock 

in simulation A than in simulation D, and therefore 

more recharge from precipitation enters the bedrock in 

simulation A than in simulation D. Conversely, less 

recharge from precipitation enters the glacial deposits 

in simulation A than in simulation D.

MIRROR LAKE GROUND-WATER BASIN

The Mirror Lake ground-water basin refers to 

the volume of subsurface through which ground water 

flows from the water table to Mirror Lake or its three 

inflow streams. The three-dimensional extent of the 

basin is determined by tracking ground-water flow 

paths in a reverse direction from Mirror Lake and its 

inflow streams back to the water table. The volume 

spanned by these flow paths delineates the Mirror 

Lake ground-water basin. In this study, flow paths are 

computed with the computer program MODPATH 

(Pollock, 1994), a particle tracking post-processor 

package for MODFLOW.

In the MODPATH simulations, the starting 

positions of particles are chosen so that ground-water 

flow paths entering the sides and bottom of Mirror
36 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

Model Layer 1

Model Layer 2

EXPLANATION

Hydraulic-head residual 

in calibration A equals:

• -4 to -2 meters

• -2 to 0 meters

A 0 to 2 meters

A 2 to 4 meters 

A 4 to 6 meters 

/\ 6 to 10 meters

Figure 16. Distribution of hydraulic head residuals in model layers 1 through 5 in 

calibration A.

Negative residual indicates that the observed 

value is less than the simulated value; 

positive residual indicates that the observed 

value is greater than the simulated value.
SIMULATION OF GROUND-WATER FLOW 37

Model Layer 3

Model Layer 4

Model Layer 5
38 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

FS3

FS1

FS2

FSE4

50 100 METERS 

J_____I

I I F

0 150 300 FEET

No Vertical Exaggeration

EXPLANATION

-240— OBSERVED HYDRAULIC-HEAD CONTOURContour interval 5 meters. Datum is sea level.

— 240— HYDRAULIC-HEAD CONTOUR SIMULATED WITH 

PARAMETER STRUCTURE A-Contour interval 

5 meters. Datum is sea level.

• HYDRAULIC-HEAD MEASUREMENT LOCATION-Shows 

location of the highest-transmissivity fracture in a 

well interval.

Figure 17. Observed and simulated hydraulic head in section A-B. (See figure 8 for location of line of 

section.)
SIMULATION OF GROUND-WATER FLOW 39

71°43' 71°42' 71°41f

43°58' —

43°57' —

43°56' —

Base from U.S. Geological Survey

1:24,000 Woodstock, New Hampshire quadrangle

EXPLANATION 

Results of Simulation A:

0 1/2 1 MILE 

CONTOUR INTERVAL 40 METERS

WATER-TABLE CONTOUR - Shows altitude of simulated water table. Interval 20 meters.

Datum is sea level.

AREA WHERE WATER TABLE LIES IN GLACIAL DEPOSITS

~^~\ Stream channel with streamflow 

"""•••. Stream channel without streamflow

Figure 18. Water table for simulation A, and areas where water table lies in glacial deposits.
40 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

71°43'

43°58' -

43°5T -

43°56' -

Base from U.S. Geological Survey

1:24,000 Woodstock, New Hampshire quadrangle

EXPLANATION 

Results of Simulation D:

0 1/2 1 MILE 

CONTOUR INTERVAL 40 METERS

WATER-TABLE CONTOUR - Shows altitude of simulated water table. Interval 20 meters.

Datum is sea level.

AREA WHERE WATER TABLE LIES IN GLACIAL DEPOSITS

Stream channel with streamflow 

Stream channel without streamflow

Figure 19. Water table for simulation D, and areas where water table lies in glacial deposits.
SIMULATION OF GROUND-WATER FLOW 41

Lake and its inlet streams are tracked backwards. 

Twenty-five particles are placed in each finite-difference cell representing Mirror Lake and in each cell 

containing a reach of Stream W, NW, or E. Five 

particles are placed on each side face and on the 

bottom face of these cells. One particle is placed in 

the center of the face, and one is placed near each 

corner of the face.

The horizontal extent of the Mirror Lake 

ground-water basin for simulation A is shown in 

figure 20, and the three-dimensional extent is shown in 

figure 21. Corresponding illustrations for simulation 

D are shown in figures 22 and 23. In both simulations, 

the ground-water basin extends far up the hillside in 

the northwestern part of the study area (fig. 20 and 22), 

and ground water flows across the topographic divide 

that forms the northwestern boundary of the surfacewater basin. At the southern end of Mirror Lake, the 

area of subsurface seepage from the lake is not 

included in the Mirror Lake ground-water basin 

because ground water in this area flows towards 

Hubbard Brook.

In figures 20 and 22, the horizontal extent of the 

Mirror Lake ground-water basin is divided into a 

"recharge zone" (heavily shaded area) and two 

"underflow zones" (lightly shaded areas). The 

recharge zone is an area where the top of the groundwater basin is coincident with the water table. 

Recharge entering this area eventually drains to Mirror 

Lake or its inlet streams. In both simulations A and D, 

the size of the recharge zone is about 1.5 times the size 

of the surface-water basin. In other words, the area 

from which ground-water recharge drains to Mirror 

Lake or its inlet streams is about 1.5 times the area of 

surface-water drainage to the lake. The underflow 

zones shown on figures 20 and 22 are areas where the 

top of the ground-water basin lies below the water 

table. In these areas, the Mirror Lake ground-water 

basin lies beneath overlying ground-water basins.

In figures 21 and 23, the heavily shaded surface 

on the three-dimensional rendering of the groundwater basin is the water table. This surface 

corresponds to the heavily shaded area (recharge zone) 

on the overlying map. Flow paths in both simulations 

A and D show that a portion of recharge occurring on 

the northwestern part of the recharge zone flows at 

depth under Norris Brook to discharge into Mirror 

Lake or its inlet streams. In the shallow saturated zone 

around Norris Brook, recharge to the water table flows 

to Norris Brook to form a separate ground-water

basin. Therefore, part of the Mirror Lake groundwater basin lies beneath the adjacent ground-water 

basin that drains into Norris Brook. The lightly 

shaded surface on figures 21 and 23 is the boundary 

below the water table that separates the Mirror Lake 

ground-water basin and the overlying Norris Brook 

ground-water basin. No flow occurs across this 

surface. Below this surface, ground water flows to 

Mirror Lake and its inlet streams. Above this surface, 

ground water flows to Norris Brook. The horizontal 

extent of this surface corresponds to the larger of the 

two lightly shaded areas (underflow zones) on the 

overlying map. A smaller underflow zone occurs on 

the north side of the basin in each simulation (fig. 20 

and 22). The surface corresponding to this underflow 

zone is hidden from view in figures 21 and 23.

Two sets of additional MODPATH simulations 

are performed to examine the sensitivity of the 

position and extent of the Mirror Lake ground-water 

basin to variations in initial particle placement. First, 

the number of particles placed on the faces of cells 

containing Mirror Lake and Streams W, NW, and E is 

reduced, and second, particles are placed in a lattice 

arrangement within the cells instead of on the cell 

faces. In these simulations, the location and size of the 

Mirror Lake ground-water basin do not significantly 

differ from the location and size of the basin 

delineated in the simulation with five particles on each 

of the bottom and side faces.

GROUND-WATER BUDGET

The ground-water budget for the Mirror Lake 

ground-water basin gives an accounting of recharge to 

the basin, discharge from the basin, and flow between 

hydrogeologic units in the basin. These budget 

components are illustrated in figure 24. Recharge to 

the ground-water basin includes: 

RI = recharge from precipitation to bedrock (L3/T), 

R2 = recharge from precipitation to glacial deposits

(L3/T), and

/?3 = recharge from streams to glacial deposits (L3/T), 

which occurs along losing reaches of 

Streams W and NW.

Discharge from the ground-water basin includes: 

£>! = discharge from glacial deposits to streams

(LVD,

DI = discharge from glacial deposits to Mirror Lake 

(LTD,
42 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

71°43'

43°58' -

43°57' —

43°56' -T

Boundary 

MkrorLake 

.surfac^wa 

basin

Base from U.S. Geological Survey

1:24,000 Woodstock, New Hampshire quadrangle

Land surface elevations converted from feet to meters

EXPLANATION

Results of Simulation A:

0 1/2 1 MILE 

CONTOUR INTERVAL 40 METERS

RECHARGE ZONE - Area where top of Mirror Lake ground-water basin coincides with the water table 

UNDERFLOW ZONE - Area where top of Mirror Lake ground-water basin lies beneath the water table

Boundary of lateral extent of Mirror Lake ground-water basin 

*^**^» Boundary between recharge zone and underflow zone 

^^ Stream channel with streamflow

Stream channel without streamflow

Figure 20. Horizontal extent of Mirror Lake ground-water basin in simulation A. Mirror Lake surface-water 

basin is included for comparison.
SIMULATION OF GROUND-WATER FLOW 43

Figure 21. Three-dimensional 

extent of Mirror Lake groundwater basin in simulation A. 

See figure 20 for explanation.

D3 = discharge from lacustrine mud sediments to

Mirror Lake (L3/T), and 

D4 = discharge from bedrock to Mirror Lake (L3/T),

which occurs where bedrock crops out on

the east side of the lake

Flow between hydrogeologic units in the groundwater basin includes 

BI = flow from glacial deposits to bedrock (L3/T),

82 - flow from bedrock to glacial deposits (L3/T), 

and

83 = flow from bedrock to lacustrine mud sediments

(L3/T).

Under steady state, total recharge to the ground-water 

basin, defined as

(11)

equals total discharge from the basin, defined as

(12)

Budget components are computed by the computer 

program ZONEBUDGET (Harbaugh, 1990), which 

calculates subregional water budgets using simulation 

results from MODFLOW. Results are given in 

table 10.

In both simulations A and D, computed total 

recharge to the Mirror Lake ground-water basin (RT) 

is approximately 300,000 m /year. Nearly all the 

recharge is derived from precipitation (R^ + /?2); only

a small portion is derived from infiltration of stream 

water along losing reaches of streams (R>$). The 

recharge from precipitation enters the basin either in 

areas where the simulated water table lies in the 

bedrock (R{) or in areas where the simulated water 

table lies in the glacial deposits (/?2). Although in 

simulation A the amount of recharge entering the 

bedrock is slightly greater than that entering the glacial 

deposits, in both simulations A and D roughly half of 

the recharge enters each of the two hydrogeologic 

units.

Of the total discharge from the Mirror Lake 

ground-water basin (DT\ in both simulations A and D 

slightly more than half flows into Streams W, NW and 

E (Dj), slightly less than half flows into Mirror Lake 

from the glacial deposits (D2), and minute amounts 

enter Mirror Lake from the lacustrine mud (D3) and 

from the bedrock where it crops out on the east side of 

the lake (D4). In both simulations A and D, computed 

total ground-water seepage to Mirror Lake (D2 + D3 + 

D4) is approximately 130,000 m3/year. Although this 

value is greater than Rosenberry and Winter's (1993) 

estimate of 80,000 m3/year, it is within the range of 

uncertainty (±100 percent) for their estimate.

Within the Mirror Lake ground-water basin, 

flows between hydrogeologic units are strongly 

controlled by the direction of the vertical flow 

component. Flow from glacial deposits to bedrock 

(B{) occurs primarily beneath upper hillsides, where 

the vertical component of flow is predominantly 

downward. Because much of the glacial deposits on
44 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

43°58'

43°57

43°56'

Base from U.S. Geological Survey

1:24,000 Woodstock, New Hampshire quadrangle

Land surface elevations converted from feet to meters

EXPLANATION 

Results of Simulation D:

1/2 1 MILE 

CONTOUR INTERVAL 40 METERS

RECHARGE ZONE - Area where top of Mirror Lake ground-water basin coincides with the water table 

I UNDERFLOW ZONE - Area where top of Mirror Lake ground-water basin lies beneath the water table

** «* Boundary of lateral extent of Mirror Lake ground-water basin 

^""•••"w Boundary between recharge zone and underflow zone 

^"^ Stream channel with streamflow 

""- Stream channel without streamflow

Figure 22. Horizontal extent of Mirror Lake ground-water basin in simulation D. Mirror Lake surface-water basin 

is included for comparison.
SIMULATION OF GROUND-WATER FLOW 45

Table 10. Ground-water budget of Mirror Lake ground-water basin in simulations A and D. 

[recharge, discharge, and flow in cubic meters per year; —, not computed]

Budget component

Value in simulation Approximate, individual, 

95% confidence interval1

Precipitation to bedrock2 (R{) 

Precipitation to glacial deposits3 

Streams to glacial deposits4 (/?3) 

Total Recharge (RT = Rl+R2+R3)

Recharge to ground-water basin

172,000 147,000

129,000 145,000

6,000 5,000

307,000 297,000

161,000 to 183,000

120,000 to 138,000

5,000 to 7,000

Discharge from ground-water basin

Glacial deposits to streams (D])

Glacial deposits to Mirror Lake (D2)

Lacustrine mud sediments to Mirror Lake (D3)

Bedrock to Mirror Lake5 (D4)

Total Discharge (DT = 0^2+0^+0^

171,000

133,000

1,000

2,000

307,000

169,000

125,000

1,000

2,000

297,000

161,000 to 181,000

126,000 to 140,000

-

-

--

Flow between hydrogeologic units in ground-water basin

Glacial deposits to bedrock (Bi)

Bedrock to glacial deposits (B2)

Bedrock to gyttja (B3)

18,000

187,000

1,000

34,000

178,000

1,000

13,000 to 23,000

174,000 to 200,000

-

^or simulation A only, 

in areas where simulated water table lies in bedrock.

o

in areas where simulated water table lies in glacial deposits. 

4along losing reaches of Streams W and NW. 

5where bedrock crops out at east side of Mirror Lake.

-600

Figure 23. Three-dimensional 

extent of Mirror Lake groundwater basin in simulation D. 

See figure 22 for explanation.
46

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

of V (y) is the standard deviation, s , for the Ith

precipitationl I streams Mirror Lake

«,

\

R2

i ' \

i

i

1 t , 1

D2

1 1 \

D3 D

glacial deposits lacustrine 

mud

D ^L 1 A T f? A^r If

°1 f | D2 | D3

bedrock

I .

• Mirror Lake ground-water basin I

Figure 24. Components of ground-water budget 

for the Mirror Lake ground-water basin.

the upper hillsides are unsaturated, flow from glacial 

deposits to bedrock is relatively small. By contrast, 

flow from bedrock to glacial deposits (52) occurs 

primarily beneath lower hillsides and valleys, where 

the vertical flow component is predominantly upward. 

Flow from bedrock to the lacustrine mud or to the 

outcrop on the east side of Mirror Lake is minimal. 

Instead, nearly all flow from the bedrock is into the 

overlying glacial deposits. Therefore, flow from 

bedrock to glacial deposits is relatively large.

Uncertainties in the budget components given in 

table 10 can be assessed by calculating the approximate, individual, confidence interval for the 

components. The i* budget component is denoted by 

y ., and all budget components are contained in the 

vector y. To calculate the confidence intervals, it is 

first necessary to calculate the variance-covariance 

matrix for the budget components, which is (Hill, 

1994, p. 28, eq. 12)

(13)

where V(b) is the variance-covariance matrix for the 

model parameters (see equation 5), and X is a matrix 

containing derivatives of budget components, y, with 

respect to model parameters, b, evaluated at the

XV

optimal parameters b . These derivatives are numerically calculated using a central-difference scheme. 

The square root of the /th element on the main diagonal

budget component. The approximate, individual, 

100(1-a) percent confidence interval for the / budget 

component is given by (Hill, 1994, p. 28, eq. 11)

( a 

y.± s x t\ n-p, 1 -- 

yi y V ^ 2

(14)

where t(n-p,\-a/2) is the 100(1-0/2) percentage point 

of the Student-t distribution with n-p degrees of 

freedom.

Approximate, individual, 95-percent confidence intervals for budget components in simulation A 

are shown in table 10. These confidence intervals are 

relatively narrow, suggesting that the simulated budget 

components are relatively well constrained. However, 

caveats on confidence intervals for model parameters 

apply equally to confidence intervals for budget 

components. Because flow from lacustrine mud to 

Mirror Lake (D3) and flow from bedrock to Mirror 

Lake (D4) are both minute, their confidence intervals 

are not calculated because they are probably not 

meaningful. Confidence intervals for budget 

components in simulation D have similar widths as the 

corresponding ones in simulation A.

FLOW PATHS THROUGH THE 

GROUND-WATER BASIN

Of the total recharge to the Mirror Lake groundwater basin, a portion travels through the basin along 

flow paths that remain in glacial deposits, while the 

remainder travels along flow paths that involve 

movement in bedrock. These two types of flow paths 

are illustrated in figure 25. For simplicity, all ground 

water is shown to discharge into Mirror Lake; flow to 

the inlet streams is not illustrated. The stippled zone 

shows flow paths that remain in glacial deposits. 

Outside the stippled zone, flow paths involve 

movement in bedrock. The distinction between the 

two types of flow paths might be of importance to the 

chemical budget for Mirror Lake, because ground 

water that flows exclusively in glacial deposits might 

contain different amounts of dissolved chemicals than 

ground water that contacts bedrock. Note, however, 

that the travel distance in bedrock varies from 

relatively short (for example, flow path P-P') to 

relatively long (for example, flow path Q-Q')-
SUMMARY 47

Not to scale

Figure 25. Flow paths through the Mirror Lake 

ground-water basin.

According to figure 24, the amount of groundwater flow that involves movement in bedrock is 

simply the total inflow to bedrock (R± + B{), which 

also equals total outflow from bedrock (B^ + #3 + D4). 

Calculations using budget components in table 10 

indicate that 181,000 m3/year (simulation D) to 

190,000 m3/year (simulation A) of ground-water flow 

involves movement in bedrock. This range is approximately 60 percent of total recharge to the groundwater basin. The amount of ground-water flow that 

remains in glacial deposits is the total recharge minus 

flow that involves movement in bedrock (that is, 

RT -Rl~ BI). Calculations using budget components 

in table 10 indicate 116,000 m3/year (simulation D) to

o

117,000 m /year (simulation A) of ground-water flow 

travels exclusively in glacial deposits. This range is 

approximately 40 percent of total recharge to the 

ground-water basin.

These calculations indicate that simulations A 

and D yield similar amounts of ground-water flow that 

involves movement in the bedrock and that remains 

exclusively in the glacial deposits. Simulation D is 

expected to be the more realistic of the two simulations because calibration D matches the new hydraulic 

head observation at FS7 in the upper altitudes of the 

Mirror Lake area more closely than calibration A.

However, the significant difference in water-table 

configuration between the two simulations at upper 

altitudes has a small effect on the calculated amount of 

ground-water flow through the bedrock and the glacial 

deposits over the Mirror Lake ground-water drainage 

basin.

SUMMARY

The Mirror Lake area discussed herein covers a

r\

10-km region in the eastern part of the Hubbard 

Brook valley in central New Hampshire. This area 

includes Mirror Lake, the three streams that flow into 

the lake, Leeman's Brook, Paradise Brook, and parts 

of Hubbard Brook and the Pemigewasset River. The 

topography of this area is characterized by steep 

hillsides and relatively flat valleys. Major hydrogeologic units include glacial deposits composed of till 

containing pockets of sand and gravel, and fractured 

crystalline bedrock composed of schist intruded by 

granite, pegmatite, and lamprophyre. A layer of 

lacustrine mud sediments covers much of the bottom 

of Mirror Lake.

Ground-water flow in the Mirror Lake area is 

typical of mountain-and-valley terrain of New 

England uplands. On hillsides, water from precipitation and snowmelt infiltrates to the water table, flows 

downslope through the saturated glacial deposits and 

fractured bedrock, and discharges to streams and to 

Mirror Lake. The water table generally is within 20 m 

of land surface, and in many locations, less than 10m 

below land surface. At a well drilled on a ridge near 

the highest point of the study area, the water table is 

about 15 m below land surface. The shallow water 

table on the ridge suggests that the water table lies 

close to land surface over the entire Mirror Lake area.

A computer model was developed to simulate 

the three-dimensional, steady-state (long-term 

average) flow of ground water in the Mirror Lake area. 

A streamflow routing package was included in the 

model to simulate baseflow in streams and interaction 

between streams and ground water. The model 

assumes that recharge from precipitation is areally 

uniform, riparian evapotranspiration along stream 

banks can be neglected, flow in the bedrock occurs 

primarily in the uppermost 150 m, and distribution of 

hydraulic conductivity can be represented by dividing 

the flow domain into several zones, each having
48 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

uniform properties. Because local variations in 

recharge and hydraulic conductivities are ignored, the 

simulation results characterize the generalized groundwater flow system, not local details.

The model was calibrated using nonlinear leastsquares regression to match long-term average 

hydraulic heads measured in piezometers and bedrock 

wells, and long-term average baseflow in the three 

inlet streams to Mirror Lake. Long-term average 

baseflow was assumed to be 40 percent of long-term 

average streamflow as measured in flumes. Model 

calibration yielded 26 to 28 cm/year of recharge from 

precipitation to the water table, and hydraulic conductivities of 1.7 x 10'6 to 2.7 x 10'6 m/s for glacial 

deposits, about 3 x 10~7 m/s for bedrock beneath lower 

hillsides and valleys, and about 6 x 10~8 m/s for 

bedrock beneath upper hillsides and hilltops. Analysis 

of parameter uncertainty indicates that the above 

values are well constrained, at least within the context 

of regression analysis. In the regression, several 

attributes of the ground-water flow model are assumed 

perfectly known. The value of hydraulic conductivity 

for bedrock beneath upper hillsides and hilltops was 

determined from few data, and additional hydraulichead observations are required to confirm this result. 

Model fit was not improved by introducing a 10-to1 ratio of horizontal-to-vertical anisotropy in hydraulic 

conductivity of glacial deposits, or by varying 

hydraulic conductivity with depth in the modeled part 

(uppermost 150 m) of the bedrock.

Model simulation indicates that, in significant 

parts of the model domain, the water table lies in 

bedrock. Outside these areas, the water table lies in 

glacial deposits. Where the water table lies in 

bedrock, glacial deposits are either unsaturated or 

absent; recharge from precipitation enters the bedrock 

either by infiltration through the overlying unsaturated 

glacial deposits, or directly from land surface where 

bedrock crops out. Unsaturated glacial deposits are 

mostly on upper hillsides and away from streams.

The calibrated model was used to delineate the 

Mirror Lake ground-water basin, which is defined as 

the volume of subsurface through which ground water 

flows from the water table to Mirror Lake or its inlet 

streams. Results indicate that Mirror Lake and its inlet 

streams drain an area of ground-water recharge that is 

about 1.5 times the area of the surface-water basin. 

The ground-water basin extends far up the hillside on 

the northwestern part of the study area. Ground water 

from this area flows at depth beneath Norris Brook and

discharges into Mirror Lake or its inlet streams. As a 

result, the Mirror Lake ground-water basin extends 

beneath the adjacent ground-water basin that drains 

into Norris Brook.

Ground-water budget calculations indicate that 

approximately 300,000 m3/year of precipitation 

recharges the Mirror Lake ground-water basin. About 

half the recharge enters the basin in areas where the 

simulated water table lies in glacial deposits; the other 

half enters the basin in areas where the simulated 

water table lies in bedrock. Flow from glacial 

deposits to bedrock occurs primarily under upper 

hillsides, where the vertical flow component is 

predominantly downward. Flow from bedrock to 

glacial deposits occurs primarily under lower hillsides 

and valleys, where the vertical flow component is 

predominantly upward. About 40 percent of recharge 

(120,000 m3/year) travels through the basin along flow 

paths that remain exclusively in glacial deposits; about 

60 percent (180,000 m3/year) travels along flow paths 

that involve movement in bedrock. Under steady 

state, total recharge to the ground-water basin equals 

total discharge from the basin. Of the total discharge, 

slightly more than half (170,000 m /year) flows into 

the three Mirror Lake inlet streams; slightly less than 

half (130,000 m3/year) flows directly into Mirror 

Lake.

REFERENCES CITED

Boelter, D.H., 1972, Methods for analysing the hydrological 

characteristics of organic soils in marsh-ridden areas, 

in Hydrology of marsh-ridden areas, Proceedings of 

the Minsk Symposium, June, 1972: Paris, The Unesco 

Press, p. 161-169.

Bouwer, H., and Rice, R.C., 1976, A slug test for

determining hydraulic conductivity of unconfined 

aquifers with completely or partially penetrating wells: 

Water Resources Research, v. 12, no. 3, p. 423—428.

Cooley, R.L., and Naff, R.L., 1990, Regression modeling of 

ground-water flow: U.S. Geological Survey 

Techniques of Water-Resources Investigations, 

Book 3, Chapter B4, 232 p.

Cooper, H.H., Jr., Bredehoeft, J.D., and Papadopulos, I.S., 

1967, Response of a finite-diameter well to an instantaneous charge of water: Water Resources Research, v. 3, 

no. 1, p. 263-269.

Davis, M.B., and Ford, M.S.J., 1982, Sediment focusing in 

Mirror Lake, New Hampshire: Limnology and 

Oceanography, v. 27, no. 1, p. 137-150.
REFERENCES CITED 49

Demir, Zafer, 1992, Some contributions to Hvorslev tests, 

slug tests, and drill-stem tests: Berkeley, University of 

California, M.S. thesis, 107 p.

Federer, C.A., Flynn, L.D., Martin, C.W., Hornbeck, J.W., 

and Pierce, R.S., 1990, Thirty years of hydrometeorologic data at the Hubbard Brook Experimental Forest, 

New Hampshire: U.S. Department of Agriculture, 

Forest Service, Northeastern Forest Experiment 

Station, General Technical Report NE-141, 44 p.

Frape S.K., and Patterson, R.J., 1981, Chemistry of interstitial water and bottom sediments as indicators of 

seepage patterns in Perch Lake, Chalk River, Ontario: 

Limnology and Oceanography, v. 26, p. 500-517.

Freeze, R.A., and Cherry, J.A., 1979, Groundwater: 

Englewood Cliffs, N.J., Prentice-Hall, 604 p.

Haeni, P.P., Lane, J.W., Jr., and Lieblich, D.A., 1993, Use of 

surface-geophysical and borehole-radar methods to 

detect fractures in crystalline rocks, Mirror Lake area, 

Grafton County, New Hampshire, in Banks, Sheila, 

and Banks, David, eds., IAH memoires of the 24th 

congress, hydrogeology of hard rocks, As, Norway, 

June 28 to July 2, 1993: International Association of 

Hydrogeologists, p. 577-587.

Harbaugh, A.W., 1990, A computer program for calculating 

subregional water budgets using results from the U.S. 

Geological Survey modular three-dimensional finitedifference ground-water flow model: U.S. Geological 

Survey Open-File Report 90-392,46 p.

Hardin, E.L., Cheng, C.H., Paillet, PL., and Mendelson, 

J.D., 1987, Fracture characterization by means of 

attenuation and generation of tube waves in fractured 

crystalline rock at Mirror Lake, New Hampshire: 

Journal of Geophysical Research, v. 92, no. B8, 

p. 7989-8006.

Harte, P.A., 1992, Regional ground-water flow in crystalline 

bedrock and interaction with glacial drift in the New 

England uplands: Durham, University of New 

Hampshire, M.S. thesis, 141 p.

Harte, P.A., 1994, Comparison of vertical discretization 

techniques in finite-difference models of groundwater flow: Example from a hypothetical New 

England setting: U.S. Geological Survey Open-File 

Report 94-343, 25 p.

Hill, M.C., 1990, Preconditioned conjugate-gradient 2

(PCG2), a computer program for solving ground-water 

flow equations: U.S. Geological Survey WaterResources Investigations Report 90-4048,43 p.

Hill, M.C., 1992, A computer program (MODFLOWP) for 

estimating parameters of a transient, threedimensional, ground-water flow model using nonlinear 

regression: U.S. Geological Survey Open-File Report 

91-484, 358 p.

Hill, M.C., 1994, Five computer programs for testing 

weighted residuals and calculating linear confidence

and prediction intervals on results from the groundwater parameter-estimation computer program 

MODFLOWP: U.S. Geological Survey Open-File 

Report 93-481,81 p.

Hsieh, P.A., Perkins, R.L., and Rosenberry, D.O., 1996, 

Field instrumentation for multilevel monitoring of 

hydraulic head in fractured bedrock at the Mirror Lake 

site, Grafton County, New Hampshire, in Morganwalp, 

D.W., and Aronson, D.A., eds., U.S. Geological 

Survey Toxic Substances Hydrology 

Program—Proceedings of the Technical Meeting, 

Colorado Springs, Colorado, September 20-24, 1993: 

U.S. Geological Survey Water-Resources Investigations Report 94-4015, v. 1, p. 137-140.

Hsieh, P.A., and Shapiro, A.M., 1994, Analyzing cross-hole 

hydraulic tests in fractured crystalline rock: An 

example from the Mirror Lake site, New Hampshire 

[abs.]: Eos, Transactions, American Geophysical 

Union, v. 75, no. 16, supplement, 1994 Spring 

Meeting, p. 153.

Hsieh, P.A., Shapiro, A.M., Barton, C.C., Haeni, P.P., 

Johnson, C.D., Martin, C.W, Paillet, F.L., Winter, 

T.C., and Wright, D.L., 1993, Methods of characterizing fluid movement and chemical transport in 

fractured rocks, in Cheney, J.T., and Hepburn, J.C., 

eds., Field trip guidebook for the northeastern United 

States: 1993 Boston GSA: Amherst, University of 

Massachusetts, Department of Geology and 

Geography, Contribution No. 67, p. R1-R30.

Institute of Hydrology, 1980a, Low flow studies: Wallingford, United Kingdom, Institute of Hydrology Report, 

no. 1.

Institute of Hydrology, 1980b, Low flow studies: Wallingford, United Kingdom, Institute of Hydrology Report, 

no. 3.

Lane, J.W, Jr., Haeni, P.P., and Watson, W.M., 1995, Use of 

a square-array direct-current resistivity method to 

detect fractures in crystalline bedrock in New 

Hampshire: Ground Water, v. 33, no. 3, p. 476-485.

Likens, G.E., Bormann, PH., Pierce, R.S., Eaton, J.S., and 

Johnson, N.M., 1977, Biogeochemistry of a forested 

ecosystem: New York, Springer-Verlag, 146 p.

Likens, G.E., ed., 1985, An ecosystem approach to aquatic 

ecology, Mirror Lake and its environment: New York, 

Springer-Verlag, 516 p.

Likens, PC., 1991, Publications of the Hubbard Brook 

Ecosystem Study: Millbrook, N.Y., Institute of 

Ecosystem Studies, 109 p.

Lyons, J.B., Bothner, W.A., Moench, R.H., and Thompson, 

J.B., Jr., eds., 1986, Interim geologic map of New 

Hampshire: Concord, New Hampshire Department of 

Resources and Economic Development, Open-File 

Report 86-1, 1 sheet, scale 1:250,000.
50 NUMERICAL SIMULATION OF GROUND-WATER FLOW THROUGH GLACIAL DEPOSITS AND CRYSTALLINE 

BEDROCK IN THE MIRROR LAKE AREA, GRAFTON COUNTY, NEW HAMPSHIRE

Mau, D.P., and Winter, T.C., 1997, Estimating ground-water 

recharge from streamflow hydrographs for a small 

watershed in a temperate humid climate, New 

Hampshire, USA: Ground Water, v. 35 [in press].

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, Chapter Al, 586 p.

McHone, J.G., 1984, Mesozoic igneous rocks of northern 

New England and adjacent Quebec: Summary, description map, and bibliography of data sources: Geological 

Society of America, Map and Chart Series, MC-49.

Moeller, R.E., 1978, The hydrophytes of Mirror Lake: A 

study of vegetational structure and seasonal biomass 

dynamics: Ithica, N.Y., Cornell University, Ph.D. 

dissertation, 212 p.

Morganwalp, D.W., and Aronson, D.A., eds., 1996, U.S. 

Geological Survey Toxic Substances Hydrology 

Program—Proceedings of the Technical Meeting, 

Colorado Springs, Colorado, September 20-24, 1993: 

U.S. Geological Survey Water-Resources Investigations Report 94-4015, 1113 p.

Paillet, F.L., 1985, Geophysical well log data for study of 

water flow in fractures in bedrock near Mirror lake, 

West Thornton, New Hampshire: U.S. Geological 

Survey Open-File Report 85-340, 27p.

Paillet, F.L., Hess, A.E., Cheng, C.H., and Hardin, E.L., 

1987, Characterization of fracture permeability with 

high-resolution vertical flow measurements during 

borehole pumping: Ground Water, v. 25, no. 1, 

p. 28-40.

Paillet, F.L., and Kapucu, K., 1989, Fracture characterization and fracture-permeability estimates from geophysical logs in the Mirror Lake watershed, New 

Hampshire: U.S. Geological Survey Water-Resources 

Investigations Report 89-4058,49 p.

Pollack, D.W., 1994, User's guide for

MODPATH/MODPATH-PLOT, version 3: A particle 

tracking post-processing package for MODFLOW, the 

U.S. Geological Survey finite-difference ground-water 

flow model: U.S. Geological Survey Open-File Report 

94-0464, 249 p.

Prudic, D.E., 1989, Documentation of a computer program 

to simulate stream-aquifer relations using a modular, 

finite-difference, ground-water flow model: U.S. 

Geological Survey Open-File Report 88-729, 113 p.

Rosenberry, D.O., and Winter, T.C., 1993, The significance 

of fracture flow to the water balance of a lake situated 

in fractured crystalline rock terrane, in Banks, Sheila, 

and Banks, David, eds., IAH Memoires of the 24th 

Congress, Hydrogeology of Hard Rocks, As, Norway,

June 28 to July 2, 1993: International Association of 

Hydrogeologists, p. 967-977.

Rutledge, A.T., 1993, Computer programs for describing 

the recession of ground-water discharge and for 

estimating mean ground-water recharge and discharge 

from streamflow records: U.S. Geological Survey 

Water-Resources Investigations Report 93-4121,45 p.

Rutledge, AT., and Daniel, C.C., III, 1994, Testing an 

automated method to estimate ground-water recharge 

from streamflow records: Ground Water, v. 32, no. 2, p. 

180-189.

Shapiro, A.M., and Hsieh, P.A., 1991, Research in

fractured-rock hydrogeology: Characterizing fluid 

movement and chemical transport in fractured rock at 

the Mirror Lake drainage basin, New Hampshire, in 

Mallard, G.E., and Aronson, D.A., eds., U.S. 

Geological Survey Toxic Substances Hydrology 

Program—Proceedings of the technical meeting, 

Monterey, California, March 11-15, 1991: 

U.S. Geological Survey Water-Resources Investigations Report 91-4034, p. 155-161.

Shattuck, PC., 1991, Shallow water-table response to 

precipitation and evapotranspiration in an ephemeral 

stream valley, Woodstock, New Hampshire: Durham, 

University of New Hampshire, M.S. thesis, 113 p.

U.S. Department of Agriculture Forest Service, 1991, The 

Hubbard Brook Ecosystem Study: Site Description and 

Research Activities: Northeastern Forest Experimental Station NE-INF-96-91, 61 p.

Wilson, Alicia, 1991, Distribution of hydraulic conductivity 

in the glacial drift at Hubbard Brook Experimental 

Forest, West Thornton, New Hampshire: Hanover, 

N.H., Dartmouth College, Department of Earth 

Sciences, Senior honors thesis, 56 p.

Winter, T.C., 1984, Geohydrologic setting of Mirror Lake, 

West Thornton, New Hampshire: U.S. Geological 

Survey Water-Resources Investigations 

Report 84-4266, 61 p.

Winter, T.C., 1985, Physiographic setting and geologic 

origin of Mirror Lake, in Likens, G.E., ed., An 

ecosystem approach to aquatic ecology, Mirror Lake 

and Its Environment: New York, Springer-Verlag, 

p. 40-53.

Winter, T.C., Eaton, J.S., and Likens, G.E., 1989, Evaluation of inflow to Mirror Lake, New Hampshire: Water 

Resources Bulletin, v. 25, no. 5, p. 991-1008.