U.S. Department of the Interior

U.S. Geological Survey

Simulation of the Recharge Area for

Frederick Springs, Dane County,

Wisconsin

Water-Resources Investigations Report 00–4172

Prepared in cooperation with the

City of Middleton

Wisconsin Department of Natural Resources

Probability

of capture 

1

240 Row

1 160 Column

43°02'30"

89°30' 89°35'

47°07'30"

0.00

1.00

Pheasant Branch stream network

Other stream

Surface-water divide

Lake

Mendota

18

151

12

14

EXPLANATION

012

012

MILES

KILOMETERS
Simulation of the Recharge Area for 

Frederick Springs, Dane County, Wisconsin

By R.J. Hunt and J.J. Steuer

U.S. GEOLOGICAL SURVEY

Water-Resources Investigations Report 00–4172

Prepared in cooperation with the

City of Middleton

Wisconsin Department of Natural Resources

Middleton, Wisconsin

2000
CONTENTS III

CONTENTS

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

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

Purpose and scope....................................................................................................................................................... 2

Physical setting of study area...................................................................................................................................... 2

Acknowledgments....................................................................................................................................................... 5

Method of study .................................................................................................................................................................... 5

Construction of telescopic mesh refinement model.................................................................................................... 5

Parameter estimation................................................................................................................................................... 8

Stochastic Monte Carlo analyses ................................................................................................................................ 9

Geochemical investigation.......................................................................................................................................... 11

Results of model simulations and geochemical investigation .............................................................................................. 11

Telescopic mesh refinement........................................................................................................................................ 11

Parameter estimation................................................................................................................................................... 13

Stochastic Monte Carlo analyses ................................................................................................................................ 17

Geochemical investigation.......................................................................................................................................... 19

Summary and conclusions .................................................................................................................................................... 31

References............................................................................................................................................................................. 32

FIGURES

1–3. Maps showing:

1. Relation of the telescopic mesh refinement model to the Dane County Regional Model ................................. 3

2. Surface water features and high capacity wells used in the telescopic mesh refinement model........................ 4

3. Locations of hydrologic response units used in the Pheasant Branch surface water modeling ......................... 7

4. Plot of parameter sensitivity from the UCODE parameter estimation..................................................................... 10

5. Plots of residuals for model runs .............................................................................................................................. 14

6. Map showing area where the connection between layers 2 and 3 was modified in the telescopic mesh 

refinement model .................................................................................................................................................. 16

7–9. Results of stochastic MODFLOW runs showing:

7. Mean absolute error for 200 realizations............................................................................................................ 18

8. Conditioned stochastic MODFLOW runs for the upper bedrock aquifer and the lower bedrock aquifer ......... 20

9. Conditioned stochastic MODPATH runs for the upper bedrock aquifer and the lower bedrock aquifer........... 22

10. Map showing location of the Springs and geochemical sampling sites ................................................................... 25

11–15. Graphs showing:

11. Average specific conductance measurements and standard deviations for 1998–99 ........................................ 26

12. Average calcium concentrations and standard deviations for 1998–99............................................................. 26

13. Average 87Sr/86Sr ratios for March 1998 and September 1998 sampling ......................................................... 27

14. Average δ18O for 1998–99 ................................................................................................................................ 27

15. Average 87Sr/86Sr ratios versus average calcium concentrations for 1998–99 ................................................. 28

16. Results of particle tracking over the bedrock map ................................................................................................... 29

17–19. Graphs showing:

17. Average nitrate + nitrite concentrations and standard deviations for 1998–99 ................................................. 30

18. Average nitrate + nitrite concentrations versus average chloride concentrations for the Frederick Spring 

sampling, 1998–99 .............................................................................................................. ........................... 30

19. Average total filterable phosphorus concentrations and standard deviations for 1998–99 ............................... 31

TABLES

1. Parameters used in the Dane County Regional Model and telescopic mesh refinement model............................... 6

2. Results of the Pheasant Branch UCODE runs using the telescopic mesh refinement model and the 

Dane County Regional Model .............................................................................................................................. 12

3. Results of geochemical sampling, 1998–99 ............................................................................................................. 24
IV CONVERSION FACTORS, VERTICAL DATUM, AND ABBREVIATED WATER-QUALITY UNITS

CONVERSION FACTORS, VERTICAL DATUM, AND ABBREVIATED WATER-QUALITY UNITS

Multiply By To Obtain

Length

inch (in.) 25.4 millimeter

foot (ft) 0.3048 meter

mile (mi) 1.609 kilometer

Area

acre (A) 0.4047 hectare

square foot (ft2) 0.09290 square meter

square mile (mi2

) 2.590 square kilometer

Volume

cubic foot (ft3) 7.4805 gallon

Hydraulic conductivity*

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

Sea level: In this report, “sea level” refers to the National Geodetic Vertical Datum of 1929 (NGVD 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.

*Hydraulic conductivity: The standard unit for hydraulic conductivity is cubic foot per day per square foot of aquifer cross-sectional area 

(ft3/d)/ft2. In this report, the mathematically reduced form, feet per day (ft/d), is used for convenience.

Other abbreviations:

ft3/s cubic feet per second

in/yr inches per year

ft/d feet per day

Kh horizontal hydraulic conductivity

Kv vertical hydraulic conductivity

The stratigraphic nomenclature used in this report is that of the Wisconsin Geological and Natural History Survey and does not necessarily 

follow usage of the U.S. Geological Survey.
U.S. DEPARTMENT OF THE INTERIOR

BRUCE BABBITT, Secretary

U.S. GEOLOGICAL SURVEY

Charles G. Groat, Director

For additional information write to: Copies of this report can be purchased from:

District Chief U.S. Geological Survey 

U.S. Geological Survey Branch of Information Services

8505 Research Way Box 25286

Middleton, WI 53562-3586 Denver, CO 80225-0286

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

by the U.S. Government.
Abstract 1

Abstract

The Pheasant Branch watershed in Dane 

County is expected to undergo development. There 

are concerns that this development will adversely 

affect water resources, including Frederick 

Springs, a large spring complex in the watershed. 

The spring’s recharge area was delineated using a 

telescopic mesh refinement (TMR) model constructed from an existing regional-scale groundwater flow model, and further refined by adding 

nearby surface-water features, a refined recharge 

array based on a surface-water model, and increasing the vertical leakage between the deep aquifers. 

This TMR model was formally optimized using the 

parameter estimation code UCODE. The results of 

optimization demonstrated that the best fit to measured heads and fluxes was obtained by using a 

horizontal hydraulic conductivity two times that of 

the original regional model for layer 2 and 80 percent smaller for layer 3. This range of parameter 

values was formally considered using a stochastic 

Monte Carlo approach. 

Two-hundred model runs used uniformly 

distributed, randomly sampled, horizontal hydraulic conductivity values within the range given by 

the TMR optimized values and the previously constructed regional model. A probability distribution 

of particles captured by the spring, or a “probabilistic capture zone” was calculated from the realistic 

Monte Carlo results (136 runs of 200). In addition 

to portions of the local surface watershed, the capture zone encompassed distant areas in the North 

Fork of the Pheasant Branch watershed and areas 

entirely outside of the Pheasant Branch—demonstrating that the ground-watershed and surface 

watershed do not coincide. 

Analysis of samples from the springs and a 

nearby municipal well identified large contrasts in 

chemistry, even for springs within 50 feet of one 

another. The differences were stable over time, 

were present in both ion and isotope analyses, and 

showed a distinct gradation from high nitrate, high 

calcium, Ordovician-carbonate dominated water in 

western spring vents to low nitrate, lower calcium, 

Cambrian-sandstone influenced water in eastern 

spring vents. The difference in chemistry was 

attributed to distinctive bedrock geology as demonstrated by overlaying the 50 percent probability 

capture zone over a bedrock geology map for the 

area. This finding gives additional confidence to 

the capture zone calculated by the ground-water 

flow model.

INTRODUCTION

As the City of Middleton and its surroundings continue to develop, the Pheasant Branch watershed is 

expected to undergo significant urbanization. The 

watershed encompasses mixed land uses and includes 

Frederick Springs, a large spring complex near the 

downstream end of the basin. For the downstream city 

of Middleton, headwater urbanization can mean 

increased flood peaks, increased water volume and 

increased pollutant loads. More subtly, the effects of 

urbanization also may reduce ground-water recharge 

and adversely affect down-gradient ecosystems such as 

Pheasant Branch marsh. The effects of stormwater runoff and wetland loss on reduced ground-water recharge 

are complex because the surface-water system is coupled to the underlying ground-water system as in most 

parts of Wisconsin. In many cases, the movement of 

water from one system to the other varies seasonally or 

daily, depending on transient events. Therefore, it is difficult to reliably predict the effects of urbanization on 

stream baseflow and spring flows a priori. Moreover, 

mitigating any adverse effects after development has 

taken place can be expensive and administratively difficult. Overlying these concerns are issues such as land 

owners’ rights—both of those developing their land and 

those whose land is affected by this development—the 

rights of the public, and stewardship of the resource. 

With these often-contradictory goals, a scientific basis 

for assessing effects and effectiveness of mitigation 

measures is important for effective decision making. 

This scientific basis was the goal of a study completed 

Simulation of the Recharge Area for Frederick Springs, 

Dane County, Wisconsin

By R.J. Hunt and J.J. Steuer
2 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

by the U.S. Geological Survey in cooperation with the 

City of Middleton and the Wisconsin Department of 

Natural Resources.

Others have investigated portions of the Pheasant 

Branch watershed. Previous work either was included 

in a larger regional study or focused solely on the stream 

itself. Krohelski and others (2000) included the basin in 

their larger ground-water flow model for Dane County; 

however, the system was coarsely gridded and did not 

include many of the hydrologic features that are locally 

important (such as Frederick Spring, Dorn Creek). 

Grant and Goddard (1980) investigated channel erosion 

and sediment transport in the stream. Krug and Goddard 

(1986) evaluated the effects of urbanization on the 

stream system. Selbig (1996) characterized the Pheasant Branch Marsh and springs as part of a UW-Madison 

wetlands ecology course. A companion volume to this 

report describing the surface-water modeling work and 

effects of predicted stresses on the hydrologic systems 

of the Pheasant Branch watershed is in preparation.

While many hydrologic studies focus only on one 

component of the ground water-surface water continuum, the overall project includes all elements of the 

hydrologic cycle including rainfall, snowmelt, evapotranspiration, interflow, streamflow, baseflow, and 

ground-water flow. The entire hydrologic system is 

characterized quantitatively; output from surface-water 

modeling (recharge) is coupled to the ground-water 

model input. This allows more realistic scenarios (that 

is, urbanization affects surface-water storm flows and 

ground-water recharge) and allows an additional check 

for reasonableness. While both surface-water and 

ground-water modeling are being performed as part of 

the overall project, this report will focus on the groundwater modeling aspects and assessment of recharge 

rates derived from the surface-water model. 

Purpose and Scope

The purpose of the overall study was to provide a 

scientific basis for evaluating changes to the water 

resources of the Pheasant Branch watershed as the 

hydrologic system responds to changes in land use. Of 

special interest is the Frederick Springs system because 

it is an important water resource in the Pheasant Branch 

watershed and an essential source of water for a rare 

wild rice community in the Pheasant Branch Marsh. 

Identifying the source waters for the spring is the first 

step in ensuring its protection. The purpose of this 

report is to describe ground-water flow modeling and 

geochemical information used to identify the source of 

water that sustains the spring complex.

Physical Setting of Study Area

The Pheasant Branch watershed (fig. 1) consists of 

24 mi2 located on the edge of the Driftless Area in Dane 

County. The geology of the Pheasant Branch area was 

derived from the hydrogeologic framework (Bradbury 

and others, 1999) and existing regional model of the 

area (Krohelski and others, 2000) and, from base 

upward, consists of: 

(1) an impermeable Precambrian basement, 

(2) an extensive lower bedrock aquifer that consists of a Cambrian sandstone aquifer (Mt. 

Simon and Eau Claire sandstone) that ranges 

between 400 and 700 ft and averages 525 ft 

thick in the Pheasant Branch area, 

(3) a shaly confining unit (the Eau Claire Shale) 

that is absent near Lake Mendota but is as thick 

as 40 ft; the average unit thickness in the 

Pheasant Branch area is 20 ft,

(4) a upper bedrock aquifer consisting primarily of 

Cambrian sandstones and Ordovician carbonates that are absent under Lake Mendota but 

can be as thick as 625 ft in western Dane 

County; the average thickness is 320 ft in the 

Pheasant Branch area, and 

(5) a thin overlying unconsolidated sedimentary 

unit that is Quaternary in age. 

The shale confining unit is present in the western 

portions of the county, but is absent under the Yahara 

lakes and in areas in eastern Dane County. The water 

table is commonly found in the Upper Paleozoic bedrock layer in the Pheasant Branch area; high capacity 

wells often have boreholes that are open to both the 

lower bedrock aquifer and the upper bedrock aquifer.

The watershed is composed of a south fork, a north 

fork and a lower system that flows into the Pheasant 

Branch Marsh (fig. 2). The south fork is ephemeral and 

is not included in the ground-water model. At the 

marsh, flow from Pheasant Branch Creek combines 

with flows from the spring complex (Springs) and other 

ground water discharged to the marsh (to the stream 

channel as well as other springs and minor tributaries); 

this combined flow ultimately discharges into Lake 

Mendota. During present conditions, these three components (stream flow, spring flow, and ground-water 
INTRODUCTION 3

Figure 1. Relation of the telescopic mesh refinement (TMR) model to the Dane County Regional Model (DCRM). (Grid is used 

to specify perimeter boundaries for the inset TMR model. Because the spring is a regional feature, the DCRM was also used to 

coarsely define the capture zone of the spring (shown near the middle of the TMR model) and to test the efficacy of possible 

changes for the TMR model. Hydrologic features used in the DCRM and the location of the Pheasant Branch watershed are 

also shown.)

Spring capture

zone (layer 2)

43°20'

89°30' 89°

42°50'

0 5 10 KILOMETERS

0 5 10 MILES

spring

Location of study area 

in Wisconsin

151 18

151

12

12 18

14

94

39

39

94

90

90

Pheasant Branch

watershed

TMR

model
4 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

Brewery Creek 

Sugar River

constant head 

perimeter boundary

Springs

1

240 Row

1 Column 160

South Fork

North Fork

Pheasant 

Branch Marsh 

Lower

Pheasant 

Branch

Six Mile Creek

Dorn Creek

Black Earth Creek

43°02'30"

18

151

12

14

Littoral zone

Deep lake sediments

Stream nodes

EXPLANATION

Ephemeral streams

Surface water divide

High capacity well

012

012

MILES

KILOMETERS

89°30' 89°35'

47°07'30"

Lake

Mendota

Pheasant

Branch Creek

Figure 2. Surface water features and high capacity wells used in the TMR model. (Note smaller grid allows more accurate 

representation of surface water feature geometry and extent (for example, Pheasant Branch Creek). Surface water features 

important for the Pheasant Branch watershed are added to those of the DCRM (for example, Brewery Creek, Dorn Creek)).
METHOD OF STUDY 5

discharge to the marsh) are roughly equal (around 2 ft3/s 

each) during baseflow conditions (Selbig, 1996).

The drainage network of the watershed has been 

significantly modified over the last 150 years. Prior to 

the mid-1880’s, the Pheasant Branch watershed 

upstream of Highway 12 drained into a large wetland 

that occupied flat-lying land that surrounds the present 

confluence of the North and South Forks (Maher, 1999). 

The watershed was closed in most years, but in 

extremely wet years may have spilled into the Black 

Earth Creek watershed to the west. In the mid-1880’s 

the wetland was drained to Lake Mendota. Most of the 

present-day channels in the Pheasant Branch watershed 

formed after land in the watershed was converted to 

agricultural uses. The channel that extends from Highway 12 to the Pheasant Branch Marsh has a high gradient (90 ft drop over 2 mi) resulting in high rates of 

erosion that threaten bridges and sewer lines. The City 

of Middleton has spent approximately $500,000 in the 

last 30 years in an attempt to protect these structures 

from erosion. 

The Pheasant Branch Marsh is a diverse wetland 

system that contains stands of native herbaceous and 

shrub-scrub vegetation near the Springs, reed canary 

grass and stinging nettles near the Pheasant Branch outlet, and wild rice and cattails near the confluence of the 

Spring and Pheasant Branch Creek channel (Selbig, 

1996). Originally the Pheasant Branch channel flowed 

north into the Pheasant Branch Marsh and discharged to 

the northwest part of the wetland, but during the 1950’s 

the channel changed course and followed the southern 

boundary of the marsh (Grant and Goddard, 1980) as a 

result of dumping near the old outlet (D’Onofrio, 1973). 

By 1971, the stream was re-routed to the northwest outlet by the City of Middleton in an attempt to trap sediments and reduce the overall stream gradient. However, 

the stream system still conveys large amounts of sediment and associated phosphorus to the Lake Mendota 

system. For example, the Pheasant Branch system had 

the highest sediment load per unit area for all rural 

streams measured in Dane County (Lathrop and 

Johnson, 1979). Increased stormwater flows, a consequence of future development, are expected to exacerbate erosion in the stream channel and sediment 

transport; these issues have become a topic of concern 

for the citizens of Middleton (North Fork Pheasant 

Branch Watershed Committee, 1999).

The lower Pheasant Branch connecting the Pheasant Branch Springs to Lake Mendota is distinct from the 

rest of the system in a number of ways. First, a channel 

connects the large spring complex to Pheasant Branch 

Creek and has probably existed in this location for thousands of years. Secondly, the stream in this area has a 

flat gradient. Finally, whereas the rest of the stream is 

characterized by relative low ground-water discharge 

and high stormwater flows, the springs and associated 

channel have relatively high ground-water discharge 

components.

Acknowledgments

Special thanks are extended to Teresa Mansor for 

her assistance in collecting and analyzing chemical 

data. Daniel Feinstein and James Krohelski assisted 

with regional model simulations. Tom Bullen is thanked 

for providing strontium isotope analyses. The City of 

Middleton Public Works Department and the Friends of 

Pheasant Branch are acknowledged for their assistance 

throughout the project. 

METHOD OF STUDY

Construction of Telescopic Mesh Refinement 

Model

The model approach consists of modifying an 

inset model extracted from a recently developed 

MODFLOW (McDonald and Harbaugh, 1988) model 

(fig. 1) for Dane County (Krohelski and others, 2000), 

hereafter called the Dane County Regional Model 

(DCRM). The DCRM grid was used to assign constant 

head boundary conditions along the perimeter of the 

inset model (fig. 2) and the inset model grid was refined 

using a telescopic mesh refinement (TMR) approach 

following the methodology of Ward and others (1987). 

The TMR routine was performed using Groundwater 

Vistas (Rumbaugh and Rumbaugh, 1996). The original 

grid spacing was 1,312.4 ft on a side; the refined grid is 

one-fourth the spacing, or 328.1 ft on a side. Thus, every 

DCRM model node is represented by 16 nodes in the 

TMR model. As a result of the smaller grid spacing, surface-water features and the hydraulic head distribution 

were represented more accurately in the TMR model. 

Therefore, stream-aquifer interactions and three-dimensional effects near the streams were simulated more 

accurately. A second refinement was the conversion of 

all streams in the TMR model domain from the MODFLOW River Package (McDonald and Harbaugh, 1988) 

to the more sophisticated Stream Routing Package (Prudic, 1989). The latter allows accounting of streamflow 

in the streams, and limits the amount of water a stream 
6 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

can lose to the aquifer to the amount of water captured 

upstream. The simpler River Package was used to simulate Lake Mendota, however, as neither stream routing 

nor the ability to control infiltration quantities were 

required for this feature. The final modification to the 

DCRM model construction was the refinement of the 

representation of the existing streams and the addition 

of surface-water features not previously modeled (such 

as the Springs and Dorn Creek northeast of the Springs 

—fig. 2) that were not critical for the regional model 

calibration but were important for the TMR model calibration. The additional surface-water features were simulated using the Stream Routing Package. Added 

stream nodes used an assumed stream bottom elevation 

equal to the 5 ft below stream stage and a sediment 

thickness equal to 1 ft. Stream width was set equal to 

5 ft in headwater reaches; width was assigned according 

to stream order elsewhere, with larger widths specified 

for reaches farther from the headwaters.

In addition to changes in model construction, 

model parameters were adjusted during the TMR model 

refinement. As part of the TMR routine the properties of 

the DCRM were directly translated; a subset of these 

parameters was modified during subsequent TMR 

model runs. One of the most important was recharge 

over the Pheasant Branch watershed. Recharge rates 

determined by surface-water modeling were input to the 

ground-water flow model via the MODFLOW recharge 

array. The surface-water model hydrologic response 

units (HRUs) were grouped such that hydrologic characteristics of an area (soil infiltration capacity, evapotranspiration, rainfall-runoff properties) were averaged 

into one representative value over the area; this grouping resulted in 21 recharge zones over the Pheasant 

Branch Watershed (fig. 3). Recharge rates ranged from 

a high of 9.5 in/yr to a low of 2.2 in/yr (table 1). This 

new recharge array corresponded to an area-weighted 

average recharge over the Pheasant Branch watershed 

Table 1. Parameters used in Dane County Regional Model (DCRM) and telescopic mesh refinement 

(TMR) model

Horizontal Hydraulic Conductivity (feet/day)

TMR base Head optimized Head and flux optimized

Layer 1 zone 1 0.3 0.3 0.3

zone 2 0.5 0.5 0.5

zone 3 1 1 1

zone 4 3 3 3

zone 5 5 5 5

zone 6 7 7 7

Layer 2 5 15.0 9.6

Layer 3 10 0.7 2.0

Nodal Recharge Rate (inches/year)

Original DCRM maximum 12.4

minimum 0.0

area weighted average 4.8

TMR model maximum 9.5

minimum 2.2

area weighted average 8.0

Leakage Underneath the Springs MODFLOW VCONT (1/day)

Original DCRM for Pheasant Branch Springs

DCRM node representing larger springs area 0.00001817

TMR node directly below Springs 0.00001817

TMR model in Pheasant Branch Springs

DCRM node representing larger springs area 0.00446

TMR node directly below Springs 0.14
METHOD OF STUDY 7

Sugar River

1

240 Row

1 Column 160

Six Mile Creek

43°02'30"

18

151

12

14

EXPLANATION

Surface-water divide

Hydrologic response unit

Springs

Highway 12 gaging station

012

012

MILES

KILOMETERS

89°30' 89°35'

47°07'30"

Pheasant Branch Creek outlet

Century Avenue measuring point

Head calibration targets

Figure 3. Locations of the hydrologic response units (HRUs) used in the Pheasant Branch surface water

modeling. (Head and flux targets used in the TMR model calibration are also shown.)
8 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

of 8.0 in/yr. It should be noted that this value is significantly higher than that used for the Pheasant Branch 

watershed in the regional Dane County Model (equal to 

4.8 in/yr). The original recharge array was considered 

theoretical, however (Krohelski and others, 2000; 

Swanson, 1996), and did not simulate baseflows for 

Pheasant Branch Creek very closely (see below). 

Other significant changes to parameters included 

changes to horizontal hydraulic conductivity (Kh) of the 

bedrock aquifers (layers 2 and 3), and increase of vertical leakage between the Springs and aquifers underneath the Springs (table 1). Leakage is defined as the 

flow of water from one hydrogeologic unit to another, 

and is input to the model through the specification of the 

VCONT array (McDonald and Harbaugh, 1988). The 

VCONT is a function of the vertical hydraulic conductivity and thickness of aquifers and confining units; 

higher values of VCONT indicate higher degrees of 

leakage. This refined model was then calibrated to the 

heads and fluxes used in the DCRM and additional flux 

data collected for the project. In particular, the calibration was aimed at more accurately simulating stream 

and spring discharge measurements while not significantly degrading the head calibration. 

Although most of the discussion focuses on the 

TMR model, many of the refinements discussed here 

were input to the DCRM before the TMR perimeter 

boundaries were extracted. In addition to giving insight 

into the type of modifications needed in the TMR model 

to adequately simulate the hydrologic system, this procedure ensured that the most appropriate perimeter 

boundary was specified. That is, the Springs themselves 

are a regional feature (discharging more on a daily basis 

than any single municipal well in the county); thus the 

Springs could conceivably have measurable influence 

at a given TMR perimeter. Conceptual model refinements, albeit coarsely discretized, to the Dane County 

Regional Flow model before the model grid was refined 

helped ensure that boundary effects would be minimized. 

Parameter Estimation

The TMR model was calibrated using parameter 

estimation techniques. The use of parameter estimation, 

or “inverse” models, for calibration is a relatively new 

advancement for the science, and numerous publications describe the advantages of inverse models (such 

as, Hill, 1992; Poeter and Hill, 1997; Hill, 1998). 

Briefly, the primary benefit of a properly constructed 

inverse model is its ability to automatically calculate 

parameter values (such as hydraulic conductivity) that 

are a quantified best fit between simulated model output 

and data measured in the real world (such as head and 

stream baseflow). Other benefits are also realized, such 

as the quantification of the quality of the calibration and 

a statistically rigorous measure of the uncertainty (that 

is, confidence interval) of predictions made using the 

optimized model. In addition, parameter correlation 

(such as hydraulic conductivity and recharge) and 

parameter sensitivity can be quantified and assessed. In 

this work, the TMR model in the Pheasant Branch area 

was coupled to the inverse code UCODE (Poeter and 

Hill, 1998). 

One of the most important operations in parameter 

estimation is the selection of observations and associated weight given to these observations. There can be 

subjective aspects to the assignment, where the ultimate 

goal is obtaining an optimization evaluation that reflects 

the modeler’s judgment of goodness of fit. Some have 

proposed using only measurement error to assign observation weights; in practice, however, the measured 

value can have additional, less quantifiable, sources of 

error. Primarily, these errors relate to how representative the measurement is to the condition being modeled. 

This model is simulating steady-state conditions when 

the hydrologic system is rarely at steady state. For 

example, measured head values can vary over tens of 

feet between years, and the location and elevation of the 

wells is often imprecisely known. Yet measurement 

error of head at a given well at a given time is roughly 

0.01 ft. Therefore, the weights assigned to the head data 

are intended to reflect their representativeness, and take 

into account the fact that all head data were not collected at the same time and that there are uncertainties 

in the location and elevation of measurement. In addition, stream flows that are point measures in time may 

have an “excellent” (± 5 percent accuracy) measurement rating, but knowledge of how well this measurement represents average conditions is unknown without 

a determination of flow duration. In addition, as shown 

in figure 3 there are many more measured head observations (454) in the model domain than the number of 

measured flux observations (4). Therefore, in practice, 

observation weights are not assigned using set criteria, 

but rather in such a way that the assessment of what 

constitutes a “best fit” in the parameter estimation routine (that, is the selection of an “optimal” model) is similar to the decision a modeler might make in trial-anderror calibration.
METHOD OF STUDY 9

The TMR model described above was optimized 

using the following criteria for head and flux targets in 

the Pheasant Branch area:

(1) The average flow (Q50) at the gaging station at 

Pheasant Branch Creek at Highway 12 was the 

most highly weighted observation because of 

the availability of the best flow duration information for the study period. The weight reflects 

the long flow record (continuously monitored 

from July 1974 through the study period) used 

to determine the flow duration at the site. The 

measured Q50 of 1.8 ft3

/s for 1974–98 was 

given a coefficient of variation = 0.01 

(= 95-percent confidence interval that spans 

± 1 percent around the measured value). This 

value is more accurate than that reported by 

Holmstrom and others (1999) for the gage (fair, 

or ± 15 percent) because the weight reflects its 

much higher quality on a relative scale to the 

other data used to optimize the model. That is, 

the weight represents the fact that we are willing to trade better results in other targets to 

have the optimization routine match this target 

well.

(2) A much smaller number of discharge measurements were made at the Pheasant Branch Creek 

at Century Avenue and at the Springs locations 

than at the Highway 12 gaging station. Therefore, these targets were given an intermediate 

weight, reflecting their shorter period of record 

and fewer discharge measurements (coefficient 

of variation = 0.3 and 0.2, respectively). 

(3) A relatively long record of flow duration is 

available for Pheasant Branch Creek at the 

Lake Mendota Outlet, but these data were not 

collected contemporaneously with this study. 

Moreover, it is difficult to measure discharge 

accurately at this location because of lake 

backwater effects (D. Graczyk, USGS, oral 

commun., 1999). Therefore, these measurements were given less weight (coefficient of 

variation = 0.5).

(4) Head measurements in the TMR model domain 

were given a low weight due to the uncertainty 

regarding their representativeness for the conditions simulated during calibration. These 

head data are the sum of all measured water 

levels for the area, but were not collected contemporaneously. Rather, these data span 

40 years. In addition, the head data are also 

less precisely located, both horizontally and 

vertically. Therefore, the resulting unfiltered 

head target data set often had multiple head 

values for a single node. Moreover, the measured head values for a single node might differ 

by over 100 ft. Clearly a finite difference 

model (such as one head value calculated per 

node) cannot simulate these data. The measured head data were filtered for use in the 

TMR model so that nodes with multiple measured values were replaced with the average of 

all head data for the node. Given these uncertainties, layer 2 head targets had a standard 

deviation equal to 10 ft, which equate to a 

95-percent confidence interval of ± 20 ft 

around the measured value. Layer 3 head targets were weighted using a standard deviation 

equal to 25 ft, representing the additional 

uncertainty resulting from the wells being open 

to multiple aquifers and the potential unsteady 

nature of water levels in high capacity wells. 

The parameters initially chosen for optimization 

included the hydraulic conductivity (6 zones in layer 1, 

1 zone in layers 2 and 3), vertical leakage of nodes 

immediately beneath the Springs, and the conductance 

of Lake Mendota littoral and deep lake sediments. Initial runs on parameter sensitivity showed the model was 

insensitive to changes in layer 1 hydraulic conductivity 

and lake bed conductance (fig. 4); that is, the measured 

observations used in the optimized calibration did not 

contain enough information to constrain these parameters. As a result, all subsequent runs used fixed values 

for these parameters based on the Dane County 

Regional Model and optimized only the remaining sensitive parameters. 

Stochastic Monte Carlo Analyses

The effect of parameter uncertainty can be more 

formally addressed using stochastic approaches. While 

detailed discussion of stochastic techniques is beyond 

the scope of this work, a brief discussion follows. A 

Monte Carlo approach was used as a means to obtain the 

probability distribution of the capture zone of the 

Springs. This approach allows calculation of the probability of a certain occurrence, in this case the probability 

that the Springs will capture water from different parts 

of the model domain, given the uncertainty that exists 
10 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

for a discrete set of parameters. In this approach, a large 

number of MODFLOW model runs using the TMR 

model were performed using Stochastic MODFLOW 

(Ruskauff and others, 1998) while randomly varying 

the horizontal hydraulic conductivity of the upper bedrock aquifer (layer 2) and lower bedrock aquifer 

(layer 3) over reasonable ranges of values based on the 

results of UCODE optimization. In the case of layer 2, 

this range was uniformly distributed (all values are 

equally likely) between 1 and 15 ft/d; layer 3 was uniformly varied between 0.7 and 10 ft/d. Each run is 

called a “realization,” and reflects one possible set of 

parameters for the model.

Because combinations of reasonable parameter 

values may yield unrealistic results, the head calibration 

statistics of the realizations were evaluated (or “conditioned”) and unreasonable realizations were removed. 

Particle tracking was performed using Stochastic 

MODPATH (Ruskauff and others, 1998) in the Pheasant 

Branch watershed. This code utilizes output from the 

Stochastic MODFLOW realizations to delineate the 

recharge area of the Springs using a probability distribution where 1 is contributing in 100 percent of the realFigure 4. Plot of parameter sensitivity from the UCODE parameter estimation. (Horizontal hydraulic conductivity zones 

representing layer 2, and to a lesser extent layer 3, were the most sensitive parameters. Additional work with optimization

and stochastic runs focused on modifications to these two parameters.)

2E-03 3E-03 6E-02 6E-02 0.1 0.2

6.8

0.6

6E-04 4E-03 2E-02 1E-02

0

1

2

3

4

5

6

7

8

COMPOSITE SCALED SENSITIVITIES (DIMENSIONLESS)

layer 1 Kh 1layer 1 Kh 2layer 1 Kh 3layer 1 Kh 4layer 1 Kh 5layer 1 Kh 6

layer 2 Kh

layer 3 Kh 

Mendota littoral conductance

Mendota deep lake conductance

Spring layer 1-2 leakageSpring layer 2-3 leakage
RESULTS OF MODEL SIMULATIONS AND GEOCHEMICAL INVESTIGATION 11

izations, and 0 is contributing in 0 percent of the 

realizations.

Geochemical Investigation

Geochemical investigation focused on Springs and 

a nearby municipal well. The Springs area was divided 

up into 8 areas in the main spring complex, one additional spring located 1,300 ft west of the main spring 

complex, and one ephemeral stream. Samples were collected for analysis of both ion and isotope chemical constituents. Major ions and nutrients were measured 

periodically from the Springs during March 1998 

through April 1999. The municipal well was sampled 

once in August 1998. Springs and stream were sampled 

using a peristaltic pump; the municipal well was sampled from the pump well head. Unfiltered samples were 

used for field measurements of conductivity, temperature, dissolved oxygen and pH, and lab measurements 

of alkalinity. Filtered samples (0.45 μm cellulose nitrate 

filter) were collected for determination of major ions, 

nitrate + nitrite (NO3 + NO2), ammonia (NH4), total N 

and total dissolved inorganic phosphorus (DIP) analyses. 

Because water and strontium isotope chemistry are 

not widely used in hydrological investigations, a short 

description, taken from Hunt and others (1998), is given 

here. Water isotopes (oxygen and deuterium) are ideal 

conservative tracers of water sources because they are 

part of the water molecule itself. Stable isotopes of 

water are conservative in aquifers at low temperature, 

but fractionate on the surface at less than 100 percent 

humidity (Gat, 1970). Because the vapor pressure of 

H2

16O is greater than H218O, the residual liquid is characterized by a higher H2

18O content after evaporation. 

Hydrogen and deuterium also fractionate, but to a 

greater extent due to larger percent mass difference. 

Thus, characteristic 18O/16O and 2H/1H ratios can fingerprint water sources. Strontium-87 is produced from 

radioactive decay of rubidium-87. Assuming a given 

initial 87Sr/86Sr, minerals that have high Rb/Sr concentration ratios will attain higher 87Sr/86Sr than minerals 

that have low Rb/Sr concentration ratios. Biological or 

low-temperature abiotic processes do not significantly 

fractionate Sr isotopes; the isotopic composition of Sr 

(unlike that of the lighter elements) is entirely controlled by the mixing of Sr from geologic regimes having different isotopic composition (Graustein, 1989). Sr 

isotopes, when considered together with ion chemistry, 

can thus distinguish weathering reactions from cation 

exchange processes. This ability makes them suitable 

for water-rock interaction investigations (such as Bailey 

and others, 1996; Bullen and others, 1996) and discerning sources of water from isotopically distinct source 

areas (Eastin and Faure, 1970; Fisher and Steuber, 

1976). 

The major ions (Ca, Mg, Na, K, SO4), total P, 

nitrate+nitrite, NH4, and total N were analyzed by the 

University of Wisconsin Soil and Plant Analysis Laboratory and State Laboratory of Hygiene during the periods March 1998–October 1998 and November 1998–

April 1999, respectively. Alkalinity was analyzed by 

gran titration at the University of Wisconsin Water 

Chemistry Department. Chloride analyses were performed using liquid chromatography at the USGS, Middleton, Wisconsin. Analyses of water and strontium 

isotopes were performed at the U.S. Geological Survey 

National Research Program Laboratory in Menlo Park, 

California. Oxygen-18 values were measured using 

CO2-H2O equilibration (Epstein and Mayeda, 1953) on 

a Finnigan-Mat 251 mass spectrometer. Oxygen-18 values are reported in standard delta notation relative to 

Vienna Standard Mean Ocean Water (VSMOW). Analytic error (2σ) is estimated at ± 0.1 per mil. Strontium 

isotope analyses were performed using the methodology described by Hunt and others (1998).

RESULTS OF MODEL SIMULATIONS AND 

GEOCHEMICAL INVESTIGATION

Telescopic Mesh Refinement

TMR modifications to the DCRM improved head 

and flux calibration for the Pheasant Branch watershed. 

The inclusion and more accurate depiction of nearby 

surface-water features (such as Brewery Creek, Dorn 

Creek, headwaters of Black Earth Creek) improved the 

DCRM head calibration in the area encompassed by the 

TMR model (table 2). The TMR area of the DCRM 

domain was slightly better calibrated than the overall 

DCRM model (root mean squared difference = 36.1 ft 

in TMR area, 37.4 ft for the overall domain—Krohelski 

and others, 2000). More notably, the recharge rate 

derived by the surface-water modeling significantly 

improved the flux calibration while only minimally 

degrading the head calibration (first TMR column versus last TMR column in table 2). Given the uncertainty 

associated with the accuracy of some of the head targets 

(such as significantly different observed values reported 

for the same model node) and the significant improve-
1

2SimulationoftheRechargeAreaforFrederickSprings,DaneCounty,Wisconsin

Table 2. Results of the Pheasant Branch UCODE runs using telescopic mesh refinement (TMR) model and Dane County Regional Mode[NA, not applicable; residual = measured - simulated; VCONT, MODFLOW VCONT (McDonald and Harbaugh); L2/L3 Kh, layer 2, layer 3 horizontal hydraulic conduPackage (Prudic, 1989); PRMS, Precipitation Runoff Modeling System (Leavesley and others, 1983)]

Sensitivity Runs Optim

Modified DCRM TMR model DCRM R and K TMR model new base TMR model head targets only 

Weighted sum of squares NA 5,276 5,659 5,864 Unweighted head residuals

Residual mean -1.1 4.6 1.0 5.1 Root mean squared difference 36.1 33.7 34.7 33.2 Mean absolute error 26.5 23.8 25.2 23.3 

Unweighted flux residuals

All flux targets—mean absolute error (ft3/s) NA 1.3 0.3 0.7 Flux residual @ Highway 12 (ft3/s) 0.6 1.1 .3 .5 Flux residual @ Springs (ft3/s) 1.0 .2 -.01 .4 

Optimization parameters

Layer 2 Kh (ft/d) 5.0 5.0 5.0 15.0 Layer 3 Kh (ft/d) 10.0 10.0 10.0 .7 

Model design Orig DCRM VCONT

Orig DCRM streamsDCRM Recharge arrayDCRM L2/L3 KhSprings added to DCRM

Spring fracture zoneTMR modified STR pkgDCRM Recharge arrayDCRM L2/L3 Kh

Spring fracture zoneTMR modified STR pkgPRMS Recharge arrayDCRM L2/L3 Kh

Spring fracture zoneTMR modified STR pkPRMS Recharge arrayUCODE L2/L3 Kh
RESULTS OF MODEL SIMULATIONS AND GEOCHEMICAL INVESTIGATION 13

ment in simulated flux, the slight degradation in head 

calibration was considered to not invalidate the 

recharge array derived from the surface-water modeling. 

The basin recharge rate is considered relatively 

well known (based on the surface-water modeling 

work), thus the recharge rate was fixed and horizontal 

hydraulic conductivity (Kh) became the driver for head 

calibration. The DCRM used a sedimentologically 

based method for assigning conductivities (Swanson, 

1996; Krohelski and others, 2000) that resulted in a 

large number of horizontal hydraulic conductivity (Kh) 

values to represent the conductivity field of layer 1. A 

large number of values is not suitable for parameter estimation because each value requires three model runs to 

perform sensitivity analysis or optimization, which in 

turn results in unacceptably long run times and timeintensive data handling. As a means to assess the effects 

of changes to hydraulic conductivity, the Kh values 

were grouped into eight Kh zones (table 1). Differences 

between the original DCRM Kh field and the Kh field 

created from the zones were negligible (that is, calibration statistics were identical to two significant figures). 

A sensitivity analysis was run on the eight Kh zones and 

demonstrated that only the zone representing the Kh of 

the upper and lower bedrock (layers 2 and 3) significantly affected the model results (fig. 4). This is not 

surprising in that much of layer 1 is dry in the Pheasant 

Branch basin, and all head targets are located in layers 2 

and 3. 

Parameter Estimation

The smaller number of Kh zones allowed investigation of the TMR model using the parameter estimation 

program UCODE. During the optimization process a 

difficulty was encountered. The automatic routines of 

UCODE that perturb parameters and run MODFLOW 

in batch mode have an unfortunate artifact—the perturbation required by UCODE can cause the stable base 

MODFLOW model to become unstable and not converge. This in turn will result in a failed UCODE run. If 

the modeler chooses to continue on without convergence (an option in Groundwater Vistas) the UCODE 

routine will successfully complete. However, the 

UCODE optimization evaluation may include a nonconverged model run that ended with an iteration that 

yields an unrealistic model (such as poor mass balance). 

Because the model is re-run and model output files 

overwritten, it can be difficult to discern if non-convergence occurred. In the TMR model runs, non-convergence was caused by location and high conductance of 

a subset of fixed-stage stream nodes in the Stream Routing Package (STR). High values of stream conductance 

effectively transmit the effects of the headwater stream 

nodes to the ground water (and vice versa), which can 

cause oscillations in the ground-water solution. This 

can cause segments of headwater streams to dry and wet 

sporadically. The model solution at the end of the maximum number of iterations often had inaccurate simulated stream flows that were used by UCODE to 

formulate the new (poorly based) parameter values. 

Removing a small set of headwater STR nodes overcame the instability at higher Kh values. It should be 

noted that the problems would likely be worse if constant flux or general head boundary conditions are used 

instead of constant head nodes at the model perimeter 

because constant head boundary conditions fix heads at 

the perimeter, that, in turn, tend to dampen the oscillation within the model domain.

When only head targets were used to obtain the 

optimal calibration, the head residual was minimized 

(fig. 5) and resulting Kh values were three times higher 

than the DCRM for layer 2, and over ten times lower for 

layer 3 (table 2). As might be expected when flux targets are not considered, the measured Q50 fluxes are 

poorly simulated (fig. 5). When heads and flux targets 

are considered the optimization routine obtains a minimum in the residual sum of square error (table 2, fig. 5) 

but increases for Kh of layer 2 to nearly twice the 

DCRM value, and decreases Kh of layer 3 to roughly 

one-fifth the DCRM value (table 2). This yields the 

“optimal” solution for these targets and weights, but 

should be considered with the following caveats.

Unweighted calibration results might indicate 

biases in the head data used for calibration. For example, the mean error values (table 2) indicate simulated 

head values are consistently too low (the residual sign is 

positive) in all TMR model runs. While the result could 

be interpreted as the need for increased recharge and/or 

lower Kh, it could also be a result of biases in the measured head data and inconsistencies in how well they 

represent the regional water levels. The latter might be 

expected because the TMR model domain contains a 

large extent of the end moraine and Driftless areas. 

These areas are expected to have more perched ground 

water (Krohelski and others, 2000), and are likely 

responsible for the large range of water levels reported 

for a single model node. In addition, perched water level 

errors are always biased towards higher heads (that is, 
14 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

Figure 5b. Plot of residuals from all head calibration targets and the flux residual at Pheasant Branch 

Creek at Highway 12.

Figure 5a. Plot of sum of squares residual for UCODE optimization. (Note that although there are 

many more head targets than flux targets, the weight assigned the flux target at Highway 12 results

in an overall worse fit when only head targets are considered.)

5000

5100

5200

5300

5400

5500

5600

5700

5800

5900

6000

SUM OF SQUARES

Dane County

Regional Model

Optimized to 

heads and 

fluxes

Optimized to 

heads only 

Dane County 

Regional Model and 

new recharge array

Head Absolute Residual Mean, in feet x 5

Highway 12 residual, in cubic feet per second x 100

0

50

100

150

200

250

Optimized to 

heads and 

fluxes

Optimized to 

heads only 

HEAD AND FLUX RESIDUAL

Dane County

Regional Model

Dane County 

Regional Model and 

new recharge array
RESULTS OF MODEL SIMULATIONS AND GEOCHEMICAL INVESTIGATION 15

perched water levels will always be above the regional 

water table). While not within the scope of this work, 

critically evaluating and filtering the head calibration 

targets is a subject of ongoing work. The parameter estimation routines use weighted statistics based on a sum 

of squares residual calculation that is more resistant to 

the biases described above. 

There is other evidence that suggests that the Kh of 

layer 2 and 3 is similar to that estimated by UCODE in 

the model area. Layer 2 bedrock includes a carbonate 

aquifer that outcrops throughout the area and is characterized by extensive fracturing and dissolution; conductivities for Cambrian sandstones elsewhere in Wisconsin have also been on the order of that reported by the 

UCODE optimization (Young, 1992). However, it is 

also possible that the optimization might have been 

overly controlled by the high weight given to the Highway 12 gaging station data. For example, layer 3 has 

only five head targets, and they are dominated by water 

levels from high capacity wells. The heads reported are 

uncertain (thus given lower weight in the optimization) 

due to problems with representativeness of water levels 

obtained from pumping wells (Driscoll, 1986). These 

problems are a result of confounding factors such as 

well-loss effects, and the fact that the wells are often 

open to multiple aquifers (Krohelski and others, 2000). 

In this case, the optimization of layer 3 is being affected 

more by how changes in Kh affect the distribution of 

flow between locations measured by the flux targets and 

deep underflow to Lake Mendota than the head calibration in layer 3. 

Flux simulation was also refined in the TMR modeling. Preliminary investigation was performed using 

the DCRM and included both changes to the recharge 

array and the degree of vertical leakage between 

layers 2 and 3 near the Springs. The use of the DCRM 

was advantageous because the Springs are a regional 

hydrologic feature, thus are suitable for the larger grid 

spacing of the DCRM. Any changes that were required 

could be first tested for efficacy and reasonableness 

without the concern of boundary violations and added 

work of translating all changes to the TMR model. 

Based on this work, modifications to the recharge array 

derived from the surface-water model were deemed reasonable and simulated flows at the Highway 12 gaging 

station were much improved. Spring discharge was 

improved by the new recharge array although it was 

noted that Spring flows could not be simulated using 

only layer 2 capture. It was concluded that maintaining 

a tight confining unit between the lower bedrock aquifer 

and upper bedrock aquifer in the area of the Springs 

could not yield sufficient water to simulate measured 

spring discharge given the present understanding of the 

recharge distribution. An increase to the degree of vertical leakage between the deep sandstone aquifer 

(layer 3) and the upper Paleozoic bedrock (layer 2) was 

investigated using the DCRM; this work indicated that 

a breach in the Eau Claire confining unit could explain 

high Spring discharge. Modifications to the degree of 

leakage between layers 2 and 3 were further investigated using the TMR model.

Spring discharge was improved by the new 

recharge array and modifications to the location of the 

Eau Claire Shale confining unit. In the TMR model, 

modifications to layer 2–3 vertical hydraulic conductivity (Kv) were restricted to nodes near the Spring 

(table 1, fig. 6). In the larger springs area, 16 TMR 

nodes were assigned a vertical leakage equal to a value 

typical for the Lake Mendota area where the Eau Claire 

shale is absent (table 1, fig. 6). In addition, one TMR 

node immediately beneath the Spring was assigned a 

higher Kv (such as might be expected from a prominent 

vertically orientated fracture zone); the vertical leakage 

value assigned to this node (VCONT = 0.14, table 1 and 

fig. 6) roughly corresponds to assigning the Kv of 

layers 2 and 3 a value five times that used in the DCRM. 

While the entire area modified is restricted to the size of 

one DCRM node (1,312 ft by 1,312 ft), the exact location where the Eau Claire confining unit is absent is 

poorly understood in this area. It is conceivable that this 

area extends from the Springs to the area near Lake 

Mendota where the Eau Claire was missing in the 

DCRM (fig. 6). This is a potential topic for investigation, and may be addressed in future work in the Pheasant Branch.

Overall simulation of fluxes in the TMR model was 

significantly improved from the unmodified DRCM 

(table 2). Flux at the Highway 12 gaging station 

improved from being 60 percent under-simulated to 

within measurement error (11 percent under-simulated) 

in the optimized case. Further improvement in simulated fluxes by UCODE optimization was limited by the 

parameters chosen for optimization. That is, recharge 

was not varied in the UCODE runs, and optimization 

proceeded by changing only hydraulic conductivities of 

layers 2 and 3. This approach is not as powerful as 

changing recharge rates, and global optimization is 

obtained by balancing the improved head calibration 

against degrading the flux calibration at Highway 12. 

Moreover, the simulated flux in the Springs was more 
16 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

Ephemeral stream

node directly below

spring, VCONT = 0.14 1/day 

larger Springs area,

VCONT = 4.46e-3 1/day 

Surface water divide

Middleton Well #4 

Head calibration point

EXPLANATION

Area where Eau Claire shale is absent

12

L a

k e

M e n d o t a

Figure 6. Area where the connection between layers 2 and 3 was modified in the telescopic refinement model is shown 

above. (The “larger Springs area” represents one DCRM node with vertical leakage equal to areas where the Eau Claire 

confining unit is absent (for example, Lake Mendota). The “node directly below Springs” corresponds to one TMR node 

that represents a vertical fracture conduit enhanced by dissolution. Larger VCONT values correspond to larger leakage.)
RESULTS OF MODEL SIMULATIONS AND GEOCHEMICAL INVESTIGATION 17

closely simulated in the runs that used the DCRM 

higher Kh of layer 3, but simulated flow at Highway 12 

was worse. This case was not considered the optimal 

model because (1) the flow duration record at the Highway 12 gaging station (24 years) was of better quality 

than the Spring discharge record (2 months), and (2) the 

surface-water modeling work and resulting recharge 

array focused on the areas upstream of Highway 12; the 

Springs capture zone encompasses a significant area 

outside the area studied in the surface-water modeling 

work (see below). These areas outside of the Pheasant 

Branch surface-water basin were, by default, assigned 

the lower flux rates of the original DCRM recharge 

array. 

Given the model results presented above, it is 

apparent that there is not a single set of parameters that 

can best simulate all facets of the Pheasant Branch system. The regional DCRM had an acceptable calibration 

to head and flux over the entire county. Calibration of 

the TMR model to heads alone indicates that layer 2 

should have higher conductivity and layer 3 should 

have lower conductivity for the Pheasant Branch area 

than specified in the DCRM. Calibration to heads and 

fluxes suggests that the DCRM values should be altered 

less than suggested by a calibration to heads alone, but 

the values are still different than the DCRM. None of 

the values can be easily dismissed as unreasonable, and 

different parameter sets might be explained by differences in scale, preferential flowpaths/spatial heterogeneity of hydraulic conductivity, and uncertainty in the 

measured data (such as wide variety in measured head 

for a single node, short flow duration records from some 

flux targets). While such areas will be topics of future 

work, stochastic methods can be used with the present 

understanding to address the effects of this parameter 

uncertainty on predictions. 

Stochastic Monte Carlo Analyses

Stochastic methods are ideal for incorporating the 

uncertainty in Kh for layers 2 and 3 because they are 

designed around the concept of probability. Similar to 

flipping a coin a large number of times to obtain probability of heads versus tails, a large number of model 

simulations can be run varying Kh to assess the probability of different head and flow configurations. By 

analyzing the set of model results, an assessment of 

the results and the uncertainty in predictions can be 

obtained. In this application the Kh of layer 2 and 

layer 3 were allowed to vary within the range specified 

by the UCODE optimization—between 5 and 15 ft/d for 

layer 2 and between 0.7 and 10 ft/d for layer 3. Because 

there are no data to suggest which of the values are more 

likely, the simulations were set up such that all values 

within the range were equally likely (that is, a uniform 

distribution). There are no set rules for the number of 

runs needed to adequately characterize the problem, 

however an analysis of changes to the mean absolute 

error (MAE) in head showed that the problem had stabilized at 200 realizations (fig. 7a). A check of the 

model using 300 realizations and a different random 

number seed also confirmed that 200 realizations were 

sufficient. Examination of the variability in absolute 

mean head difference for 200 realizations shows the 

presence of outlier combinations that clearly do not 

yield a calibrated model (fig. 7b). These outliers point 

out the need for filtering the 200 runs such that only reasonable models are left for predictions.

This filtering (or “conditioning” in Monte Carlo 

parlance) was performed using the head MAE. The 

UCODE results in table 2 are considered reasonable 

model results—of those results the highest acceptable 

MAE was 25.2 ft. Therefore, this was the cutoff value 

for conditioning; realizations with an MAE less than or 

equal to 25.2 were passed and included in the predictive 

runs. Of the original 200 realizations, 136 realizations 

met this criterion. The automated routine in Groundwater Vistas can only condition based on head results, 

therefore the reasonableness of these 136 realizations 

with respect to flux calibration was not assessed. However, this was not considered problematic because none 

of the simulated fluxes using the surface water-derived 

recharge array in table 2 are grossly out of calibration, 

and the effects on simulated fluxes are expected to be 

less sensitive than if recharge rates were included in the 

stochastic runs.

The stochastic runs were summarized by performing statistics and particle tracking on the conditioned 

runs. Mean heads are the average heads in the 136 realizations; the standard deviation addresses the variability 

of simulated heads. These results show a reasonable distribution of head, with the largest variation in simulated 

heads away from surface-water features or near pumping wells (figs. 8a and 8b, respectively). The conditioned heads and cell-by-cell flows can also be used to 

perform particle tracking with a stochastic version of 

the USGS program MODPATH. In this application the 

program computes the probability that specified particles are captured by a boundary condition of interest. In 

the stochastic code used here, no MODPATH variables 
18 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

Figure 7b. Scatter plot of mean absolute error (heads) for each of the 200 realizations of Stochastic MODFLOW. 

(Horizontal hydraulic conductivity (Kh) of layers 2 and 3 varied between 5 and 15 ft/d and 0.7 and 10 ft/d respectively. 

The plot illustrates the presence of outlier combinations of Kh that result in an uncalibrated model.)

Figure 7a. Cumulative mean absolute error (heads) for 200 realizations of Stochastic MODFLOW. (Horizontal 

hydraulic conductivity of layers 2 and 3 varied between 5 and 15 ft/d and 0.7 and 10 ft/d, respectively.) MEAN ABSOLUTE ERROR, IN FEET

23.8

24.0

24.2

24.4

24.6

24.8

25.0

25.2

25.4

25.6

0 50 100 150 200

REALIZATION NUMBER

REALIZATION NUMBER

20

25

30

35

40

45

0 20 40 60 80 100 120 140 160 180 200

MEAN ABSOLUTE ERROR, IN FEET
RESULTS OF MODEL SIMULATIONS AND GEOCHEMICAL INVESTIGATION 19

are considered uncertain and all uncertainty is assumed 

to occur in the MODFLOW results. In this application, 

capture probability was calculated using particles 

placed near the top of layers 2 and 3 in the northern half 

of the model domain; the stream node representing the 

Springs was used to tabulate capture.

Based on the model results, the Springs’ recharge 

area includes areas outside its immediate surface watershed (figs. 9a and 9b) demonstrating that the ground 

watershed and surface watershed do not coincide. The 

majority of the Springs’ recharge area is within the 

larger Pheasant Branch Creek watershed illustrating the 

strong connection between the North Fork of the Pheasant Branch and the Springs. The fact that the DCRM 

capture zone for layer 2 (fig. 1) approximates the TMR 

capture zone for layer 2 (fig. 9a) underscores the 

regional nature of problem. The longest travel times 

(from the farthest reaches of recharge area) are on the 

order of thousands of years. The stochastic results show 

that the uncertainties in Kh of layers 2 and 3 have relatively isolated effects on the north and west edges of the 

Springs’ recharge area. Moreover, the layer 3 capture 

zone is larger than and encompasses all of the layer 2 

capture zone demonstrating the reduced effects of competing water sinks (such as nearby streams) in the 

deeper layer.

Although the global hydrologic leakage between 

layer 2 and 3 was not explicitly included in sensitivity/optimization/stochastic runs, the larger layer 3 capture zone can be thought of as a conservative estimate 

of the Springs’ capture. Changes to the vertical leakage 

between layers 2 and 3 were not formally considered 

due to a lack of data on the Eau Claire confining unit 

properties/occurrence and sparse head data in the deep 

sandstone aquifer (layer 3). As a result the TMR model 

used unmodified calibrated values from the DCRM. 

However, it is conceivable that there is a higher leakage 

between the layers. If this were the case, water could 

flow to the deeper layer instead of laterally flowing to 

competing water sinks (such as Dorn and Pheasant 

Branch Creek) in areas where there is capture in layer 3 

but not layer 2; the layer 2 capture zone shown in 

figure 9a would be more applicable if the leakage 

between the layers 2 and 3 were smaller (tighter). 

Because the layer 3 capture zone encompasses the 

whole of the layer 2 capture zone, it is more likely to 

encompass the range of leakage between layers 2 and 3.

Geochemical Investigation 

The water in the Pheasant Branch Springs area is 

calcium-magnesium-bicarbonate type (table 3), the 

most common type in natural waters (Stumm and Morgan, 1981). Locally, other ions can be important (such 

as nitrate in Frederick Springs 1 and 2, table 3). As demonstrated by the small standard deviations reported in 

table 3, the variation in water chemistry over time at a 

given sampling location in the Pheasant Branch Springs 

area is relatively small. There are, however, notable differences in chemistry sampled at the Frederick Spring 

complex, the Gate Spring, and the Culvert Crossing surface-water runoff location (fig. 10) that are described 

below.

One of the most notable findings of the geochemical sampling was the large spatial variability of the Frederick Springs throughout the sampling period; there 

were large differences observed between the west (Frederick Springs 1 and 2) and east (Frederick Springs 5 

and 6) portions of the Springs, even though the spring 

vents 1 and 6 are located within 50 ft of each other 

(table 3, figs. 11 through 14). The difference in spring 

chemistry is gradational such that Frederick Springs 3 

and 4 are intermediate between the west and east spring 

water type. Moreover, this pattern holds for total dissolved solids as measured by specific conductance 

(fig. 11), major ions (such as fig. 12), as well as strontium (fig. 13) and oxygen (fig. 14) isotopes. In most 

cases the variability within the Frederick Springs is 

larger than that between Frederick Springs and a spring 

that was sampled 1,300 ft away (Gate Spring, table 3). 

The water chemistry differences observed in the 

springs are attributed to changes in location in the 

recharge area. Ground-water flow occurs in “flow 

tubes” that extend from the recharge areas to the discharge point. These tubes do not cross, and can be delineated using flow net analysis (Freeze and Cherry, 1979) 

or particle tracking in numerical modeling. Once a contributing area is defined, the relation between recharge 

area and discharge point chemistry (that is, spring vent 

chemistry) can be made.

In the case of the Frederick Spring system, there is 

evidence for different conditions in the recharge areas 

influencing the water chemistry observed at the spring 

vent. The calcium distribution (fig. 12) and strontium 

isotope (fig. 13) plots can be combined to show a twocomponent mixing diagram where the two end members are Ordovician carbonate dominated water in the 

western spring vents (Frederick Springs 1, 2a and 2b in 
20 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

Figure 8a. Results of conditioned stochastic MODFLOW runs for the upper bedrock aquifer (a1, model layer 2) and the lower 

bedrock aquifer (a2, model layer 3). (The figures show mean heads for the 136 realizations; the mean absolute error equals 

23.78 ft.)

1

240 Row 1 Column 160

Six Mile Creek

43°02'30"

Equipotential of lead 

contour—10-foot interval

High capacity well

EXPLANATION 012

012

MILES

KILOMETERS

89°30' 89°35'

47°07'30"

1

240 Row 1 Column 160

43°02'30"

89°30' 89°35'

47°07'30"

 860 

 880 

 880 

 900 

 900 

 900 

 920 

 920 

 920 

 920 

 940 

 940 

 940 

 880 

 880 

 900 

 900 

 920 

 920 

 920 

a1 a2
RESULTS OF MODEL SIMULATIONS AND GEOCHEMICAL INVESTIGATION 21

Figure 8b. Results of conditioned stochastic MODFLOW runs for the upper bedrock aquifer (b1, model layer 2) and the lower 

bedrock aquifer (b2, model layer 3). (The figures show standard deviation contours for the head distribution shown in figure 8a. 

The largest variations in head are near pumping wells and away from the surface water.)

Six Mile Creek

43°02'30"

012

012

MILES

KILOMETERS

89°30' 89°35'

47°07'30"

43°02'30"

89°30' 89°35'

47°07'30"

Standard deviation contours

for head distribution— 

2-foot interval 

High capacity well

EXPLANATION

 2 

 2 

 2 

 2 

 2 

 2 

 2 

 2 

 2 

 4 

 4 

 4 

 4 

 4 

 6 

 6 

 6 

 8 

 2 

 2 

 2 

 2 

 2 

 2 

 2 

 2 

 4 

 4 

 4 

 4 

 4 

 6 

 6 

b1 b2
22 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

Lake

Mendota

Probability

of capture 

1

240 Row

1 Column 160

43°02'30"

89°30' 89°35'

47°07'30"

0.00

1.00

18

151

12

14 EXPLANATION

Stream

Surface-water divide

012

012

MILES

KILOMETERS

Figure 9a. Results of conditioned stochastic MODPATH runs for the upper bedrock aquifer (layer 2). (The probabilistic capture 

zone for the Frederick Springs area is calculated by tracking particles placed near the top of layer 2 and tabulating if the particles are captured by the simulated spring. A value of 1 represents capture in 100 percent of runs. The capture zone is smaller 

than figure 9b because of stronger effects of competing hydrologic sinks (for example, other surface-water features) in layer 2

than in the deeper layer 3. In addition, the TMR model compares well to the capture zone delineated by the DCRM in figure 1).
RESULTS OF MODEL SIMULATIONS AND GEOCHEMICAL INVESTIGATION 23

Lake

Mendota

Probability

of capture 

1

240 Row

1 Column 160

43°02'30"

89°30' 89°35'

47°07'30"

Stream

Surface-water divide

18

151

12

14 EXPLANATION

012

012

MILES

KILOMETERS

0.00

1.00

Figure 9b. Results of conditioned stochastic MODPATH runs for the lower bedrock aquifer (layer 3). (The probabilistic capture 

zone for the Frederick Springs area is calculated by tracking particles placed near the top of layer 3 and tabulating if the pa rticles are captured by the simulated spring. A value of 1 represents capture in 100 percent of runs.)
2

4SimulationoftheRechargeAreaforFrederickSprings,DaneCounty,Wisconsin

Table 3. Results of geochemical sampling, 1998–99 (location of sampling sites are shown in figure 10)[n, number of samples; temp., temperature in degrees Celsius; D.O., dissolved oxygen; %, percent; mg/L, milligrams per liter; mmol/L, millimoles per liter; μmol, micromo

Sampling site n Temp. (C) D.O. (mg/L) Sat. (%) Conduc- tance (μS/cm) pH Ca mmol/L Mg mmol/L Na mmol/L K mmol/L Alkalinity mmol/L CO3 Cl mmol/L SO4 mmol/L NO3-

mmol

Gate Spring 1 average 8 9.6 4.9 43 714 7.1 1.98 1.62 0.45 0.04 3.14 0.87 0.20 0.44 

std dev 0.1 0.4 4 8 0.2 0.08 0.05 .03 .01 .07 .02 .01 .01 

Frederick Spring 1 average 8 10.4 7.7 69 787 7.0 2.19 1.80 .40 .04 3.16 .99 .29 .89 

std dev .8 .4 4 8 .2 .11 .08 .02 .01 .07 .02 .004 .09 

Frederick Spring 2A average 8 9.6 7.9 69 776 7.1 2.18 1.79 .32 .04 3.10 .83 .28 .94 

std dev .1 .2 1 8 .2 .10 .07 .02 .01 .07 .01 .002 .09 

Frederick Spring 2B average 8 9.6 7.8 68 779 7.0 2.18 1.80 .33 .04 3.13 .85 .28 .93 

std dev .1 .3 3 7 .2 .07 .07 .02 .01 .03 .02 .003 .08 

Frederick Spring 3 average 8 9.7 8.0 70 687 7.1 1.96 1.62 .32 .05 3.04 .64 .23 .57 

std dev .1 .4 4 34 .2 .08 .04 .01 .01 .04 .07 .01 .05 

Frederick Spring 3C average 6 10.1 8.2 72 636 7.1 1.80 1.50 .26 .06 3.00 .45 .20 .31 

std dev .7 .4 4 14 .3 .13 .08 .05 .01 .08 .06 .01 .06 

Frederick Spring 4 average 8 10.1 7.6 68 565 7.2 1.65 1.39 .20 .04 2.91 .23 .12 .17 

std dev .8 .3 3 8 .2 .05 .04 .01 .01 .07 .03 .01 .03 

Frederick Spring 5 average 8 10.1 6.8 60 512 7.2 1.51 1.30 .14 .02 2.86 .08 .06 .04 

std dev 1.1 .3 3 8 .2 .04 .05 .01 .01 .05 .01 .003 .02 

Frederick Spring 6 average 8 10.0 6.7 59 548 7.2 1.62 1.36 .20 .03 2.88 .22 .09 .15 

std dev .9 .1 2 13 .2 .12 .13 .05 .02 .07 .05 .01 .07 

Culvert Crossing average 2 14.1 4.8 44 475 7.5 .61 .61 .44 2.21 1.10 .80 .06 .01 

std dev 6.8 .9 NA 200 .04 .38 .40 .33 1.28 NA .58 .0005 NA 
RESULTS OF MODEL SIMULATIONS AND GEOCHEMICAL INVESTIGATION 25

fig. 15) and water that was influenced by a Cambrian 

sandstone component in the eastern vents (Frederick 

Springs 5 and 6 in fig. 15). The two components identified can be used to provide additional evidence for the 

capture zone calculated by the stochastic particle tracking described above. When the 50 percent probability 

capture zone is plotted on the bedrock map for Dane 

County (Massie-Ferch and others, in press), the western 

spring vents coincide with the western portion of the 

capture zone that is often capped by Ordovician carbonate (fig. 16). Carbonate rocks can significantly influence ground-water carbonate chemistry because of their 

relatively high solubility and rapid rates of dissolution. 

The northeastern portion of the capture zone feeds the 

east vents and is characterized by an absence of Ordovician carbonate, with the uppermost bedrock consisting 

of Cambrian sandstones (fig. 16). The eastern spring 

vents also have similar chemistry to Municipal Well #4 

(figs. 11, 12, and 14)—a well that is opened only to the 

Cambrian sandstone. The western area of Ordovician 

carbonate is larger than the northeastern area where it is 

missing, corresponding to the qualitative assessment 

that the flow from western spring vents was much larger 

than the flow from eastern spring vents.

The water isotope oxygen-18 samples from the 

spring vents also showed a difference between the western and eastern vents (fig. 14). The western vents were 

characterized by lighter 18O/16O ratios (had more negative values); eastern vents had slightly heavier ratios. 

There are a number of mechanisms that could explain 

this difference including factors such as different ages 

of water thus different climate and associated meteoric 

water line during recharge, amount of arable land in the 

recharge area (such as Darling and Bath, 1988), or differing ratio of snowmelt to rain precipitation in recharge 

water (Clark and Fritz, 1997). One possibility is that the 

elevation of the recharge area of the western vents is 

higher than that of the eastern vents. Others have noted 

Figure 10. Location of the Springs and geochemical sampling sites included in table 3 and figures 11–19. (The Frederick Spring 

complex consists of eight spring vents located in an 80 ft by 120 ft area. The spring vents are numbered clockwise starting from 

the western edge.)

Gate 

Springs

Culvert

Crossing 

Well #4

1

2a

2b

3

4

5

3c

6

flow to marsh

0 10 20 FEET

Frederick

Springs 

43°05'30"

43°07'

89°30'

0 0.5 1 MILE

0 0.5 1 KILOMETER

89°28'30"
26 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

0

100

200

300

400

500

600

700

800

900

Specific conductance

Specific conductance

(standard deviation)

SPECIFIC CONDUCTANCE,

IN MICROSIEMENS PER CENTIMETER

Gate 

Spring 1

Fred.

Spring 1 

Fred.

Spring 2a 

Fred.

Spring 2b

Fred.

Spring 3 

Fred.

Spring 3c

Fred.

Spring 4

Fred.

Spring 5

Fred.

Spring 6

Middleton

well #4

0.0

0.5

1.0

1.5

2.0

2.5

Calcium

Calcium

(standard deviation)

NA

Surface

 water

CALCIUM, IN MILLIMOLES PER LITER

Culvert

Crossing

Fred.

Spring 1 

Fred.

Spring 2a 

Fred.

Spring 2b

Fred.

Spring 3 

Fred.

Spring 3c

Fred.

Spring 4

Fred.

Spring 5

Fred.

Spring 6

Middleton

well #4

Figure 11. Average specific conductance measurements and standard deviations for 1998–99. (The largest standard 

deviation is observed in Frederick Spring 3, an intermediate spring between the two dominant water types; this observation 

suggests that the source of water to the spring may vary as transient stresses affect the system. It is notable that this is 

occurring in a regional discharge area.)

Figure 12. Average calcium concentrations and deviations for 1998–99. (The distribution of high calcium agrees well with 

the distribution of high specific conductance shown in figure 11.)
RESULTS OF MODEL SIMULATIONS AND GEOCHEMICAL INVESTIGATION 27

0.7080

0.7085

0.7090

0.7095

0.7100

0.7105

0.7110

0.7115

*from Hunt et al. 1998

Monroe County

STRONTIUM-87/STRONTIUM-86

Fred.

Spring 1 

Fred.

Spring 2a 

Fred.

Spring 2b

Fred.

Spring 3 

Fred.

Spring 3c

Fred.

Spring 4

Fred.

Spring 5

Fred.

Spring 6 Ordovician 

carbonate

(rock)*

Cambrian SS

porewater*

-9.2

-9.1

-9.0

-8.9

-8.8

-8.7

-8.6

-8.5

-8.4

-8.3

-8.2

Average Delta OXYGEN-18

DELTA OXYGEN-18 PER MIL

Gate 

Spring 1

Fred.

Spring 1 

Fred.

Spring 2 

Fred.

Spring 2B

Fred.

Spring 3 

Fred.

Spring 3C

Fred.

Spring 4

Fred.

Spring 5

Fred.

Spring 6

Middleton

well #4

Figure 14. Average δ180 for 1998–99. (Error bars represent 1 standard deviation for all samples measured. The municipal 

well was sampled only once in August 1998 therefore no error bar can be calculated. Analysis precision is plus or minus 

0.1 per mil.)

Figure 13. Average 87Sr/86Sr ratios for the March 1998 and September 1998 sampling. (Error bars reflect reported 

laboratory precision. Ratios obtained from a whole rock digestion of the Ordovician carbonate and porewater sampled 

from the Cambrian sandstone in Monroe County, Wisconsin, are also shown. Samples from the Gate Spring and the 

municipal well were not submitted for strontium isotope analysis.)
28 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

0.7091

0.7092

0.7093

0.7094

0.7095

0.7096

0.7097

0.7098

0.7099

1.40 1.50 1.60 1.70 1.80 1.90 2.00 2.10 2.20 2.30

Fred. Spring 1,

2a and 2b 

Fred. Spring 

5 and 6 

Fred. Spring 4

Fred. Spring 3c

Fred. Spring 3

Ordovician carbonate 

dominated

Cambrian Sandstone influenced

STRONTIUM-87/STRONTIUM-86

CALCIUM, IN MILLIMOLES PER LITER

the relationship between increasing elevation and more 

depleted (or lighter, more negative) 18O/16O ratios; this 

has been observed even in areas of minor relief, with the 

amount of change ranging from -0.15 to -0.5 per mil for 

every 100 m (328 ft) of elevation increase (Clark and 

Fritz, 1997). A rough estimate of the average elevation 

of the recharge area was calculated using the 30 m 

(98 ft) digital elevation model (DEM) for the basin and 

the estimated capture zones for the western and eastern 

vents. This calculation showed that, on average, the 

western portion of the capture zone was about 40 ft 

higher than the eastern. While this difference is not 

enough to explain observed differences, the discrepancy 

may be due to errors in the 30 m DEM or inaccuracies 

in the exact location of the capture zone for the east and 

west spring vents. Alternatively, the difference in elevation may only be one factor, and other controls such as 

age of recharge water, land use, or slope aspect may 

play a role in the composition of water isotopes measured in the spring vents.

Nutrient chemistry among the spring samples also 

was spatially variable. Dissolved nitrate+nitrite in the 

western spring vents were significantly higher than 

eastern vents (fig. 17), likely reflecting the effects of 

different agricultural practices in the recharge area. Others (such as Gambrell and others, 1975) have noted that 

a comparison of dissolved nitrate to dissolved chloride 

can indicate if denitrification was occurring along the 

flowline from recharge area to discharge area; no evidence of denitrification was found during this study 

(fig. 18). Finally, as is true in most natural systems, only 

low concentrations of total dissolved phosphorus were 

observed in the ground water. Surface-water concentrations of total dissolved P measured in the ephemeral 

stream (Culvert Crossing in fig. 19), were much higher 

than that observed in ground water. However, there are 

not enough data to discern if the temporally short but 

high concentration surface-water input is a more important source of nutrients to the Pheasant Branch Marsh 

system than the longer duration, low concentration, 

ground-water inputs.

Figure 15. Average 87Sr/86Sr ratios versus average calcium concentrations for 1998–99. (The samples represent a gradation 

between a two-component mixing where the end members include high calcium water with a strontium isotope ratio near that 

measured in Ordovician carbonates and lower calcium water with a strontium isotope ratio that suggests a larger contribution 

from Cambrian sandstones.)
RESULTS OF MODEL SIMULATIONS AND GEOCHEMICAL INVESTIGATION 29

Draft Bedrock map: 1:100,000 Bedrock Geology of Dane County, Wisconsin 

by K.M. Massie-Ferch, R.M. Peters, and B.A. Brown, in review, Wisconsin

Geological and Natural History Survey

calculated capture zone 

for eastern springs 12

14

89°30' 89°35'

47°07'30"

Cambrian bedrock

Ordovician bedrock

Stream

Surface-water divide

High capacity well

12

14

EXPLANATION

012

012

MILES

KILOMETERS

Middleton Well #4 

Figure 16. Results of particle tracking over the bedrock map. (The capture zone for eastern springs using the optimized model 

is shown within the 50 percent probability capture zone using all conditional realizations. The model results are overlain on a

bedrock map of the Pheasant Branch vicinity. Areas where the Ordovician carbonate is absent (northeast areas of the capture 

zone) have Cambrian sandstone bedrock in the recharge areas; these areas feed eastern spring vents at Frederick Springs. 

The southwestern areas of the capture zones are characterized by Ordovician carbonates and are areas that feed the western 

spring vents.)
30 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

Figure 17. Average nitrate + nitrite concentrations and standard deviations for the Pheasant Branch 

sampling, 1998–99. (The 10 mg/L nitrate-N drinking water standard is shown for reference.)

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

10 milligrams per liter nitrate-N

Gate 

Spring 1

Fred.

Spring 1 

Fred.

Spring 2a 

Fred.

Spring 2b

Fred.

Spring 3 

Fred.

Spring 3c

Fred.

Spring 4

Fred.

Spring 5

Fred.

Spring 6

Culvert

Crossing

NITRATE PLUS NITRITE, IN MILLIMOLES PER LITER

Nitrate plus nitrite as nitrogen

Nitrate plus nitrite standard deviation

0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.90

1.00

0.00 0.20 0.40 0.60 0.80 1.00 1.20

CHLORIDE, IN MILLIMOLES PER LITER 

straight line indicates no denitrification 

NITRATE PLUS NITRITE, IN MILLIMOLES PER LITER

Figure 18. Average nitrate + nitrite concentrations versus average chloride concentrations for the Frederick 

Spring sampling, 1998–99. (The straight line indicates that denitrification (which removes nitrogen but not 

chloride) is not occurring.)
SUMMARY AND CONCLUSIONS 31

SUMMARY AND CONCLUSIONS

The Pheasant Branch watershed includes Frederick 

Springs, a large spring complex that flows into Pheasant 

Branch Marsh. As the City of Middleton and its surroundings continue to develop, the watershed is 

expected to undergo significant urbanization. The 

effects of urbanization may reduce ground-water 

recharge and adversely affect down-gradient hydrologic 

features such as Frederick Springs. The effect of stormwater runoff and wetland loss on reduced ground-water 

recharge is complex because the surface-water system 

is coupled to the underlying ground-water system. As a 

result, it is difficult to reliably predict the effects of 

urbanization on stream baseflow and spring flows 

a priori. Identifying the source waters for the spring, 

however, is the first step in ensuring its protection. 

Therefore, the objective of this work is to use groundwater flow modeling and geochemical investigation to 

identify the source of water that sustains the spring 

complex.

The notable findings of this work can be summarized as follows:

• The linking of the ground-water model to the 

surface-water model gave higher confidence in 

the results of both models than if either had 

been used independently. Values of recharge 

calculated with the surface-water model 

improved flux calibration in the ground-water 

model. By linking the two approaches, the 

entire water budget (precipitation, evapotranspiration, baseflow, stormflow and groundwater recharge) is encompassed; this helps 

ensure that reasonable values are used for all 

parameters.

• Parameter estimation sensitivity analyses on 

the ground-water flow model demonstrated 

Figure 19. Average total filterable phosphorus concentrations and standard deviations for the Pheasant Branch sampling, 

1998–99. (As is true in the most natural systems, phosphorus is more efficiently transported in the surface water (Culvert Crossing). The Culvert Crossing surface water sampling site was only wet during the spring runoff event thus it is unclear if it is more 

important than the ground-water contribution on an annual scale.)

0

5

10

15

20

25

30

35

40

Gate 

Spring 1

Fred.

Spring 1 

Fred.

Spring 2a 

Fred.

Spring 2b

Fred.

Spring 3 

Fred.

Spring 3c

Fred.

Spring 4

Fred.

Spring 5

Fred.

Spring 6

Culvert

Crossing

Total phosphorus

Total phosphorus

standard deviation

TOTAL FILTERABLE PHOSPHORUS, IN MICROMOLES PER LITER
32 Simulation of the Recharge Area for Frederick Springs, Dane County, Wisconsin

that the calibration targets used in this study 

only supported changes in the upper bedrock 

aquifer (layer 2) and the lower bedrock aquifer 

(layer 3). Other potential parameter changes 

did not have significant effects on the calibration.

• Parameter estimation optimization of the 

ground-water flow model suggested that the 

horizontal hydraulic conductivity of the upper 

bedrock layer (layer 2 in the model) might be 

higher in the Pheasant Branch area than the 

global value that represented a best fit for the 

Dane County Regional Model (DCRM). The 

parameter estimation routine also suggests that 

the lower bedrock aquifer (layer 3) may have 

lower horizontal hydraulic conductivity than 

the global value used in the DCRM. The range 

was not exceedingly large (5 to 15 ft/d for 

layer 2, 0.7 to 10 ft/d for layer 3); insight into 

which values best represent the bedrock in the 

Pheasant Branch area was not readily available, nor within the scope of the project.

• Stochastic Monte Carlo simulations using the 

ground-water flow model formally addressed 

the uncertainty in the horizontal hydraulic conductivity identified in the parameter estimation 

optimization. In the Pheasant Branch model the 

Monte Carlo simulations sampled a range of 

horizontal hydraulic conductivity values spanning from 5 to 15 ft/d in layer 2 and 0.7 to 

10 ft/d in layer 3. The ranges were specified 

using a uniform distribution; thus, all values 

between the endmembers were equally likely. 

This approach allowed calculation of a probability distribution of the capture zone for the 

Springs.

• The calculated capture zone for the Springs 

showed that they are receiving water that was 

recharged from areas inside and outside of its 

immediate surface watershed. The capture 

zone encompassed the North Fork of Pheasant 

Branch basin, areas downstream of Highway 

12 in the Pheasant Branch surface watershed, 

and an area outside of the Pheasant Branch 

watershed. This result underscored the need for 

linking the surface-water model of the North 

Fork basin to the ground-water model of the 

Springs, even though the surface-water systems are in different basins.

• Geochemical sampling of the Frederick Spring 

complex showed very large differences in 

chemistry between the spring vents that were 

located within 50 ft of each other. The differences were stable in time, were present in both 

ion and isotope analyses, and showed a distinct 

gradation from high nitrate, high calcium, 

Ordovician carbonate dominated water in western spring vents to low nitrate, lower calcium, 

Cambrian sandstone influenced water in eastern spring vents. The difference in chemistry 

was explained by different bedrock geology in 

the recharge area as demonstrated by overlaying the 50 percent probability capture zone 

over a bedrock geology map for the area. This 

result gives additional confidence to the capture zone calculated by the ground-water flow 

model.

REFERENCES

Bailey, S.W., Hornbeck, J.W., and Driscoll, C.T., 1996, Calcium inputs and transport in a base-poor forest ecosystem 

as interpreted by Sr isotopes: Water Resources Research, 

v. 32, no. 3, p. 707–719.

Bradbury, K.R., Swanson, S.K., Krohelski, J.T., and Fritz A.K., 

1999, Hydrogeology of Dane County, Wisconsin: Wisconsin Geological and Natural History Survey Open-File 

Report 1999–04, 66 p.

Bullen, T.D., Krabbenhoft, D.P., and Kendall, C., 1996, Kinetic 

and mineralogic controls on the evolution of groundwater 

chemistry and 87Sr/86Sr in a sandy silicate aquifer, northern Wisconsin: Geochimica Cosmochimica Acta, v. 60, 

no. 10, p. 1807–1821.

Clark, I.D., and Fritz, P., 1997, Environmental Isotopes in 

Hydrogeology: Lewis Publishers, Boca Raton, Fla., 

328 p.

Darling, W.G., and Bath, A.H., 1988, A stable isotope study of 

recharge process in the English Chalk: Journal of Hydrology, v. 101, p. 31–46.

D’Onofrio, A.J., 1973, Pheasant Branch Marsh environmental 

study—Phase 1: D’Onofrio, Kottke and Associates, Inc. 

Middleton, Wis., 25 p.

Driscoll, F.G., 1986, Groundwater and Wells, Second Edition: 

Johnson Division, St. Paul, Minn., 1,089 p.

Eastin, R., and Faure, G., 1970, Seasonal variation of the solute 

content and 87Sr/86Sr ratio of the Olentangy and Scioto 

Rivers at Columbus, Ohio: Ohio Journal of Science, v. 70, 

p. 170–179.

Epstein, S., and Mayeda, T., 1953, Variation of 18O content of 

water from natural sources: Geochimica Cosmochimica 

Acta, v. 4, p. 213–244.
REFERENCES 33

Fisher, R., and Stueber, A.M., 1976, Strontium isotopes in 

selected streams within the Susquehanna River basin: 

Water Resources Research, v. 12, p. 1061–1068.

Gambrell, R.P., Gilliam, J.W., and Weed, S.B., 1975, Denitrification in subsoils of the North Carolina Coastal Plain as 

affected by soil drainage: Journal of Environmental Quality, v. 4, no. 3, p. 311–316.

Gat, J.R., 1970, Environmental isotope balance of Lake Tiberias: In Isotope Hydrology. Internal Atomic Energy 

Agency, Vienna, p. 109–127.

Grant, R.S., and Goddard G.L., 1980, Channel erosion and sediment transport in Pheasant Branch Basin near Middleton, 

Wisconsin—A preliminary report: U.S. Geological Survey Open-File Report 80–161, 19 p.

Graustein, W.C., 1989, 87Sr/86Sr ratios measure the sources 

and flow of strontium in terrestrial ecosystems: In Stable 

Isotopes in Ecological Research, New York, SpringerVerlag, p. 491–512.

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

estimating parameters of a transient, three-dimensional, 

ground-water flow model using nonlinear regression: 

U.S. Geological Survey Open-File Report 91–484, 358 p.

____1998, Methods and guidelines for effective model calibration: U.S. Geological Survey Water-Resources Investigation Report 98–4005, 90 p.

Holmstrom, B.K., Olson, D.L., and Ellefson, B.R., 1999, Water 

Resources Data Wisconsin water year 1998: U.S. Geological Survey Water-Data Report WI–98–1, p. 227–236.

Hunt, R.J., Bullen, T.D., Krabbenhoft, D.P., and Kendall, C., 

1998, Using stable isotopes of water and strontium to 

investigate the hydrology of a natural and a constructed 

wetland: Ground Water, v. 36, no. 3, p. 434–443.

Krohelski, J.T., Bradbury, K.R., Hunt, R.J., and Swanson, S.K., 

2000, Numerical simulation of ground-water flow in 

Dane County, Wisconsin: Wisconsin Geological and Natural History Survey Bulletin 98, 44 p.

Krug, W.R., and Goddard, G.L., 1986, Effects of urbanization 

on streamflow, sediment loads, and channel morphology 

in Pheasant Branch Basin near Middleton, Wisconsin: 

U.S. Geological Survey Water-Resources Investigations 

Report 85–4068, 82 p.

Lathrop, R.C., and Johnson, C.D., 1979, Dane County Water 

Quality Plan: Appendix B, Water Quality Conditions, 

359 p.

Leavesley, G.H., Lichty, R.W., Troutman, B.M., and Saindon, 

L.G., 1983, Precipitation-runoff modeling system—

user’s manual: U.S. Geological Survey Water Resources 

Investigations Report 83–4238, 207 p.

Maher, L.J., 1999, The early history of the Pheasant Branch 

watershed: Appendix 3, unpublished report from the 

North Fork Pheasant Branch Watershed Committee, 

http://www.geology.wisc.edu/~maher/pheasant_branch.html.

Massie-Ferch, K.M., Peters, R.M., and Brown, B.A., in press, 

Bedrock Geology of Dane County, Wisconsin: Wisconsin 

Geological and Natural History Survey, scale 1:100,000.

McDonald, M.G., and Harbaugh, A.W., 1988, A Modular 

Three-Dimensional Finite-Difference Ground-Water 

Flow Model: U.S. Geological Survey Techniques of 

Water-Resources Investigations Report, bk. 6, chap. A1, 

576 p. 

North Fork Pheasant Branch Watershed Committee, 1999, 

unpublished “A Report From The North Fork Pheasant 

Branch Watershed Committee”,

http://www.geology.wisc.edu/~maher/pheasant_branch.html.

Poeter, E.P., and Hill, M.C., 1997, Inverse models—A necessary next step in ground-water modeling: Ground Water, 

v. 35, no. 2, p. 250–260.

Poeter, E.P., and Hill, M.C., 1998, Documentation of UCODE, 

a computer code for universal inverse modeling: U.S. 

Geological Survey Water-Resources Investigation Report 

98–4080, 116 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.

Rumbaugh, J.O., and Rumbaugh, D.B., 1996, Guide to using 

Groundwater Vistas: Advanced model design & analysis, 

Environmental Simulations, Inc., 180 p. 

Ruskauff, G.J., Rumbaugh, J.O., and Rumbaugh, D.B., 1998, 

Stochastic MODFLOW for Monte Carlo Simulation: 

Environmental Simulations Incorporated, Herndon, Va., 

58 p.

Selbig, W.R., 1996, Pheasant Branch Marsh—Interpretive 

study for Spring of 1996: Unpublished report for Institute 

of Environmental Studies 361, University of WisconsinMadison, 36 p.

Stumm, W., and Morgan, J.J., 1981, Aquatic Chemistry: John 

Wiley and Sons, New York, 780 p.

Swanson, S.K., 1996, A comparison of two methods for estimating groundwater recharge in Dane County, Wisconsin: 

M.S. Thesis, Department of Geology and Geophysics, 

University of Wisconsin-Madison, 123 p.

Young, H.L., 1992, Hydrogeology of the Cambrian-Ordovician Aquifer System in the Northern Midwest, United 

States: U.S. Geological Survey Professional Paper 

1405–B, 99 p.