In cooperation with the Menominee Indian Tribe of Wisconsin

Simulation of the Regional Groundwater-Flow System 

of the Menominee Indian Reservation, Wisconsin

Scientific Investigations Report 2014–5237

U.S. Department of the Interior

U.S. Geological Survey

89° 88°50' 88°40' 88°30'

45°10'

45°

44°50'

OCONTO COUNTY

SHAWANO COUNTY

LANGLADE COUNTY

Wolf

River

Legend Lake Round

Lake

Long

Lake

Sand

Lake

MENOMINEE

INDIAN RESERVATION

Keshena

Middle

Village

Neopit

Onekewat

Zoar

Red

Wing

West Branch

Wolf River

Tourtillotte

Creek

South Branch

Oconto River

Red Wing

Keshena

Neopit

Zoar

Neopit

Red Wing

Middle Village Middle Village

Onekewat Onekewat
Cover image. Map of the Menominee Indian Reservation with inset maps showing the probabilistic 

contributing areas to production wells.
Simulation of the Regional Groundwater-Flow 

System of the Menominee Indian Reservation, 

Wisconsin

By Paul F. Juckem and Charles P. Dunning

In cooperation with the Menominee Indian Tribe of Wisconsin

Scientific Investigations Report 2014–5237

U.S. Department of the Interior

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

SALLY JEWELL, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

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

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

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

For an overview of USGS information products, including maps, imagery, and publications, visit http://www.usgs.gov/

pubprod

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

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

Government.

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

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

Suggested citation:

Juckem, P.F., and Dunning, C.P., 2015, Simulation of the regional groundwater-flow system of the Menominee 

Indian Reservation, Wisconsin: U.S. Geological Survey Scientific Investigations Report 2014–5237, 40 p., 

http://dx.doi.org/10.3133/sir20145237.

ISSN 2328-0328 (online)
iii

Contents

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

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

Purpose and Scope ..............................................................................................................................2

Physical Setting.....................................................................................................................................2

Data Sources..................................................................................................................................................4

Modeling Methods.........................................................................................................................................4

Conceptual Model of the Groundwater-Flow System .............................................................................6

Aquifers and Confining Units ..............................................................................................................6

Sources and Sinks of Water ...............................................................................................................8

Hydrogeologic Boundaries .................................................................................................................8

Hydraulic Properties of the Groundwater-Flow System .........................................................................8

Hydraulic Conductivity.........................................................................................................................8

Recharge................................................................................................................................................9

Riverbed Resistance...........................................................................................................................10

Groundwater Withdrawals.........................................................................................................................10

Simulation of the Regional Groundwater-Flow System.........................................................................10

Model Construction and Assumptions............................................................................................12

Model Calibration................................................................................................................................15

Sensitivity and Uncertainty Analysis...............................................................................................21

Evaluation of Simulated Results of the Menominee Indian Reservation Model...............................21

Simulation of Areas Contributing Recharge to Existing Wells .............................................................23

Simulation of Flow of Wastewater From Infiltration Lagoons ..............................................................29

Assumptions and Limitations.....................................................................................................................30

Summary and Conclusions.........................................................................................................................32

References Cited..........................................................................................................................................34

Appendix 1.....................................................................................................................................................38

Figures

1. Map showing hydrologic features of the Menominee Indian 

Reservation, Wisconsin...............................................................................................................3

2. Map showing the distribution of glacial deposits in the Menominee 

Indian Reservation, Wisconsin,..................................................................................................5

3. Representations of the groundwater-flow system of the Menominee Indian 

Reservation, Wisconsin. .............................................................................................................7

4. Map showing hydrologic features and base elevation and 

hydraulic-conductivity zones simulated with analytic elements, Menominee 

Indian Reservation, Wisconsin.................................................................................................11

5. Graph showing simulated water levels in relation to target water levels, 

Menominee Indian Reservation, Wisconsin ..........................................................................18

6. Graph showing simulated base flow and gains in base flow in relation to target 

base flow and gains in base flow, Menominee Indian Reservation, Wisconsin .............19
iv

Figures (continued)

7. Map showing water-level- and base-flow-target residuals for the 

groundwater-flow model of the Menominee Indian Reservation, Wisconsin .................20

8. Map showing simulated groundwater-flow directions and water-table contours, 

Menominee Indian Reservation, Wisconsin ..........................................................................22

9. Map showing simulated water-table contours and the area contributing 

recharge to production wells for the village of Zoar, 

Menominee Indian Reservation, Wisconsin ..........................................................................24

10. Map showing simulated water-table contours and the area contributing 

recharge to production wells for the village of Neopit, 

Menominee Indian Reservation, Wisconsin ..........................................................................25

11. Map showing simulated water-table contours and the area contributing 

recharge to the production well for Middle Village, 

Menominee Indian Reservation, Wisconsin ..........................................................................26

12. Map showing simulated water-table contours and the area contributing 

recharge to production wells for the villages of Keshena and Onekewat, 

Menominee Indian Reservation, Wisconsin ..........................................................................27

13. Map showing simulated water-table contours and the area contributing 

recharge to production wells for the village of Redwing, 

Menominee Indian Reservation, Wisconsin ..........................................................................28

14. Map showing showing simulated probability of the horizontal plume extent 

from the Neopit wastewater-treatment lagoons under current (2013) conditions, 

Menominee Indian Reservation, Wisconsin ..........................................................................31

Tables

1. Pumping-test data and TGUESS calculations for production wells near Zoar, 

Menominee Indian Reservation, Wisconsin ............................................................................9

2. Groundwater withdrawal rates for production wells simulated in the groundwaterflow model of the Menominee Indian Reservation, Wisconsin..........................................10

3. Calibrated and specified parameter values and composite sensitivities for the 

groundwater-flow model of the Menominee Indian Reservation, Wisconsin .................13

4. Calibrated hydraulic-conductivity parameters and associated base elevations 

used to calculate effective hydraulic-conductivity estimates for the 

Menominee Indian Reservation, Wisconsin ..........................................................................13

5. Adjusted porosity estimates used to simulate time-dependent areas contributing 

recharge for pumped wells, Menominee Indian Reservation, Wisconsin........................14

6A. Calibration-target descriptions and results for groundwater and lake-level 

targets, and associated weights used for calibration with the parameter 

estimation program PEST..........................................................................................................16

6B. Calibration-target descriptions and results for stream base flow and 

base-flow gain targets, and associated weights used for calibration with the 

parameter estimation program PEST ......................................................................................16
v

Conversion Factors

Multiply By To obtain

Length

inch (in) 2.54 centimeter (cm)

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Area

acre 0.4047 hectare (ha)

square mile (mi2

) 2.590 square kilometer (km2) 

Flow rate

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

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

cubic foot per day (ft3

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

gallon per day (gal/d) 3.785 liter per day (L/d)

gallon per year (gal/yr) 3.785 liter per year (L/yr)

million gallons per year (Mgal/yr) 3,785 cubic meter per year (m3

/yr)

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

Hydraulic conductivity*

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

Transmissivity**

foot squared per day (ft2

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

Horizontal flux

foot squared per day (ft2

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

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

°C=(°F-32)/1.8

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

(NAVD 88).

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

and the North American Datum of 1927.

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

*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, foot per day (ft/d), is used for convenience.

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

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

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

Acknowledgments

Local Project Coordinators

Gary Schuettpelz, Former Environmental Services Director, Menominee Indian Tribe of Wisconsin, 

Keshena, Wisconsin (Wis.)

Jennifer Hill-Kelley, Former Environmental Services Director, Menominee Indian Tribe of Wisconsin, 

Keshena, Wis.

Kevin Staus, Field Engineer, Indian Health Service, Rhinelander, Wis.

Field and Model Development Support

Leah Kammel, Hydrologist, U.S. Geological Survey Wisconsin Water Science Center, 

Middleton, Wis.

Greg Mueller, Former Hydrologist, U.S. Geological Survey Wisconsin Water Science Center, 

Middleton, Wis.

James Rauman, Hydrologic Technician, U.S. Geological Survey Wisconsin Water Science 

Center, Middleton, Wis.

Jason Smith, Hydrologic Technician, U.S. Geological Survey Wisconsin Water Science Center, 

Middleton, Wis.

Publishing Support

Michelle Greenwood, Reports Specialist, U.S. Geological Survey Wisconsin Water Science Center, 

Middleton, Wis.

Technical Reviewers

Christopher Hoard, Hydrologist, U.S. Geological Survey Michigan Water Science Center, 

Lansing, Michigan.

James Stark, Director, U.S. Geological Survey Minnesota Water Science Center, 

Mounds View, Minnesota.

Editor

Mike Eberle, Technical Writer-Editor, U.S. Geological Survey Columbus Publishing Service Center, 

Columbus, Ohio

Groundwater Specialist Reviewer

Rodney Sheets, Groundwater Specialist, U.S. Geological Survey, Columbus, Ohio

Approving Official

Kevin Breen, Bureau Approving Official, U.S. Geological Survey Office of Science Quality and 

Integrity, New Cumberland, Pennsylvania.
Simulation of the Regional Groundwater-Flow System 

of the Menominee Indian Reservation, Wisconsin

By Paul F. Juckem and Charles P. Dunning

Abstract

A regional, two-dimensional, steady-state groundwaterflow model was developed to simulate the groundwater-flow 

system and groundwater/surface-water interactions within the 

Menominee Indian Reservation. The model was developed by 

the U.S. Geological Survey (USGS), in cooperation with the 

Menominee Indian Tribe of Wisconsin, to contribute to the 

fundamental understanding of the region’s hydrogeology. The 

objectives of the regional model were to improve understanding of the groundwater-flow system, including groundwater/

surface-water interactions, and to develop a tool suitable for 

evaluating the effects of potential regional water-management 

programs. The computer code GFLOW was used because 

of the ease with which the model can simulate groundwater/

surface-water interactions, provide a framework for simulating 

regional groundwater-flow systems, and be refined in a stepwise fashion to incorporate new data and simulate groundwater-flow patterns at multiple scales. Simulations made with the 

regional model reproduce groundwater levels and stream base 

flows representative of recent conditions (1970–2013) and 

illustrate groundwater-flow patterns with maps of (1) the simulated water table and groundwater-flow directions, (2) probabilistic areas contributing recharge to high-capacity pumped 

wells, and (3) estimation of the extent of infiltrated wastewater 

from treatment lagoons.

The groundwater-flow model described in this report 

simulates the major hydrogeologic features of the modeled 

area, including surficial unconsolidated aquifers, groundwater/

surface-water interactions, and groundwater withdrawals from 

existing high-capacity production wells. Areas contributing 

recharge to pumped high-capacity wells on the Menominee 

Indian Reservation were delineated by tracking simulated 

water particles from the water table to wells in combination 

with Monte Carlo techniques, and maps of the probability of 

capture for each well nest were produced. Groundwater-agebased areas contributing recharge to wells were simulated 

by using the calibrated set of parameters and porosity values 

adjusted to account for bias in simulated saturated thickness. 

Simulations were performed for current (2013) pumping rates. 

The simulations show a range in sensitivity of the simulated 

areas contributing recharge to wells given the parameters 

evaluated through the Monte Carlo analysis. The areas contributing recharge to supply wells for the villages of Zoar and 

Neopit are long and narrow, with a sharp gradation from high 

to low probability of capture. The areas contributing recharge 

to supply wells for Middle Village and the village of Keshena 

exhibit a sharp gradation from high to low probability over a 

relatively small area between the well and a local groundwater 

mound. The highest probability areas contributing recharge 

to the supply wells for the Villages of Onekewat and Redwing are in the immediate vicinity of the wells. These wells 

also have an extensive area with low probability for capturing water that is likely due to a locally low hydraulic gradient and the large degree of uncertainty associated with the 

lakebed resistance parameters that control interaction between 

groundwater and local lakes. Additional field investigations 

and associated local model refinements would facilitate further 

reductions in uncertainty associated with simulated areas contributing recharge to the wells.

The likely extent of the Neopit wastewater plume was 

simulated by using the groundwater-flow model and Monte 

Carlo techniques to evaluate the sensitivity of predictive 

simulations to a range of model parameter values. Wastewater 

infiltrated from the currently operating lagoons flows predominantly south toward Tourtillotte Creek. Some of the infiltrated 

wastewater is simulated as having a low probability of flowing 

beneath Tourtillotte Creek to the nearby West Branch Wolf 

River. Results for the probable extent of the wastewater plume 

are considered to be qualitative because the method only 

considers advective flow and does not account for processes 

affecting contaminant transport in porous media. Therefore, 

results for the probable extent of the wastewater plume are 

sensitive to the number of particles used to represent flow 

from the lagoon and the resolution of a synthetic grid used for 

the analysis. Nonetheless, it is expected that the qualitative 

results may be of use for identifying potential downgradient 

areas of concern that can then be evaluated using the quantitative “area contributing recharge to wells” method or traditional 

contaminant-transport simulations. 
2 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

Introduction

A primary objective for resource managers of the 

Menominee Indian Reservation (fig. 1) is sustainable management of water resources (http://www.menominee-nsn.

gov/MITW/DepartmentDetails.aspx?departmentID=2100, 

accessed November 3, 2013). Though both surface-water and 

groundwater resources are generally plentiful, most domestic 

and municipal needs are provided by groundwater. Several 

local-scale drilling investigations (described under the heading 

“Data Sources”) were completed over a nearly 10-year period 

on the Menominee Indian Reservation, each targeting specific 

areas and resource issues. This report describes the data collected as part of those investigations and subsequent reservationwide groundwater-flow simulations.

The study described in this report was conducted by 

the U.S. Geological Survey (USGS) in cooperation with the 

Menominee Tribe of Wisconsin. The objectives of the study 

were to improve understanding of the hydrogeology of the 

Menominee Indian Reservation, estimate areas contributing recharge to pumped high-capacity production wells, and 

provide a tool for effective management of surface-water and 

groundwater resources. These objectives were achieved by 

developing a computer model using the program GFLOW to 

simulate the groundwater-flow system of the reservation. This 

computer model embodies the regional characteristics of the 

groundwater-flow system and includes additional hydrogeologic detail in areas of particular interest to answer local-scale 

questions. This model is a tool that can be used to improve the 

overall understanding of the hydrology of a region by testing alternative conceptual models of the groundwater-flow 

system. Additionally, this model can be used to identify areas 

where more hydrogeological or water-quality data are needed 

to address a specific question. Finally, this model serves as a 

foundation from which future local or site-specific models can 

be developed. The benefits of having this foundation include 

reducing construction time of future local-scale models in the 

area, reducing data-collection and interpretation efforts, and 

providing a simulated connection to the regional flow system.

On the Menominee Indian Reservation, the groundwater 

and surface-water systems are believed to be hydraulically 

well connected, and as a result the groundwater-flow model 

has been constructed to include many aspects of the surfacewater network. The simulation of groundwater flow and its 

interaction with the surface-water network is a necessary foundation for understanding and protecting the water resources of 

the reservation. By improving the understanding of the hydrology of the Menominee Indian Reservation, this study provides 

a basis for interpreting previously collected water-quality data 

and managing water resources for the future. 

The hydrogeologic framework of the Menominee Indian 

Reservation model draws from previous geologic and hydrologic studies. The bedrock geology of Wisconsin was mapped 

by Mudrey and others (1982), with a subsequent digital 

representation by Cannon and others (1997). Other reports 

present descriptions of the soil and unconsolidated glacial 

deposits of the area (Milfred and others, 1967; Hadley and 

Pelham, 1976; Farrand and others, 1984; Mickelson, 1986; 

Attig and Ham, 1999; Hooyer and Mode, 2007). The general 

hydrology of the region is described in two USGS Hydrologic 

Investigations Atlases (Olcott, 1968; Oakes and Hamilton, 

1973). Krohelski and others (1994) describe water resources 

within the Menominee Indian Reservation. 

Purpose and Scope

The purpose of this report is to describe the hydrogeology 

of the Menominee Indian Reservation and the development 

and applications of a computer model that simulates regional 

groundwater flow for the reservation. Geological and hydrological data used for this study consisted of previously published reports and interpretive maps, historical streamflow and 

water-level measurements, and hydrological and physical data 

collected during recent years by scientists of the USGS and 

Menominee Tribe. The report includes a summary of selected 

hydrogeological data, conceptualization of the hydrogeologic 

setting of the reservation, and details on the construction 

and calibration of a one-layer, steady-state, analytic element 

model that simulates groundwater flow and its interaction 

with surface-water features at a regional scale. The purpose of 

the model is to simulate regional groundwater-flow patterns, 

delineate areas contributing recharge to pumped wells, and 

estimate the extent of wastewater infiltration from a wastewater-treatment lagoon. Results of simulations conducted with 

the calibrated analytic element model are presented with maps 

and associated descriptions that illustrate (1) groundwaterflow directions, (2) areas contributing recharge to pumped 

high-capacity wells, and (3) the flow of infiltrated water from 

a wastewater-treatment lagoon.

Physical Setting

The Menominee Indian Reservation is in northeastern 

Wisconsin, and the reservation boundaries approximate those 

of Menominee County, except near Middle Village, which is 

in Shawano County. The reservation encompasses a hydrologically complex landscape of approximately 360 square miles 

(mi2

) (230,400 acres) containing numerous rivers, lakes, and 

wetlands (fig. 1). The Wolf River flows from north to south 

through the reservation and, with its tributaries, drains most of 

the reservation. The eastern part of the reservation is drained 

by the South Branch of the Oconto River as it flows eastward 

into Oconto County. Another important hydrological feature 

is the area of lakes east of Keshena, which includes the largest 

of the 44 lakes on the reservation—Legend Lake (1,304 acres). 

The majority of the other lakes (32) on the reservation are less 

than 50 acres in size; the smallest named lake is Red Springs 

Lake (1 acre).

The Menominee Indian Reservation is in part of two 

different physiographic regions of Wisconsin (Martin, 1965). 

The northwestern two-thirds is in the Northern Highlands 
Introduction 3

89° 88°50' 88°40' 88°30'

45°10'

45°

44°50'

Base map from U.S. Census Bureau, 2013, and Wisconsin

Department of Natural Resources, 2009

0 3 6 MILES

0 3 6 KILOMETERS

Menominee Indian

Reservation

WISCONSIN

EXPLANATION

Reservation boundary

Lake or pond

River or stream

State highways

County and town roads

Production well

Extent of conceptual model 

in figure 3

OCONTO COUNTY

SHAWANO COUNTY

LANGLADE COUNTY

MENOMINEE COUNTY

Wolf

River

Legend Lake Round

Lake

Long

Lake

Sand

Lake

MENOMINEE INDIAN RESERVATION

Keshena

Middle

Village

Neopit

Onekewat

Zoar

Red

Wing

West Branch

Wolf River

Tourtillotte

Creek

South

Branch

Oconto River

A

A'

A

Figure 1. Hydrologic features of the Menominee Indian Reservation, Wisconsin.
4 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

region, which is characterized by mixed glacial till of the 

most recent glacial epoch. The drainage system in this area 

is characterized by abundant rapids in the streams along with 

large undrained interstream areas where lakes and wetlands 

are common. The southeastern one-third is in the Central Plain 

region, where glacial deposits are typically sandier. The northeastern part of Wisconsin, which includes the Menominee 

Indian Reservation, receives about 30 inches (in.) of precipitation annually (Wisconsin State Climatology Office, 2013). 

Daily high temperatures typically range from an average 

of around 24 degrees Fahrenheit (°F) in January to an average of around 81 °F in July (Midwestern Regional Climate 

Center, 2013). Land cover in the reservation is dominated by 

roughly 223,500 acres of heavily forested lands (97 percent of 

total land area), representing the largest single tract of virgin 

timberland in Wisconsin. The remaining 3 percent of land area 

consists of developed land and open water (Wisconsin Department of Natural Resources, WISCLAND, 1998). 

The Menominee Indian Reservation is underlain by Precambrian crystalline rocks known as the Wolf River Granite 

(Krohelski and others, 1994). The elevation of the Precambrian crystalline bedrock surface varies by hundreds of feet 

over the study area as a result of erosion by late Pleistocene 

glaciers, but in general the surface slopes irregularly to the 

southeast at about 26 feet per mile. Bedrock elevation is about 

1,250 feet (ft) in the northwestern corner of the reservation and 

decreases to less than 700 ft in the southeastern corner (Krohelski and others, 1994). The crystalline bedrock crops out 

locally on the reservation but is most often overlain by unconsolidated glacial deposits up to about 200 ft thick. The sediments which compose these deposits were carried by glacial 

ice and meltwater across most of the reservation area during 

the Wisconsin glaciation (10,000–25,000 years ago) and are 

generally pitted outwash and moraine deposits comprised of 

sand, gravel, and till. The distribution of glacial deposits, soils, 

and vegetation was mapped by Milfred and others (1967), and 

an interpretation of the glacial sediments is provided by Farrand and others (1984) (fig. 2).

Data Sources

Background geologic information used in this study 

was based primarily upon a Quaternary geologic map of 

the Lake Superior quadrangle, United States and Canada, 

by Farrand and others (1984). Supplemental mapping and 

interpretation were garnered from maps by Milfred and others (1967), Hadley and Pelham (1976), Mickelson (1986), 

Attig and Ham (1999), and Hooyer and Mode (2007).

Test borings were drilled at three locations between 

2001 and 2008 on the Menominee Indian Reservation to 

gather localized data on depth to bedrock, static water level, 

and lithology. Borehole investigations were conducted in 

the vicinity of the villages of Neopit and Zoar from July 9 to 

July 17, 2001. A similar investigation was conducted in the 

vicinity of Keshena from August 7 to September 11, 2002. The 

Neopit area was revisited for additional borehole investigations from May 29 to July 1, 2008. Data were collected from 

boreholes augered with either a CME75 drill rig or Geoprobe 

66DT. Field operations included the following steps: (1) set 

up the drill rig or Geoprobe and determine map coordinates 

by using a handheld GPS unit, (2) auger until bedrock refusal 

or until all available augers had been added to the auger string, 

(3) collect and roughly describe the lithology of auger returns 

(Zoar and Keshena only), (4) temporarily install a monitoring 

well in the borehole by using a 5-ft-long slotted screen placed 

so as to intersect the water table, and (5) measure and record 

the static water level in the monitoring well upon stabilization. 

For boreholes in which a monitoring well was not installed, 

the borehole was permanently abandoned by backfilling with 

bentonite grout chips to about 5 ft of land surface, followed 

by emplacement of auger cuttings to the surface. These data 

increased the understanding of the buried topography of the 

crystalline bedrock and the thickness and saturated thickness 

of the unconsolidated aquifer at select locations. Details on 

these borehole investigations are provided in appendix 1.

River water-level elevations derived from USGS topographic maps were used for model development and are 

expected to represent long-term average base-flow conditions 

to within a few feet. Although water levels on topographic 

maps represent elevations observed during the mapmaking period, the rivers are assumed to be at a geomorphic 

steady state in terms of erosion and deposition, with little 

systematic trend in river base flow water level over time. All 

surface-water-level measurements are reported as feet and 

are referenced to the North American Vertical Datum of 1988 

(NAVD 88). River flow data consist of gaged and miscellaneous discharge measurements archived in the USGS National 

Water Inventory System (NWIS) database (Dempster, 1990) 

and analyzed with methods developed by Gebert and others 

(2011) to estimate base flows. Groundwater level data were 

obtained from NWIS, well-construction reports submitted to 

the Wisconsin Department of Natural Resources (2010), and 

measurements made during this study (appendix 1). 

Modeling Methods

An analytic element groundwater-flow model was 

developed to simulate the regional groundwater-flow system 

and its interaction with surface-water features; the computer 

program GFLOW (Haitjema, 1995) was used to construct the 

model. The model simulates groundwater flow by using the 

Dupuit-Forchheimer approximations, which simplify a threedimensional flow system into a two-dimensional, areal flow 

system. These approximations are well suited for the Menominee Indian Reservation, where the ratio of horizontal to vertical dimensions of the aquifer is very large (Haitjema, 1995); 

that is, the Menominee Indian Reservation extends a distance 

from east to west of about 25 mi, or 132,000 ft, whereas the 

total aquifer thickness does not exceed about 200 ft. A complete description of analytic element modeling is beyond the 
Data Sources 5

89°20' 89° 88°40' 88°20'

45°20'

45°

44°40'

Base map compiled from digital data sources 

and Farrand and others (1984)

0 8 16 MILES

0 8 16 KILOMETERS EXPLANATION

Outwash sand and gravel

Ice-contact sand and gravel

Lake (open water)

Lake silt and clay

Clayey till

Calcareous sandy loamy till

Noncalcareous sandy loamy till

Sandy loamy till

Loamy till

Loamy till (pre-Wisconsinan)

Glacial deposits

Reservation boundary

County boundary

(Lithologic code from Farrand and others, 1984)

gg

kg

lca

tc

td

tdb

tde

tl

tlq

Figure 2. The distribution of glacial deposits in the Menominee Indian Reservation, Wisconsin, after Farrand and others (1984).
6 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

scope of this report, but a brief description is provided in the 

following paragraph. In addition, Hunt (2006) gives a review 

of applications of the analytic element method, and Haitjema 

(1995) discusses the underlying concepts and mathematics of 

the method in detail. 

A horizontally infinite aquifer is assumed in analytic 

element modeling. The problem domain (model area) does 

not require a grid or involve interpolation between cells as 

required for finite difference or finite element methods. To 

construct an analytic element model, the modeler identifies 

features that are important for controlling groundwater flow 

(for example, wells, surface-water features, and zones of 

uniform transmissivity) and represents them in the model as a 

point element or line element with defined physical properties. 

The closer a hydrologic feature is to an area of interest, the 

greater the amount of detail required to adequately simulate 

that feature. For example, greater detail for a river would 

result in a linesink string with many individual segments that 

follow the meandering channel or variations in downstream 

slope. Analytic elements included in the Menominee Indian 

Reservation model included specified-discharge wells, linesinks (river), and inhomogeneities (zones of transmissivity). 

Each analytic element adds an equation to the model, and the 

effects of these individual features are added together to form 

a solution for any location in the simulated groundwater-flow 

system. Because the solution is not confined to a grid, hydraulic heads and hydrologic flows can be computed at any location in the model domain without interpolating between grid 

cells. In the GFLOW model developed for the Menominee 

Indian Reservation, the analytic elements are two dimensional 

and only simulate steady-state conditions (water levels do not 

vary with time). A primary value of regional-scale analytic element modeling is to identify the hydrological features which 

are most important in controlling groundwater flow. These 

features can then be used to better define local conditions, test 

hypotheses, and answer site-specific questions as specific local 

data are incorporated into the model. The analytic element 

method and comparisons of analytic element to finite-difference numerical modeling techniques have been discussed by 

others (Haitjema, 1995; Hunt and others 1998; and Hunt and 

others, 2003).

The GFLOW model of the Menominee Indian Reservation was calibrated by means of parameter-estimation techniques. Numerous publications detail the advantages of parameter estimation (for example, Poeter and Hill, 1997; Kelson 

and others, 2002). Briefly, the primary benefit of a properly 

prepared and executed parameter-estimation calibration over 

typical trial-and-error calibration is the ability to automatically 

calculate parameter values (for example, hydraulic conductivity and recharge) that are a quantified best fit between simulated model output and observed data (typically groundwater 

level and streamflow data). In addition, parameter sensitivity can be quantified to assess model parameters that have 

the greatest influence on simulated results. In this study, the 

GFLOW model was coupled with the parameter estimation 

code PEST (Doherty, 2012).

Conceptual Model of the GroundwaterFlow System

Conceptualization of the hydrologic system forms the 

framework for mathematical model development and simplifies the groundwater system into important components. 

Simplification of the natural system is necessary because 

inclusion of all its complexities into a computer model is not 

feasible given existing knowledge of the subsurface. Steps 

in the development of the conceptual model are (1) characterization of aquifers and confining units, (2) identification 

of sources and sinks of water, and (3) identification and 

delineation of hydrologic boundaries encompassing the area 

of interest. The first two of these steps were accomplished by 

review and interpretation of available geologic and hydrogeologic data. The third step was accomplished through the 

model design (see “Hydrogeologic Boundaries,” below). The 

conceptual model of the groundwater-flow system is shown in 

figure 3. The vertical scale of the conceptual model diagram 

has been greatly exaggerated to illustrate the geologic relationships (fig. 3A); regional flow through the aquifers is primarily 

horizontal (fig. 3B).

Aquifers and Confining Units

The principal aquifer of the Menominee Indian Reservation is composed of laterally extensive unconsolidated glacial 

deposits, represented generally by pitted outwash and moraine 

deposits. Outwash deposits are composed of stratified sand 

and gravel that occur as layers, lenses, terrace deposits, and 

valley fills that range in thickness across the study area from 

0 to 200 ft (Krohelski, and others, 1994). Moraine deposits 

are primarily composed of poorly sorted tills that incorporate 

lithologies ranging from clay and silt to sand, gravel, and 

cobbles, although most tills in the study area are mapped as 

sandy tills (Hadley and Pelham, 1976; Farrand and others, 

1984; Mickelson, 1986; Attig and Ham, 1999; Hooyer and 

Mode, 2007). Confining units of finer grained glacial till are 

believed to be only locally important.

The reservation is underlain by crystalline bedrock composed of the Wolf River Granite (Krohelski, and others, 1994). 

Although the crystalline bedrock can be a productive aquifer 

and is often included in the completion interval of wells, its 

hydraulic conductivity is dependent largely on the presence of 

fractures and is much less than that of the overlying unconsolidated sediments (Krohelski, and others, 1994). Therefore, 

the highly transmissive glacial sediments are considered to be 

bounded below by low-transmissivity bedrock in the model. 

However, where present, fractures in the upper crystalline 

bedrock may contribute locally to the groundwater resource.

The irregularity of the bedrock surface and the hilly 

terrain across the Menominee Indian Reservation result in a 

glacial aquifer with variable depths to groundwater and saturated thickness across the study area (Krohelski, and others, 

1994). The glacial aquifer has a saturated thickness ranging 
Conceptual Model of the Groundwater-flow System 7

Vertical dimension is greatly exaggerated.

Vertical dimension is greatly exaggerated.

Legend Lake 

area

Unsaturated

Zone

Glacial aquifer

Crystalline bedrock

Recharge

General direction of groundwater flow

B. Conceptual model

Wolf River Legend Lake 

area

Unsaturated

Zone

Wolf River

West Branch 

SOUTHEAST

A. Idealized cross section

Crystalline bedrock

EXPLANATION

NORTHWEST

A A'

SOUTHEAST

A'

NORTHWEST

A

Wolf River

Wolf River

West Branch 

3 2 6 1

ZONE ZONE ZONE ZONE

Glacial aquifer

Figure 3. The groundwater-flow system of the Menominee Indian Reservation, Wisconsin. A, Idealized cross section. 

B, Conceptual model developed to guide construction of the groundwater-flow model. GFLOW, the groundwater-flow modeling 

code, requires a horizontal aquifer base (top of crystalline bedrock). Zone numbers in figure 3B correspond with zone numbers 

in figure 4.
8 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

from zero (in areas near Precambrian bedrock outcrops) to 

approximately 150 ft. The regional groundwater-flow direction generally matches that of the rivers, flowing from the 

northwest toward the southeast as more complex local flow 

patterns interact with local rivers and lakes (plate 3 of Krohelski and others, 1994). 

Sources and Sinks of Water

The primary source of water to the groundwater-flow 

system on the Menominee Indian Reservation is precipitation that infiltrates through the ground to the water table—a 

process known as recharge and expressed as a rate in inches 

per year. Groundwater moves from higher to lower hydraulic 

potential (areas of higher groundwater levels to areas of lower 

groundwater levels). As a result, water generally enters the 

groundwater system areally throughout the reservation and 

discharges to surface-water features. These features behave 

as groundwater sinks in the model, where groundwater leaves 

the aquifer to become surface water. Accurate locations of 

surface-water features are, therefore, critical to simulating the 

groundwater system and groundwater/surface-water interaction within the reservation.

Pumped wells are another type of hydrological sink, 

capturing groundwater that would otherwise discharge to 

rivers, lakes, and wetlands. In areas of significant groundwater withdrawal by wells, rivers that otherwise would receive 

groundwater as base flow may instead locally recharge the 

groundwater system (Feinstein and others, 2012). Only a 

small amount of the total recharge to the water table within 

the reservation is captured by wells. However, areas with large 

annual withdrawals (for example, the city of Keshena) can 

capture a substantial portion of the local recharge. 

Hydrogeologic Boundaries

The groundwater-flow model consists of two domains: 

the near field and the far field. The near field is the area of 

interest, which for this study is the entire Menominee Indian 

Reservation. The near-field rivers and lakes are represented by 

linesink networks (Haitjema, 1995). A linesink is a mathematical representation of a hydrologic sink (such as a stream), and 

a network of such linesinks will route base flow (the portion of 

total flow derived from groundwater) downstream. The degree 

of groundwater/surface-water interaction at a particular river 

segment depends on the riverbed sediment resistance (thickness divided by vertical hydraulic conductivity) and the difference in elevation between the river stage and the adjacent 

water table. The headwaters of certain river systems surrounding the reservation were simulated as near-field linesinks in 

order to allow the simulated river to “go dry” if the water table 

is simulated as being below the riverbed elevation, which 

removes the linesink from the solution. This is an important 

consideration near groundwater divides where ephemeral 

streams tend to be prevalent.

The model far field is the area surrounding the near field 

and contains hydrologic features that control the groundwater 

flow toward or away from the near field. These features are 

rivers and lakes that border the reservation and are simulated 

with coarsely defined linesink networks having little or no 

resistance between the surface-water features and the groundwater system. The function of the far field is to establish the 

groundwater divides near the edge of the model that to a large 

extent determine the flow of water to sinks (rivers, lakes, 

wetlands, and wells) within the near field. The importance 

of properly simulated model boundaries is to ensure that 

simulated responses to stresses applied to the near field (for 

example, pumping or other scenario simulations) are not artificially affected by conditions at the boundaries of the model 

domain. An efficient way to design boundaries is to have the 

groundwater-flow model directly simulate hydraulic divides 

between the area of interest (near field) and adjacent rivers 

(Anderson and Woessner, 1992; Haitjema, 1995).

Hydraulic Properties of the 

Groundwater-Flow System

Initial estimates of aquifer hydraulic conductivity, 

recharge, and riverbed resistance for the regional groundwaterflow model were based on available geologic and hydrologic 

data. The following is a brief description of these estimates. 

Final values were derived during the parameter-estimation 

process.

Hydraulic Conductivity

Glacial lithologies within the Menominee Indian Reservation were mapped by Farrand and others (1984) as glacial 

tills (loamy and calcareous sandy loamy), outwash sand and 

gravel, and minor deposits of glacial lakebed silts and clays. 

Whereas these descriptions suggest a potentially large range 

in lithology and associated hydraulic conductivities, descriptions from reports and maps by Hooyer and Mode (2007), 

Mickelson (1986), and Hadley and Pelham (1976) describe 

relatively modest variation among most tills in the area, with 

most tills being composed primarily of sand. Moreover, glacial maps rarely consider the degree of connectivity among 

clay units, a factor that is important for making regional-scale 

hydraulic-conductivity estimates. Hydraulic conductivity of 

the surficial glacial aquifer ranged from 0.1 to 55 feet per day 

(ft/d) for 11 slug tests conducted within the reservation, as 

reported in Krohelski and others 1994. 

Pumping tests conducted on the Zoar community wells in 

the 1970s provide a further estimate for hydraulic conductivity 

of the glacial aquifer in this region of the Menominee Indian 

Reservation. The pumping rates and measured drawdowns 

for wells Zoar 1 and Zoar 2 (Menominee Tribe of Wisconsin, 

Department of Environmental Services, written commun., 
Hydraulic Properties of the Groundwater-Flow System 9

2001) were first used to calculate the specific capacity of the 

aquifer (table 1). The specific-capacity values are the basis for 

estimating hydraulic conductivity by means of the program 

TGUESS (Bradbury and Rothschild, 1985). The value of the 

aquifer storage coefficient is an input to the TGUESS program; 

a range of values (0.01 to 0.0001) was used to provide a range 

of estimated hydraulic conductivities. On the basis of these 

storage coefficients and the pumping-test data, locally calculated hydraulic conductivities range from 17.3 to 19.5 ft/d. The 

results were very similar among the three tests conducted on the 

two Zoar wells, and the hydraulic-conductivity estimates fall 

within the range of the values estimated previously from slug 

tests (Krohelski and others, 1994).

The hydraulic conductivity of the crystalline bedrock is 

estimated to be 0.3 ft/d (Krohelski, and others, 1994). This is 

between 1 and 2 orders of magnitude less than that of the glacial 

aquifer. In the model, the crystalline bedrock is simplified to 

represent an impermeable boundary even though locally it may 

have some secondary permeability. 

Although these ranges and estimates are useful for characterizing the system, the model requires specific values for 

the hydraulic conductivity across regions, or zones, in the flow 

system. Thus, on the basis of available data, values associated 

with zones of locally uniform hydraulic conductivity were 

treated as calibration parameters. These zones were delineated 

by considering multiple geologic maps and descriptions of glacial material along with topographic features. Values used in the 

modeling described here were initialized by using reasonable 

values, with the final values determined during calibration by 

use of PEST (Doherty, 2012) and constrained to a reasonable 

range based on available measurements and estimates. 

Recharge

Recharge takes place nearly everywhere on the reservation except in groundwater discharge areas associated with 

surface-water bodies. Recharge in the Menominee Indian Reservation and upstream areas ranges from about 5 to 12 inches 

per year, as determined from base flows in gaged streams 

divided by the gaged area (Gebert and others, 2011). The 

magnitude of estimated recharge can vary spatially as a result 

of the characteristics of the watershed (precipitation, lithology, 

topography, land cover and use, evapotranspiration, and so 

forth), and also from stream capture of groundwater from adjacent surface watersheds that are not coincident with groundwater basins. These local differences tend to average out over 

large areas such as the Menominee Indian Reservation and 

show little spatial trend within or upstream of the reservation. 

Thus, an average recharge rate was applied uniformly across 

the regional model and was treated as a calibration parameter, 

with the value constrained to reasonable ranges and determined by use of the parameter estimation code PEST to match 

observed water levels and base flows. 

Table 1. Pumping-test data and TGUESS calculations for production wells near Zoar, Menominee 

Indian Reservation, Wisconsin.

[(gal/min)/ft, gallons per minute per foot; ft2/d, feet squared per day; ft/d, feet per day]

Well parameter Zoar well 1 Zoar well 2

[low pump rate]

Zoar well 2

[high pump rate]

Well diameter (inches) 6.0 8.0 8.0

Total depth (feet) 86.0 82.0 82.0

Static water level (feet) 12.0 18.2 18.2

Saturated thickness (feet) 75.5 64.8 64.8

Open interval (feet) 74.0–86.0 72.0–82.0 72.0–82.0

Test date 12/04/75 05/03/77 05/03/77

Test duration (hours) 3 2 2

Pump rate (gallons per minute) 60 60 73.3

Final depth to water (feet) 47.8 55.9 65.7

Well loss (assumed) 0.8 0.8 0.8

Storage coefficient (range) 0.01–0.0001 0.01–0.0001 0.01–0.0001

Specific capacity [(gal/min)/ft)] 1.68 1.59 1.54

Transmissivity (ft2

/d) 1,339–1,460 1,149–1,270 1,115–1227

Hydraulic conductivity (ft/d)* 17.7–19.3 17.8–19.5 17.3–18.9

*

 Calculated by TGUESS.
10 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

Riverbed Resistance

Estimates of riverbed and lakebed resistance are needed 

to simulate the interaction between surface water and groundwater. Bottom resistance for a water body is equal to the 

thickness of a riverbed or lakebed divided by its vertical 

hydraulic conductivity. In this model, riverbed resistance 

was set at 0.5 day. (The unit “day” is a mathematical reduction resulting from thickness, in feet, being divided by the 

hydraulic conductivity, in feet per day.) Resistance was set 

at 1 day for lakes that outlet to rivers and 10 days for lakes 

with no surface inflow/outflow. These values correspond to a 

1-ft sediment thickness and vertical hydraulic conductivities 

of 2, 1, and 0.1 ft/d, respectively, and are within commonly 

observed ranges for bed sediment (Calver, 2001). Because 

of the relative insensitivity of the regional flow system to 

riverbed resistance (see “Sensitivity Analysis” section), this 

parameter was not adjusted during calibration but was evaluated during the Monte Carlo analysis. 

Groundwater Withdrawals

Groundwater withdrawals from the glacial aquifer were 

incorporated into the GFLOW model according to the conceptual model design. Pumping rates for high-capacity production 

wells were computed as the annual average rate provided by 

Indian Health Service (Kevin Staus, written commun., 2013). 

Pumping from small, individual private wells in the reservation is not included in the model because the discharge from 

these wells is widely distributed and relatively small (especially when including return flow from septic systems), so it 

has a negligible effect on the overall regional water table.

Total high-capacity well withdrawal in the Menominee 

Indian Reservation is about 117.6 million gallons per year 

(Mgal/yr) distributed across 13 different wells on the reservation. The largest municipal withdrawals in the model were in 

Keshena (60.3 Mgal/yr) and Neopit (32.9 Mgal/yr). These and 

other production wells are listed in table 2. The groundwaterflow model included three production wells from outside of 

the reservation—the Dumke C wells of the Stockbridge Munsee Tribe, which withdrew water at about 27.4 Mgal/yr total. 

Simulation of the Regional 

Groundwater-Flow System

An analytic element groundwater-flow model of the 

Menominee Indian Reservation was developed by using the 

computer program GFLOW (Haitjema, 1995). The model 

simulates the groundwater-flow system and its interaction with 

surface-water features. The model consists of one layer and 

simulates steady-state conditions (no change in water levels 

over time). Simulated rivers and lakes, and zones of hydraulic 

conductivity are shown in figure 4.

Table 2. Groundwater withdrawal rates for production wells simulated in the groundwater-flow model of the Menominee 

Indian Reservation, Wisconsin.

[gal/yr, gallons per year; gal/d, gallons per day; ft3/d, cubic feet per day]

Well Name Pumping rate

(gal/yr)

Pumping rate

(gal/d)

Pumping rate

(ft3

/d)

Total pumping 

by community

(gal/yr)

Keshena 3 20,088,750 55,000 7,352

Keshena 4 20,088,750 55,000 7,352 60,266,250

Keshena 5 20,088,750 55,000 7,352

Zoar 3 2,995,050 8,200 1,096

Zoar 4 2,995,050 8,200 1,096 8,985,150

Zoar 5 2,995,050 8,200 1,096

Neopit 3 10,957,500 30,000 4,010

Neopit 4 10,957,500 30,000 4,010 32,872,500

Neopit 5 10,957,500 30,000 4,010

Middle Village well 7,852,875 21,500 2,874 7,852,875

Redwing well 1 2,483,700 6,800 909

4,967,400 Redwing well 3 2,483,700 6,800 909

Onekewat well 2,629,800 7,200 963 2,629,800

Stock-Bridge Dumke C well 1 9,131,250 25,000 3,342

Stock-Bridge Dumke C well 2 9,131,250 25,000 3,342 27,393,750

Stock-Bridge Dumke C well 3 9,131,250 25,000 3,342
Simulation of the Regional Groundwater-Flow System 11

89°20' 89° 88°40' 88°20'

45°20'

45°

44°40'

Zone with base elevation 

of 1,100 feet and hydraulic 

conductivity parameter 

value of 4.6 feet per day

Base map compiled from digital data sources

6

1

2

3

5 4

7

County boundary

400 ft, 0.5 ft/d, 4.4 ft/d

XXX ft, X.X ft/d, X.X ft/d

400 ft; 8 ft/d, 21.3 ft/d

500 ft; 2.9 ft/d, 21.7 ft/d

600 ft; 4.6 ft/d, 21.5 ft/d

700 ft; 18.5 ft/d, 64.9 ft/d

700 ft; 12.3 ft/d, 43.1 ft/d

1,100 ft; 4.6 ft/d, 16.1 ft/d

Nonrouted near-field

 linesink

Routed near-field 

 linesink

Far-field linesink

EXPLANATION

Base elevation (in feet above NAVD 88)

 Hydraulic conductivity parameter (in feet per day)

Zone Effective hydraulic conductivity estimate (in feet per day)

number

X EXAMPLE

6

1

2

3

5

4

7

0 8 16 MILES

0 8 16 KILOMETERS

Figure 4. Hydrologic features and base elevation and hydraulic-conductivity zones simulated with analytic elements, 

Menominee Indian Reservation, Wisconsin.
12 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

Model Construction and Assumptions

Initial model development included estimating the 

elevation of the base of the groundwater system, a recharge 

rate, and a horizontal hydraulic conductivity. The base of the 

model is set at 400 ft above NAVD 88 for the southeastern 

part of the reservation (zone 1). The model base increased to 

500 ft for the central part (zone 2) and 600 ft for northwestern 

part of the reservation (zone 3), representing the increasing 

elevation of the crystalline bedrock. Smaller areas that were 

expected to have different geologic lithologies were separated into individual zones (zones 4, 5, and 6) but maintained 

the simulated bedrock elevation from the surrounding area. 

A local bedrock high near the village of Zoar (zone 7) that was 

identified through borings was assigned a base elevation of 

1,100 ft above NAVD88. The regional base elevations are hundreds of feet below the mapped bedrock elevations in the area 

(Krohelski and others, 1994) but generally match the regional 

bedrock dip. Lower-than-actual base elevations were used as 

a means to improve performance of the model and stability of 

the conjunctive groundwater and surface-water solution for 

the model; hydraulic-conductivity adjustments were applied to 

ensure proper representation of the bulk aquifer transmissivity; 

that is, the combination of zones with flat bases (a requirement 

of analytic element methods) and the sloping streams resulted 

in thin saturated thicknesses in the model when more realistic 

base elevations were used. This thin saturated thickness amplifies problems associated with nonlinearity of the conjunctive 

groundwater and surface-water solution. Beyond adjusting 

solver settings (such as numerical relaxation), increasing the 

simulated saturated thickness in the model by using low base 

elevations reduced nonlinearity problems and stabilized the 

model solution, resulting in an improved simulation.

In two-dimensional areal models, where transmissivity 

(hydraulic conductivity multiplied by aquifer saturated thickness) of a single layer represents the flow system, the base 

elevation is correlated with hydraulic conductivity; that is, the 

same aquifer transmissivity can be represented with a thick 

aquifer and low hydraulic conductivity values as could be 

represented by a thin aquifer and high conductivities. Therefore, parameter calibration focused on hydraulic conductivity 

in each zone rather than the aquifer base elevations. Use of 

lower-than-actual base elevations, or overestimating saturated 

thickness, resulted in estimates of hydraulic-conductivity 

parameter values that were lower than computed hydraulic 

conductivities from pumping tests. As a result, the calibrated 

hydraulic-conductivity parameter values listed in table 3 are 

only reasonable for the base elevations used in this model. 

More realistic effective hydraulic conductivity estimates for 

simulated zones in the model are listed in table 4. These calculated effective hydraulic-conductivity values were computed 

by adjusting for more realistic saturated thicknesses based on 

mapped bedrock elevations. The effective hydraulic-conductivity and adjusted saturated-thickness values generally match 

values computed from pumping-test results near the village of 

Zoar (table 1).

This representation of transmissivity has two noteworthy implications for model simulations. First, the range in 

simulated transmissivity within an inhomogeneity is less than 

would have been simulated by using shallow base elevations. This result is viewed as a positive attribute of the model 

because the actual bedrock surface is sloped, and not horizontal with discrete step changes in elevation, as must be simulated with analytic elements; that is, the relatively more subtle 

simulated changes in transmissivity within and across zones 

is likely more representative of the natural system than would 

have been simulated with zones of very shallow bedrock surfaces. Second, simulated groundwater traveltimes are affected 

by saturated thickness (for example, Vogel, 1967; Haitjema, 

1995). That is, groundwater velocity is computed according to 

equation 1, as modified from Anderson and Woessner (1992):

 (1)

dl

dh

ne neb

K

v = − dl

T dh = −

where

v is the velocity of groundwater flowing 

through connected pore spaces,

K is the aquifer hydraulic conductivity,

ne is the aquifer effective porosity (the 

porosity available for groundwater 

movement), which is referred to simply as 

porosity in this report,

dh/dl is the hydraulic gradient, or change in water level 

divided by change in distance,

b is the saturated thickness of the aquifer, and

T is the aquifer transmissivity, which is 

equal to K*b.

 Similar to the adjustment of hydraulic conductivity as a 

means to maintain reasonable aquifer transmissivities in the 

model due to the exaggerated saturated thickness simulated in 

the model, porosity was locally adjusted to maintain realistic groundwater velocity and traveltime simulations around 

pumped wells. For example, doubling the simulated saturated 

thickness would require halving the aquifer porosity in order 

to produce an “adjusted” porosity value for use in the model 

that would maintain appropriately simulated traveltimes. The 

estimated aquifer porosity (0.2) and corresponding adjusted 

porosity values used in the model are provided in table 5.

The river network for the Menominee Indian Reservation 

(fig. 1) is represented in the GFLOW model as a series of linesinks (fig. 4). Multiple linesinks are joined into linesink strings 

representing river segments with uniform downstream slope 

(change in water-surface elevation over distance). The river 

slope assigned to linesink strings in the Menominee Indian 

Reservation model was based on data from 7.5-minute topographic quadrangle maps. GFLOW solves for the exchange 

between groundwater and surface water at the center of each 

linesink. 
Simulation of the Regional Groundwater-Flow System 13

Table 3. Calibrated and specified parameter values and composite sensitivities for the groundwater-flow model of the Menominee 

Indian Reservation, Wisconsin.

[in/yr, inches per year; ft/d, feet per day]

Parameter name Calibrated or specified

parameter value

Sensitivity 

(no units)

Description 

(from Farrand and others, 1984)

R_reservation 8.4 in/yr; calibrated 0.08 Uniform recharge to the entire model area.

K_southeast 8.0 ft/d; calibrated 0.13 Hydraulic-conductivity parameter for the southeastern 

part of the model area, representing sandy loamy till 

and outwash sand and gravel glacial deposits.

K_central 2.9 ft/d; calibrated 0.07 Hydraulic-conductivity parameter for the central part 

of the model area, representing sandy, loamy till and 

outwash sand and gravel glacial deposits.

K_northwest 4.6 ft/d; calibrated 0.12 Hydraulic-conductivity parameter for the northwestern 

part of the model area, representing sandy, loamy till 

and outwash sand and gravel glacial deposits.

K_antigo_outwash_mult 1.5; calibrated 0.006 Multiple of the K_sandy_till parameter representing 

outwash sand and gravel in the Antigo area of the 

model.

K_sandy_till 12.3 ft/d; calibrated 0.07 Hydraulic-conductivity parameter for the sandy till 

deposit in the northwestern part of the model area 

east of Antigo.

K_loamy_till 0.5 ft/d; calibrated 0.07 Hydraulic-conductivity parameter for loamy till deposits in the southeastern part of the model area.

C_drainage_lakes 1.0 day; specified 0.003 Lakebed resistance for drainage lakes (have an outlet 

stream) simulated with routed linesinks.

C_seepage_lakes 10 days; specified 0.005 Lakebed resistance for seepage lakes (no outlet stream) 

simulated with nonrouted linesinks.

C_rivers 0.5 day; specified 0.02 Riverbed resistance for perennial rivers simulated with 

routed and nonrouted near-field linesinks.

Table 4. Calibrated hydraulic-conductivity parameters and associated base elevations used to calculate effective hydraulic-conductivity 

estimates for the Menominee Indian Reservation, Wisconsin.

[ft, feet; ft/d, feet per day; NAVD 88, North American Vertical Datum of 1988]

Parameter 

name

Estimated 

actual 

base 

elevation

(ft above 

NAVD 88)1

Base 

elevation 

in

GFLOW 

(ft above 

NAVD 88)

Estimated 

average water 

level from the 

calibrated 

GFLOW model 

(ft above 

NAVD 88)

Simulated 

saturated

thickness 

(ft)

Calculated 

actual 

saturated 

thickness 

(ft)

Calibrated 

hydraulicconductivity 

parameter 

(ft/d)

Calculated 

effective 

hydraulic 

conductivity 

(ft/d)

Geologic 

description 

(from Farrand 

and others, 

1984)

K_southeast 700 400 880 480 180 8.0 21.3 Mix of outwash 

and till

K_central 950 500 1,020 520 70 2.9 21.7 Mix of outwash 

and till

K_northwest 1,150 600 1,300 700 150 4.6 21.5 Mix of outwash 

and till

K_shallow_base 1,150 1,100 1,170 70 20 4.6 16.1 Mix of outwash 

and till

K_antigo_outwash 1,250 600 1,510 910 260 18.5 64.9 Sand and gravel 

outwash

K_sandy_till 1,250 600 1,510 910 260 12.3 43.1 Sandy till

K_loamy_till 850 400 910 510 60 0.5 4.4 Loamy till

1

Actual base elevation for each zone was estimated from Krohelski and others (1994).
14 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

Table 5. Adjusted porosity estimates used to simulate time-dependent areas contributing recharge for pumped wells, Menominee 

Indian Reservation, Wisconsin. Associated values for estimated and simulated bedrock elevations and saturated thicknesses are 

included for reference. Effective porosity of the aquifer was estimated at approximately 0.2, which was locally adjusted within the model 

simulations to account for simulated biases in the simulated saturated thickness. (Actual base elevations were estimated from wellconstruction reports as available for individual wells. Base elevation was estimated from Krohelski and others, 1994, for wells that did 

not extend to bedrock.)

[ft, feet; ft/d, feet per day]

Simulated 

area

Estimated 

actual 

base 

elevation

(ft above 

NAVD 88)

Base 

elevation 

in

GFLOW 

(ft above 

NAVD 88)

Estimated 

average water 

level from the 

calibrated 

GFLOW model 

(ft above 

NAVD 88)

Simulated 

saturated

thickness 

(ft)

Calculated 

actual 

saturated 

thickness (ft)

Calibrated 

hydraulicconductivity 

parameter 

(ft/d)

Calculated 

effective 

hydraulic 

conductivity 

(ft/d)

Adjusted 

porosity 

used in 

GFLOW (no 

units)

Redwing wells 700 400 845 445 145 8.0 24.6 0.065

Onekewat wells 725 400 846 446 121 8.0 29.5 0.054

Keshena wells 780 400 850 450 70 8.0 51.4 0.031

Keshena upgradient 

area

800 400 905 505 105 8.0 38.5 0.042

Middle Village 910 400 994 594 84 0.5 3.6 0.028

Neopit wells 990 600 1,063 463 73 4.6 29.2 0.032

Neopit immediately 

upgradient area

1,043 600 1,095 495 52 4.6 43.8 0.021

Neopit distant 

upgradient area

1,100 600 1,145 545 45 4.6 55.8 0.017

Zoar wells 1,125 600 1,182 582 57 4.6 47.1 0.020

Zoar upgradient area 1,175 600 1,240 640 65 4.6 45.4 0.020

In the near field of the model—that is, within and immediately adjacent to the reservation—each linesink is assigned 

a width based on stream order and field observations; widths 

ranged from about 5 ft for headwaters to 300 ft along the Wolf 

River near Keshena. Each near-field linesink is also assigned 

a resistance term that is equal to the thickness of the riverbed 

divided by its hydraulic conductivity. This resistance term is 

multiplied by the river width and the difference between the 

fixed river level and calculated water-table elevation adjacent 

to the river to compute groundwater flow across the riverbed. 

In the regional model, a single value of resistance equal to 

0.5 day was applied to all rivers. Initial parameter sensitivities demonstrated that the model results were not sensitive to 

changes in riverbed resistance when varied over reasonable 

ranges; therefore, the values for all rivers were fixed for all 

model runs. 

Near-field linesinks are linked so that base flow is routed 

from headwaters at higher elevations through tributaries to the 

main trunk of rivers at lower elevations. By routing base flow 

through the river network, the amount of water gained from 

or lost to the aquifer system by the river is tabulated. This 

accounting allows the amount of water simulated in the river 

at any point to be compared to flows recorded at streamgages. 

In general, streamflow consists of (1) overland flow derived 

mostly from storms and (2) base flow derived from groundwater discharge. Only the base-flow component of streamflow is 

simulated with the GFLOW model. 

Drainage lakes, or lakes that have a surface-water outlet, 

were simulated as linesinks with resistance along their shorelines and assigned widths based on the approximate distance 

from shoreline to shoreline (Haitjema, 2005) and a uniform 

resistance of 1 day. Seepage lakes, or lakes that have no surface inlet or outlet, in areas of interest near production wells 

(Round, Sand, and Long Lakes near the Villages of Onekewat and Redwing) also were simulated with linesinks along 

their shorelines. These seepage lakes were similarly assigned 

widths based on the approximate distance from shore to shore 

and assigned a uniform resistance of 10 days. Simulating 

seepage lakes with linesinks allowed for uncertainty associated with their relative degree of connection with the water 

table to be evaluated during Monte Carlo uncertainty simulations. Seepage lakes outside of local areas of interest that were 

expected to be well connected with the groundwater system 

were represented in the model as water-level targets.
Simulation of the Regional Groundwater-Flow System 15

Linesinks also represent water bodies in the model “far 

field” area outside of the reservation (fig. 4). However, these 

far-field elements are assigned no resistance or width. The 

assigned stage for each far-field linesink is equivalent to the 

water-table elevation along the linesink. The far-field linesinks, 

therefore, act as fixed water-level conditions that serve as 

important groundwater sinks or sources outside the near field. 

In this manner, the far-field water bodies help to define the 

groundwater divides around the outer perimeter of the model. 

This design ensures that simulated responses to stresses 

applied to the near field are not artificially affected by conditions at the boundaries of the model domain. 

Other inputs to the GFLOW model include pumping wells and recharge zones. Wells are assumed to be fully 

penetrating (screened from the water table to the model base) 

and to have constant pumping rates, as described previously 

(see the “Groundwater Withdrawals” section). Recharge was 

simulated uniformly across the model. Although some local 

variation in recharge is expected, such variation “averages 

out” at the regional scale of the model (see “Recharge” section 

above). Additional recharge was added to an area in the model 

representing the seepage cells of the wastewater-treatment 

facility near the community of Neopit. The steady-state rate 

of wastewater added to the facility is 86,500 gallons per day 

(gal/d), as reported by Indian Health Service (Kevin Staus, 

written commun., 2013), and the model adds this as recharge 

to the water table.

Model Calibration

Ground-water model calibration is a process whereby 

simulated values of groundwater levels and stream base flows 

are compared to measured values. The GFLOW model was 

calibrated with the aid of the parameter estimation program 

PEST (Doherty, 2012). The PEST program automatically 

adjusts parameter values within a specified reasonable range 

and compares simulated groundwater levels and stream base 

flows to measured water levels and base flows after each run 

of the GFLOW model. With traditional trial-and-error methods, model calibration is considered complete when simulated 

and observed water levels and flows match “reasonably” well 

and values for parameters (in this model, hydraulic conductivity and recharge) are considered “reasonable.” By contrast, the 

primary benefit of a properly constructed parameter-estimation 

routine is the capacity to automatically calculate parameter 

values that are a quantified best fit between simulated and 

observed data while also ensuring that parameter values are 

reasonable. 

Although a steady-state model was used (in which 

groundwater levels do not change with time), measured water 

levels used for calibration spanned many years, and the location of the data points can be somewhat uncertain. Groundwater-level targets (water levels in wells) were binned into four 

categories (table 6A) and assigned weights for the calibration 

process based on their locational accuracy and the amount and 

type of historical measurements. Assessing the quality and 

priority of calibration targets and assignment of corresponding 

weights is important because targets with higher weights are 

honored by the calibration process more closely than targets 

with low weights. Although the relative difference among 

target weights determines their respective influence on the 

calibration, it is useful to compute weights on the basis of 

measureable or inferred estimates of accuracy or uncertainty. 

Five wells with long-term continuous or periodic water-level 

measurements recorded in the USGS National Water Inventory System (NWIS) database (Dempster, 1990) were assigned 

the highest weight (1.2), reflecting a relatively low uncertainty 

(in terms of location, elevation, well construction, and waterlevel measurement precision) associated with these wells. This 

weight roughly equates with the 95-percent confidence level 

that the true steady-state water level at each well is within 

1.6 ft of the target value (Hill, 1998). Water levels measured 

in monitoring wells in the area that were installed for site-specific investigations were assigned the second highest weight 

(0.39) for water-level targets. These targets incorporated accurate locational and well-construction information, but the duration and frequency of water level monitoring varied. Thus, the 

true water levels for historical monitoring wells were expected 

to be within 5 ft of the target value at the 95-percent confidence level. Water levels for seepage lakes expected to be well 

connected with the aquifer also were used as targets, with the 

target water level derived from topographic maps. These seepage lake water levels are subject to uncertainty associated with 

seasonal and long-term lake fluctuations and the degree of 

connection with the underlying aquifer. Thus, water levels for 

seepage lake targets were expected to be within about 11 ft of 

the true value. Lastly, well-construction information was 

obtained from Well Construction Reports (WCRs; Wisconsin 

Department of Natural Resources, 2010), which contain limited information on the location and elevation of the well site. 

On the basis of WCR information, the well locations were 

estimated to the nearest quarter-quarter section for wells that 

lacked latitude and longitude coordinates. For instances where 

multiple WCRs were reported for the same location, target 

water levels were computed as the average water level for all 

wells. Weights assigned to WCR targets varied (0.05 to 0.18) 

depending upon reported and assumed locational and elevation data. Thus, water levels for WCR targets were expected to 

be accurate to within about 11 to 39 ft. The weighted residuals 

between the target and simulated values were used by PEST to 

determine the calibrated best fit. 

In addition to water-level targets, base-flow estimates 

were computed for 40 streamflow sites; the gain in flow 

between nine sets of streamflow sites also were computed and 

used as targets for the model calibration (table 6B). Gebert and 

others (2011) estimated base flow at 123 long-term gaging stations and 1,495 partial-record stations in Wisconsin, of which 

4 gaging stations and 23 partial-record stations were available 

for calibrating the Menominee Indian Reservation groundwater-flow model. Gebert and others (2011) estimated base flow 

at gaging stations by means of the Base-Flow Index method 
16 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

Table 6A. Calibration-target descriptions and results for groundwater and lake-level targets, and associated weights 

used for calibration with the parameter estimation program PEST. The equivalent standard-deviation-calculated 

95-percent confidence interval represents the expected accuracy of an individual target value, in the form of “plus or 

minus (the specified number of) feet.”

[std, standard deviation; CI, confidence interval]

Target type Number of 

targets

Mean 

error 

(feet)

Mean 

absolute 

error 

(feet)

Root 

mean 

square 

error 

(feet)

Weight 

(1/std)

Equivalent 

standard-deviationcalculated

95% CI

(feet)

Measured water levels 

for long-term wells in 

the USGS groundwater 

observation network

5 3.3 6.0 6.4 1.2 1.6

Historical measured water 

levels in monitoring 

wells in the model area

82 7.7 10 13.8 0.39 5

Well Construction Report 

(WCR) data with latitude 

and longitude data

650 –3.4 11.8 15.6 0.005 to 0.18 11 to 39

Seepage lake water level 

reported on topographic 

maps

34 –0.6 9.9 13.5 0.18 11

Table 6B. Calibration-target descriptions and results for stream base flow and base-flow gain targets, and associated weights 

used for calibration with the parameter estimation program PEST. The equivalent standard-deviation-calculated 90-percent 

confidence interval represents the expected accuracy of an individual target value, in the form of “plus or minus (the specified 

number) percent.”

[ft3/s, cubic feet per second]

Stream target type

Range in target 

base

flow (ft3

/s)

Number of

targets

Mean error

(ft3

/s)

Mean 

percent 

error

Equivalent coefficient of 

variation-calculated 

90% confidence interval

(used to compute target weight)

Historical gaging station 10 to 620 4 16.0 8% 5%

Partial-record gaging station 1.9 to 121 23 4.6 −64% 10%

Miscellaneous measurement site 1.3 to 145 13 1.5 6% 20%

Gain in flow between select 

streamflow targets

0.6 to 59.2 9 0.8 −220% Average of target types used to 

compute the gain in flow
Simulation of the Regional Groundwater-Flow System 17

(Wahl and Wahl, 1995; Institute of Hydrology, 1980a; 1980b) 

for daily streamflows from water year 1970 to water year 

1999. A relation-line approach, specific to each partial record 

station and paired gaging station, was used by Gebert and 

others (2011) to estimate base flow at partial-record stations. 

The estimated base flow values computed by Gebert and others (2011) for both gaging stations and partial-record stations 

were used directly as targets for the groundwater-flow model 

calibration. In addition, 13 miscellaneous streamflow measurements within the reservation were used to estimate base flow 

with the statewide regression equation developed by Gebert 

and others (2011). Weights assigned to base flow targets in 

PEST (table 6b) were based on the type of base-flow target 

(gage, partial-record, or miscellaneous measurement site) and 

the amount of flow so that the influence of both small and 

large streams on the parameter-estimation process was roughly 

similar to the influence of water-level targets; that is, the accuracy of each group of base-flow targets and their representation of long-term average conditions was estimated by using a 

90-percent confidence level (Hill, 1998) such that target values 

were as follows: (1) gaging stations were expected to be within 

5 percent of the “true” long-term base flow, (2) partial-record 

stations were expected to be within 10 percent of the “true” 

long-term base flow, and (3) miscellaneous measurement sites 

were expected to be within 20 percent of the “true” long-term 

base flow. The weight assigned to individual base-flow targets 

were then determined according to the target flow value so that 

targets of the same type had approximately equal influence on 

the calibration process.

Only a subset of all possible parameters was estimated 

by PEST. Parameters were excluded or modified if they 

were insufficiently sensitive for automated calibration (for 

example, riverbed resistance) or highly correlated with another 

parameter (for example, outwash deposits in the Antigo 

area). In these cases, the parameter was fixed at a reasonable value (riverbed resistance) or tied to another parameter 

(K_antigo_outwash_mult in table 3 was tied to K_sandy_till). 

Hydraulic-conductivity parameter values used in the calibrated 

model are listed in table 3 and shown on figure 4; effective 

hydraulic conductivities are listed in table 4 and also included 

on figure 4 after adjusting for the simulated saturated thickness, as described under the heading “Model Construction 

and Assumptions.” As described earlier, calibrated hydraulicconductivity parameter values used in the model are biased 

low compared with effective hydraulic-conductivity estimates 

(table 4), owing to the exaggerated saturated thickness used 

to improve model performance. Thus, the following discussion of aquifer properties focuses on the more transferable 

effective hydraulic-conductivity estimates that better align 

with actual saturated thicknesses. The average effective 

hydraulic conductivity for the majority of the modeled area 

representing mixed glacial outwash and till extending from 

the northwest to the southeast parts of the reservation (zones 

1, 2, 3, and 7) was estimated to be about 16 to 22 ft/d (table 

4, fig. 4). Higher average effective hydraulic conductivities 

(43 to 65 ft/d) were estimated for areas of sandy till (zone 4) 

and a sand and gravel outwash plain in the northwestern part 

of the model (zone 5). This area of relatively higher hydraulic 

conductivity corresponds with mapped geologic units (Farrand 

and others, 1984) and areas with a relatively lower density of 

perennial streams. Low stream density is associated with high 

transmissivity, in that recharge to the water table is effectively 

transmitted through the aquifer, thereby limiting groundwater 

mounding and discharge to the land surface, which in turn limits land-surface flooding and the generation of surface streams. 

Glacial deposits of loamy till were simulated with effective 

hydraulic conductivities of about 4 ft/d (zone 6). Recharge to 

the Menominee Indian Reservation was estimated at 8.4 in/yr 

(table 3). 

Following calibration, model-simulated groundwater 

level and streamflow show a close fit to target values (figs. 5 

and 6). Unweighted statistics relating all target water levels 

to all simulated levels include a mean difference of −2.0 ft 

(negative indicates that target values are, on average, less than 

simulated values), a mean absolute difference of 11.5 ft, and 

a root mean squared difference of 15.3 ft. Highly weighted 

targets generally had smaller residuals than targets given lower 

weight (table 6A). Simulated water levels generally matched 

measured (target) water levels over the entire 745-ft range of 

measured water levels (fig. 5) with little spatial bias (fig. 7). 

All simulated base flows at the gaging stations were within 

21 percent of estimated base flows, with a mean percent error 

of 8 percent. The mean percent error for partial-record stations (−64 percent) and miscellaneous sites (6 percent) further 

demonstrate a close fit to target base-flow values (table 6B and 

fig. 6). The mean percent error for gains in base flow between 

upstream and downstream target sites (−220 percent; table 6B) 

is dominated by a few targets with very low gains in base flow 

between sites. However, simulated base flow at the upstream 

and downstream targets were generally in close agreement 

with target values (fig. 6). This result illustrates the limited 

ability of the regional-scale GFLOW model to simulate smallscale gains in base flow without modifying properties, such 

as riverbed resistance, for site-specific elements or the local 

hydraulic conductivity. Given that the objectives of this modeling effort were regional in scale, site-specific parameterization was not pursued. Regardless, perfect agreement between 

measured and simulated values from groundwater-flow 

models is not expected, because of uncertainties associated 

with some target data (location, elevation, seasonal variability, 

well construction, and so on) and simplifications inherent in 

constructing models of complex natural systems.
18 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

EXPLANATION

Target water level, in feet above NAVD 88

Simulated water level, in feet above NAVD 88

700

800

900

1,000

1,100

1,200

1,300

1,400

1,500

1,600

700 800 900 1,000 1,100 1,200 1,300 1,400 1,500 1,600

Water levels from drillers’ reports (WCRs)

Lake water levels from topographic map

Historical monitoring-well data

USGS Observation Network wells

Figure 5. Simulated water levels in relation to target water levels, Menominee Indian Reservation, Wisconsin.
Simulation of the Regional Groundwater-Flow System 19

0.1

1

10

100

1,000

EXPLANATION

Target base flow or base-flow gain, in cubic feet per second

Simulated base flow or base-flow gain, in cubic feet per second

Gain in flow between base-flow targets

Miscellaneous measurement sites

Partial-record stations

Long-term gaging stations

0.1 1 10 100 1,000

Figure 6. Simulated base flow and gains in base flow in relation to target base flow and gains in base flow, Menominee 

Indian Reservation, Wisconsin.
20 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

GF

GF

GF

GF

GF

GF

GF

GF

GF

*#

*#

*# #*

*#

*#

*# #*

*#

*#

*#

*#

#*

*#

*#

#*

*#

*#

*#

*#

#*

#*

*#

*#

*#

#* #*

#*

#*

*#

#*

*#

*#

*#

*#

*#

*#

#*

*#

*#

89°10' 89° 88°50' 88°40' 88°30'

45°10'

45°

44°50'

0 4 8 MILES

0 4 8 KILOMETERS

Base map compiled from digital data sources

EXPLANATION

Simulated water-level residual

Greater than 40 feet too high

11−40 feet too high

Within 10 feet of target value

11−40 feet too low

Greater than 40 feet too low

Simulated base-flow residual (triangle)

and base-flow gain residual (cross)

*# Greater than 40 percent too high

*# 10−40 percent too high

*# Within 10 percent of target value

*# 10−40 percent too low

*# Greater than 40 percent too low

GF

GF

GF

GF

GF

Inhomogeneity boundary

County boundary

Well Construction Report data

Seepage lake

Historical monitoring well

Long-term observation well

Water-level target type Miscellaneous measurement site

Partial-record gaging station

*# Historical gaging station

Base flow target type

#*

*#

Nonrouted near-field linesink

Routed near-field linesink

Far-field linesink

Figure 7. Water-level- and base-flow-target residuals for the groundwater-flow model of the Menominee Indian 

Reservation, Wisconsin. The size of base-flow gain targets depends on the base flow target types from which the 

gains were computed.
Evaluation of Simulated Results of the Menominee Indian Reservation Model 21

Sensitivity and Uncertainty Analysis

Some uncertainty about the accuracy of models is inevitable because the model parameter values are never exactly 

known. However, the importance of each input parameter and 

its effect on simulation results can be evaluated through sensitivity tests in which the value of a parameter, such as hydraulic 

conductivity, is adjusted above or below the calibrated value 

and the magnitude of changes in simulated groundwater levels 

and base flows are quantified. In this study, PEST was used to 

calculate the sensitivity of all water-level and base-flow targets to changes in each parameter value during the calibration 

process. For the final calibrated parameter values, sensitivities 

computed by PEST (table 3) indicate that water levels and 

streamflows were generally most sensitive to the hydraulicconductivity parameters representing the glacial outwash and 

till zones covering the majority of the modeled area. Moderately sensitive parameters generally included recharge applied 

to the model and the hydraulic conductivity of the smaller 

zones representing local till deposits. Relatively insensitive 

parameters included sediment resistance for rivers and lakes. 

Initial sensitivity analyses showed similar results and were 

used to guide the selection of parameters for estimation.

Although sensitivities are useful for identifying parameters that are important for matching calibration targets, it 

is not always the case that predictions, such as the area that 

contributes recharge to wells, will be sensitive to these same 

calibration parameters; that is, the area contributing recharge 

may also be sensitive to parameters, such as the sediment 

resistance for rivers and lakes, which were identified as being 

of low sensitivity to water-level and base-flow targets during 

the calibration process. To address this potential confounding 

factor, a Monte Carlo technique (Starn and Bagtzoglou, 2012) 

was employed for scenario testing in this report whereby input 

parameter values were systematically adjusted to generate a 

range of simulated outcomes. Parameter covariance matrices 

computed by PEST during the calibration process were used 

to guide the range around the calibrated values that were 

tested for all calibrated parameters. For parameters with low 

sensitivity that were not calibrated, such as lakebed and riverbed resistance, a large variance was specified so that the range 

of tested values extended more than an order of magnitude 

from the specified value used in the calibrated model; that is, 

the range of simulated parameter values was higher for insensitive parameters than for more sensitive parameters, because 

the calibration data provided more information to constrain 

the highly sensitive parameters (hydraulic conductivity and 

recharge). The range of model parameter values were sampled 

using a Latin Hypercube approach (Starn and Bagtzoglou, 

2012) to construct about 800 parameter realizations (combinations of parameter values) for each Monte Carlo simulation. 

For example, to estimate the area contributing recharge to a 

pumped well, each parameter in table 3 was assigned about 

800 different values that typically ranged higher and lower 

than the calibrated value by one to two times the calibrated 

value. For each realization, the GFLOW model was solved 

and particle traces were evaluated to identify the area contributing recharge to selected wells for that particular realization. (A description of particle tracking is included later under 

the heading “Simulation of Areas Contributing Recharge to 

Existing Wells.”) Finally, a map of the probability that a well 

would capture water recharged in adjacent areas was produced 

by evaluating the area contributing recharge for every realization and was summarized in terms of the percentage of model 

solutions for which an individual particle of water that started 

at the water table flowed through the aquifer and was captured 

by the well. This type of visual representation of probability 

accounts for uncertainty in the parameter values as they relate 

to the particular scenario of interest, the model design, and the 

information content contained by the calibration targets; that 

is, the Monte Carlo uncertainty analysis is predicated upon the 

regional model design and calibration and is amenable to local 

refinements that could modify and improve future uncertainty 

analyses for local areas. A similar uncertainty analysis using 

Monte Carlo techniques was performed to estimate the probability of wastewater presence in the aquifer downgradient 

of wastewater-treatment lagoons near the village of Neopit. 

Application and results of the Monte Carlo simulations are further described in the subsequent sections “Simulation of Areas 

Contributing Recharge to Existing Wells” and “Simulation of 

Flow of Wastewater from Infiltration Lagoons.”

Evaluation of Simulated Results of the 

Menominee Indian Reservation Model

Results from the regional GFLOW model illustrate simulated water levels and groundwater-flow directions. Figure 8 

shows contours of the simulated water table and the path of 

several simulated particles of water flowing downgradient, 

starting at the water table and flowing to discharge locations 

such as wells and streams. The regional groundwater-flow pattern is generally parallel to the direction of flow in most rivers 

in the reservation; that is, groundwater west of the Wolf River 

generally flows from the northwest toward the southeast, as 

do most of the tributary rivers. A groundwater divide forms 

between the Wolf and Oconto Rivers (fig. 1) in the eastern part 

of the reservation. East of the Wolf River, groundwater flows 

both west toward the Wolf River and east toward the Oconto 

River. Imprinted on top of the regional flow system are local 

flow patterns that are affected by local aquifer properties and 

river geometries. As evident by the simulated particle traces in 

figure 8, much of the water that enters the aquifer as recharge 

at the water table within the reservation discharges to nearby 

rivers. 

The GFLOW model for the Menominee Indian Reservation integrates information over a large area and was designed 

to address regional-scale issues. In addition, the model forms 

a basis for local refinement or extraction to a three-dimensional model, as necessary. For example, zones of differing 

hydraulic conductivity are in the model to represent mapped 
22 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

960

980

940

920

900

1020

1000

1040

880

1080

1100

1120

1140

1160

860

1200

1220

1240

1260

1280

1300

1320

1360

820

1400

840

980

940

920

940

88°50' 88°40' 88°30'

45°

44°50'

Nonrouted near-field linesink

Routed near-field linesink

Far-field linesink

Inhomogeneity boundary

County boundary

EXPLANATION

Simulated groundwater flow path

Simulated water-table contours— 

Contour interval 20 feet. Datum is 

NAVD 88

0 3 6 MILES

0 3 6 KILOMETERS

Base map compiled from digital data sources

Figure 8. Simulated groundwater-flow directions and water-table contours, Menominee Indian Reservation, Wisconsin.
Simulation of Areas Contributing Recharge to Existing Wells 23

geologic changes, such as relatively lower permeability loamy 

till deposits in the southeastern part of the reservation. This 

heterogeneity, together with refined linesinks near groundwater withdrawal wells, facilitates application of the model to 

evaluate the sources of water to existing wells. Possible future 

studies of specific parts of the reservation or evaluation of 

sources of water to possible new wells would benefit from an 

evaluation of the appropriate level of local detail in the model.

Simulation of Areas Contributing 

Recharge to Existing Wells

Areas contributing recharge to pumped wells on the reservation were delineated by using forward and backward particle tracking (Haitjema, 1995). Particle tracking is performed 

by using the simulated water levels and groundwater-flow 

patterns within the aquifer that are produced as part of the 

GFLOW model solution to compute the path (and rate) that 

a mathematical particle of water follows through the aquifer. 

Because this is a mathematical method, water can be tracked 

either forward (with the direction of groundwater flow) or 

backward (opposite the direction of groundwater flow) from 

multiple specified locations. Both forward and backward 

tracking were applied for this study. 

Forward tracking was used as part of the Monte Carlo 

uncertainty analysis by first generating a grid of mathematical 

particles at the water table upgradient of the wells of interest. 

Next, approximately 800 different combinations of parameter 

values spanning a reasonable range around the calibrated values (table 3) were generated, with the groundwater-flow model 

re-solved for each realization. For every realization, all particles were forward-tracked and evaluated as to whether each 

mathematical particle ultimately was or was not captured by a 

pumped well, regardless of the simulated traveltime for each 

particle. Upon completion of all realizations, the percentage 

of realizations for which each particle was captured by a well 

was computed. The starting location for each particle was then 

plotted on a map, with the color code indicating the percentage of realizations, or probability, that the particle would be 

captured by a nearby well (figs. 9–13). Only particles that had 

a probability of 1 percent or greater are shown in figures 9–13; 

additional particles beyond the colored areas were simulated 

but are not shown because the particles were never captured 

by a well during any of the approximately 800 Monte Carlo 

simulations.

Backward particle tracking was used to estimate the agebased areas contributing recharge to each pumped well (for 

example, the 5-year or 10-year area contributing recharge) for 

only the calibrated set of parameters. Age-based groundwater 

traveltimes surrounding pumped wells are often used as a 

management tool for regulating approved land uses upgradient of wells, where contamination of the aquifer could prove 

particularly problematic for the desired use of pumped water. 

Backward particle tracking involved starting approximately 

20 mathematical particles along the circumference of the 

pumped well, with each particle starting at the base of the 

aquifer in order to ensure the mapped age-based contributing 

areas were not truncated. The rate at which particles travel 

through the aquifer was computed by GFLOW as a function 

of the local groundwater-flow patterns, saturated thickness, 

recharge rate, and porosity of the aquifer. A uniform porosity of 20 percent was assumed for the entire glacial aquifer. 

Although this value is on the lower end of the typical range 

of total porosity for sand-dominated aquifers (Freeze and 

Cherry, 1979), it is a reasonable estimate of effective porosity 

for sand-dominated aquifers. Effective porosity is the proportion of the aquifer through which fluids can flow and is often 

lower than the total porosity, which is simply the ratio of void 

space to total aquifer volume, including the solid particles. As 

described earlier, the saturated thickness of the aquifer was 

exaggerated to improve performance of the GFLOW solution. Therefore, local zones with adjusted porosity values were 

added to the model in the expected upgradient areas contributing recharge to pumped wells. Adjusted porosity values 

used in the model were computed by multiplying 20 percent 

porosity by the ratio of the expected actual saturated thickness 

(based on local borings and geologic maps) and the simulated 

saturated thickness. This simple adjustment ensured that simulated groundwater-flow rates approximated rates that would 

have been simulated had more representative bedrock elevations been used for the GFLOW model. The adjusted porosity 

values and associated aquifer values are listed in table 5. Areas 

contributing recharge were delineated at 5-, 10- and 20-year 

intervals. 

The results of simulations show a range in the pattern 

of simulated areas contributing recharge for individual wells 

given the parameters evaluated through the Monte Carlo 

analysis. For example, the area contributing recharge to some 

wells, such as the village of Keshena wells (fig. 12), display a 

relatively sharp spatial contrast from areas of high probability 

of capture to low probability of capture. This sharp gradation 

in probability is likely a result of the simulated steady-state 

conditions and the relatively steep hydraulic gradient along 

the area contributing recharge. This steep hydraulic gradient 

is the result of the aquifer properties (relatively low hydraulic 

conductivity and relatively high recharge) and the relatively 

steep slope of the surface-water bodies surrounding the 

contributing area; that is, the headwater streams are several 

tens of feet higher than the downstream Wolf River near 

Keshena. These surface-water features have strong control 

over groundwater-flow directions, and reasonable changes to 

simulated aquifer properties have only moderate capacity to 

modify the flow patterns between these surface-water features. 

Conversely, the areas contributing recharge for other wells, 

such as the village of Onekewat well (fig. 12), display a diffuse gradation from areas of high probability of capture to low 

probability of capture. This diffuse gradation in probability is 

likely a result of the relatively low hydraulic gradient near the 

well due to the proximity of the well to Legend Lake and other 

lakes with similar water levels. The sensitivity of hydraulic 
24 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

>>

1220

1250

1230

1240

1210

1200

1260

1190

1270

1180

1280

1290

1170

1300

1160

1310

1150

1320

1140

1130

1330

1120

88°58' 88°56' 88°54' 88°52'

45°4'

45°2'

Adjusted porosity = 0.02

Adjusted porosity = 0.02

Probability of capture

100 percent

1 percent

EXPLANATION

5-year area contributing recharge

10-year area contributing recharge

20-year area contributing recharge

Adjusted porosity inhomogeneity

Simulated water table using calibrated 

 parameters, in 10-foot intervals above NAVD 88

Change in aquifer permeability or base 

 simulated with an inhomogeneity

Streams and lakes simulated with routed linesinks

> Pumping well

0 4,000 8,000 FEET

0 1,000 2,000 METERS

Base map compiled from digital data sources with 

aerial photography from U.S. Dept. of Agriculture (2008)

“Model Construction and Assumptions” section.

Figure 9. Simulated water-table contours and the area contributing recharge to production wells for the village of Zoar, 

Menominee Indian Reservation, Wisconsin. Calculations of the adjusted porosity values are described in Table 5 and the 
Simulation of Areas Contributing Recharge to Existing Wells 25

Figure 10. Simulated water-table contours and the area contributing recharge to production wells for the village

of Neopit, Menominee Indian Reservation, Wisconsin. Calculations of the adjusted porosity values are described 

in Table 5 and the “Model Construction and Assumptions” section.

"

""" """ """ """ """ """ """ """" """" """" """"" """"" """"" """""" """""" """""" """""" """""" """""" """""" """""" """"" """""" """""" """""" """"" """"" """""" """"" """" """" """" """" """" """" """" """" """" """ """ """""" """ """""

""

"""""""""""""""" """ """""

""

""" ""

""" """ """ """ """ """ """ """" " " "" " """"""

>>>

1120

1110

1130

1140

1150

1100

1160

1170

1090

1180

1080

1070

1190

1060

1200

1050

1210

1040

1220

88°52' 88°50' 88°48'

45°2'

45°

0 3,000 6,000 FEET

0 750 1,500 METERS

EXPLANATION

5-year area contributing recharge

10-year area contributing recharge

20-year area contributing recharge

Adjusted porosity inhomogeneity

Simulated water table using calibrated 

 parameters, in 10-foot intervals above NAVD 88

Change in aquifer permeability or base 

 simulated with an inhomogeneity

Streams and lakes simulated with routed linesinks

> Pumping well

Base map compiled from digital data sources with 

aerial photography from U.S. Dept. of Agriculture (2008)

Adjusted porosity = 0.017

Adjusted porosity = 0.021

Adjusted porosity = 0.032

Probability of capture

100 percent

1 percent
26 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

960

950

940

970

980

930

990

88°46' 88°45' 88°44' 88°43'

4°56'

4°55'

4°54'

0 1,740 3,400 FEET

0 500 1,000 METERS

Adjusted porosity = 0.028

Probability of capture

100 percent

1 percent

EXPLANATION

5-year area contributing recharge

10-year area contributing recharge

20-year area contributing recharge

Simulated water table using calibrated 

 parameters, in 10-foot intervals above NAVD 88

Change in aquifer permeability or base 

 simulated with an inhomogeneity

Streams and lakes simulated with routed linesinks

> Pumping well

Base map compiled from digital data sources with 

aerial photography from U.S. Dept. of Agriculture (2008).

4

4

4

Figure 11. Simulated water-table contours and the area contributing recharge to the production well for 

Middle Village, Menominee Indian Reservation, Wisconsin. Calculations of the adjusted porosity values are 

described in Table 5 and the “Model Construction and Assumptions” section.
Simulation of Areas Contributing Recharge to Existing Wells 27

850

860

870

840

880

890

830

900

910

930

920

820

840

88°38' 88°36' 88°34'

44°54'

44°52'

0 3,000 6,000 FEET

0 750 1,500 METERS

Adjusted porosity = 0.042

Adjusted porosity = 0.031

Adjusted porosity = 0.054

Village of 

Onekewat

well Village of 

Keshena wells

Probability of capture

100 percent

1 percent

EXPLANATION

5-year area contributing recharge

10-year area contributing recharge

20-year area contributing recharge

Adjusted porosity inhomogeneity

Simulated water table using calibrated 

 parameters, in 10-foot intervals above NAVD 88

Change in aquifer permeability or base 

 simulated with an inhomogeneity

Streams and lakes simulated with routed linesinks

> Pumping well

Base map compiled from digital data sources with 

aerial photography from U.S. Dept. of Agriculture (2008)

Figure 12. Simulated water-table contours and the area contributing recharge to production wells for the villages 

of Keshena and Onekewat, Menominee Indian Reservation, Wisconsin. Calculations of the adjusted porosity values 

are described in Table 5 and the “Model Construction and Assumptions” section.
28 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

850

860

840

870

830

88°33' 88°32' 88°31' 88°30'

44°56'

44°55'

44°54'

0 2,000 4,000 FEET

0 400 800 METERS

Adjusted porosity = 0.065

Probability of capture

100 percent

1 percent

EXPLANATION

5-year area contributing recharge

10-year area contributing recharge

20-year area contributing recharge

Adjusted porosity inhomogeneity

Simulated water table using calibrated 

 parameters, in 10-foot intervals above NAVD 88

Streams and lakes simulated with routed linesinks

> Pumping well

Base map compiled from digital data sources with 

aerial photography from U.S. Dept. of Agriculture (2008).

Figure 13. Simulated water-table contours and the area contributing recharge to production wells for the village

of Redwing, Menominee Indian Reservation, Wisconsin. Calculations of the adjusted porosity values are described 

in Table 5 and the “Model Construction and Assumptions” section.
Simulation of Flow of Wastewater From Infiltration Lagoons 29

gradients to small changes in parameter values is most apparent for the Onekewat well (fig. 12) and Redwing wells (fig. 

13). For these wells, not only is there a diffuse gradation from 

high to low probability of capture, but the area with low capture probability is relatively extensive and diffuse.

For the Zoar and Neopit wells, the highest probability 

of capturing recharge from the water table occurs immediately upgradient of the wells, followed by a narrow band of 

medium probability extending thousands of feet upgradient 

(figs. 9 and 10). The distant, low-probability portions of the 

areas contributing recharge for the Zoar and Neopit wells also 

intersect a few headwater river linesinks that may contribute 

a minor amount of water to the wells. It is important to note 

that the small tributary streams immediately upgradient of the 

Zoar and Neopit wells were simulated as being dry for the 

calibrated parameter set and appear to have little influence on 

the map of probability of capture for the wells. Field reconnaissance to identify whether these headwater streams are 

perennial may be useful and, if combined with additional local 

hydrogeologic investigation, could be used to refine representation of the local groundwater-flow system upgradient of the 

Zoar and Neopit wells.

The probable area contributing recharge to the Middle 

Village well is concentrated in a small area to the south and 

east of the well (fig. 11). The relatively small size of the area 

contributing recharge and the sharp gradient from high to 

low probability of capture reflects the relatively small distance between the well and a local groundwater mound. The 

presence of a simulated groundwater mound is supported by 

target water levels used to calibrate the model. Nonetheless, 

additional local water-level measurements would be needed to 

confirm the magnitude and specific location of the mound and, 

in turn, further reduce uncertainty associated with the area 

contributing recharge to the Middle Village well.

The probable area contributing recharge to three wells 

northwest of the village of Keshena is concentrated in a small 

area north (upgradient) of the wells (fig. 12). Similar to the 

area contributing recharge to the Middle Village well, this 

contributing area reflects the relatively small distance between 

the wells and a local groundwater mound, all of which (the 

three wells and the mound) are enclosed within a single 

hydraulic-conductivity zone. Little in the way of calibration 

data is available to evaluate the extent of the mound, although 

the natural water bodies and lithologic texture in the area support the simulation of a zone of lower hydraulic conductivity, 

the simulated flow pattern, and the associated area contributing recharge. Refinement of the local geology through field 

investigations, in conjunction with refinement of the model, 

would further reduce uncertainty associated with local areas 

contributing recharge to wells.

The highest probability of capture for the area contributing recharge to the village of Onekewat well is in the area 

immediately surrounding the well and extending hundreds 

of feet to the southwest (fig. 12). An area of relatively lower 

probability extends to the northwest and southeast of the well. 

This more diffuse low-probability area is associated with the 

low pumping rate, which induces a relatively small gradient 

toward the well that is insufficient to substantially modify the 

natural hydraulic gradient in the vicinity of the well. Moreover, 

the difference in water levels among surface-water bodies 

in the Onekewat area is small, leading to a naturally low 

water-table gradient that may be sensitive to changes in the 

model parameters. For example, the lakebed resistance, which 

controls the hydraulic connection between the aquifer and 

surrounding lakes, was varied over several orders of magnitude during the Monte Carlo simulation. Order-of-magnitude 

changes in lakebed resistance can effectively eliminate the 

lake’s influence on the local water table and induce potentially 

large changes in the local groundwater-flow pattern among 

individual Monte Carlo realizations. Simulated changes in 

the local groundwater-flow pattern among realizations are 

expected to be more pronounced in areas with low water-table 

gradients compared with high-gradient areas. High uncertainty 

associated with the riverbed and lakebed resistance parameters 

is related to the limited information provided by water-level 

and streamflow targets; that is, investigation of the groundwater/lake-water interaction near Onekewat would be useful for 

reducing uncertainty associated with local lakebed resistance 

and, in turn, the area contributing recharge to the well.

The highest probability of capture for the areas contributing recharge to the village of Redwing wells is in the area 

immediately surrounding the wells and extending hundreds of 

feet upgradient of the wells toward Long Lake (fig. 13). Areas 

of relatively lower probability extend to the north and west of 

the wells. The more diffuse low-probability areas are associated with the naturally low hydraulic gradient in the immediate vicinity of the wells and high uncertainty associated with 

the lakebed resistance. For example, the low-probability area 

between Legend and Long Lakes is likely the result of a few 

realizations for which the lakebed resistance reached values 

several orders of magnitude greater than the specified value. 

Direct investigations of groundwater/lake-water interaction 

near Redwing would be useful for reducing uncertainty associated with local lakebed resistance and the areas contributing 

recharge to the wells.

Simulation of Flow of Wastewater 

From Infiltration Lagoons

The Menominee Tribe and U.S. Indian Health Service 

recently modified the wastewater-treatment facility for the 

village of Neopit and are interested in understanding the 

probable extent of infiltrated wastewater in the downgradient 

groundwater-flow system. The current operation plans call for 

rotating wastewater discharge among each of three infiltration 

lagoons, or seepage cells, on an approximately weekly basis 

(note that the lagoons were redesigned subsequent to the aerial 

photograph in fig. 14). Considering the nearly continuous rotation among the infiltration lagoons, the total daily effluent of 

86,500 gal/d (Kevin Staus, Indian Health Service, oral commun., 
30 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

2013) was distributed among the entire area of the three infiltration lagoons on a steady-state basis in the model. Infiltrated 

wastewater from the lagoons was simulated using a specified 

recharge inhomogeneity in the groundwater-flow model.

The horizontal extent of the wastewater plume downgradient from the infiltration lagoons was simulated by using 

forward particle tracking, with particles tracked from the water 

table to a terminal surface-water body. The method considered 

only advective flow from the sources; no attempt was made to 

evaluate contaminant transport and the various processes that 

affect it. Instead, results were posed in terms of the probability 

that one or more particles would flow through a synthetic grid 

cell in the area downgradient of the lagoons. Monte Carlo 

techniques were used to estimate the probable horizontal 

extent of the wastewater plume given a range of parameter 

values that were informed by covariance matrices from the 

PEST calibration, in the same manner as was performed to 

compute probabilistic areas contributing recharge to pumped 

wells as described previously. The range of model parameter 

values were sampled by using a Latin Hypercube approach 

(Starn and Bagtzoglou, 2012) to generate about 800 parameter 

realizations for the Monte Carlo simulations. 

The horizontal extents of the probable wastewater plumes 

were evaluated according to whether one or more particles 

passed through individual cells that formed a synthetic grid 

covering the aquifer in the area of interest. A synthetic grid is 

required for this plume extent problem because the analysis 

aims to answer the question, “Does one mathematical particle 

or more pass through a specific area (as defined by individual 

grid cells)?” Hence, the entire path of each particle is of 

interest for this analysis rather than simply an evaluation of 

whether each particle ultimately arrives at a specific location, 

such as a pumped well for the areas contributing recharge 

examples described above. The synthetic grid was constructed 

with square cells, each 16 ft by 16 ft, as guided by the optimized particle-tracking step size estimated by GFLOW for the 

evaluation window (Haitjema, 1995). One particle was started 

in each synthetic grid cell (156 particles with 16-ft separation 

between each particle) over the area occupied by the infiltrating wastewater lagoons and forward tracked through the 

aquifer until it discharged to a linesink or well. 

The method is sensitive to the level of detail used to represent the synthetic grid and the number of particles released 

from the wastewater source (Juckem and Fienen, 2013). Thus, 

results for source tracking of the wastewater plume are considered to be qualitative only. Nonetheless, it is expected that the 

qualitative results will be of use for evaluating management 

alternatives. For example, areas with extreme probability of 

capture (values near 100 percent or 1 percent) are expected 

to be less sensitive to the number of particles and level of discretization than areas with more uncertain probabilities (near 

50 percent computed probability). Moreover, this qualitative 

approach could be used to identify wells or surface-water bodies in areas of concern under simulated hydrologic conditions. 

For these areas of concern, the problem could be reformulated 

into an “area contributing recharge” analysis for a well or 

surface-water body of interest, which is not sensitive to userspecified settings (discretization of synthetic grid or number of 

tracked water particles). 

Wastewater infiltrated from the current wastewater infiltration lagoons flows nearly linearly south from the lagoons 

toward Tourtillotte Creek and the West Branch Wolf River 

(fig. 14). Most of the infiltrated wastewater eventually discharges to Tourtillotte Creek. However, some of the infiltrated 

wastewater is simulated as having a low probability of flowing 

beneath Tourtillotte Creek and discharging to the West Branch 

Wolf River. Although results for the probable extent of the 

wastewater plume are considered to be qualitative, the nearly 

linear flow pattern is expected to be relatively less sensitive to 

user-specified settings than the radial flow patterns simulated 

by Juckem and others (2014).

Assumptions and Limitations

As is the case with all groundwater-flow models, the 

GFLOW model of the Menominee Indian Reservation is a 

simplification of the physical system and has corresponding limitations in model precision and in how the model can 

be used. For example, local complexities, such as wetlands, 

springs, and some irregular bedrock topographic features, 

were not explicitly simulated in the regional model. Although 

such local features were included if data were available near 

local areas of interest, such as near some pumped wells, use 

of the model to answer local-scale questions in other locations 

would benefit from local refinement, including refinement 

and (or) addition of linesinks and possibly recalibration of the 

model with a focus on newly collected local data.

Steady-state conditions were assumed to be appropriate 

for this system because hydraulic conductivity is high and 

distances between surface-water features are relatively small; 

these characteristics help dampen the effects of periodic transient stresses applied to the system (Haitjema, 1995). Steadystate assumptions, which ignore groundwater release from 

storage, are expected to simulate a greater system response to 

hydrologic stresses (such as drought or pumping) than transient simulations. In addition, the base of the aquifer system, 

which corresponds to the crystalline bedrock surface, was 

specified below mapped elevations in order to improve model 

performance and solution stability. As a consequence, adjustment of hydraulic-conductivity parameter values was necessary to ensure proper simulation of the bulk transmissivity of 

the aquifer, groundwater-flow patterns, and groundwater/surface-water interactions. Similarly, simulated porosity values in 

the model were adjusted upgradient of pumped wells to ensure 

that simulated groundwater velocities and traveltimes were not 

biased because of the higher-than-actual saturated thickness 

employed in the model.
Simulation of Flow of Wastewater From Infiltration Lagoons 31

"

""""

""""

" " "" " " " " " "" ""

"" """" " """"" " """"" "" """""" " """"""" """""""""" " """""""" """""""""" """"""""""" """"""""""" """""""""""" """"""""""""" """"""""""""" """""""""""""" """""""""""""""" """""""""""""""""" """"""""""""""""""" """""""""""""""""""" """"""""""""""""""""" """""""""""""""""""""" """"""""""""""""""""""" """"""""""""""""""""""" """""""""""""""""""""""" """"""""""""""""""""""""" """""""""""""""""""""""""" """"""""""""""""""""""""""" """"""""""""""""""""""""""" """""""""""""""""""""""""""" """"""""""""""""""""""""""""" """""""""""""""""""""""""""""" """"""""""""""""""""""""""""""" """""""""""""""""""""""""""""""" """""""""""""""""""""""""""""""" """"""""""""""""""""""""""""""" """"""""""""""""""""""""""""""" """""""""""""""""""""""""""""" """""""""""""""""""""""""""""" """""""""""""""""""""""""""""" """""""""""""""""""""""""""""" """""""""""""""""""""""""""""" """"""""""""""""""""""""""""" """"""""""""""""""""""""""""" """"""""""""""""""""""""""""" "" "" "" "" "" "" " "" """""""""""""""""""""""""""" " " " " " " " " """""""""""""""""""""""""""" """""""""""""""""""""""""""" " " " " "

"

""""""""""""""""""""""""" """""""""" """""""""" """""""" """"""" """""" """" """" "" ""

1060

1070

1050

1040

88°50' 88°49'30"

44°59'30"

44°59'

Base map compiled from digital data sources with aerial 

photography from U.S. Dept. of Agriculture (2008).

0 500 1,000 FEET

0 150 300 METERS

Probability of plume extent

High

Low

EXPLANATION

Infiltration lagoon

Simulated water table using calibrated

 parameters, in 10-foot intervals above NAVD 88

Streams and lakes simulated with routed linesinks

Figure 14. Simulated probability of the horizontal plume extent from the Neopit wastewater-treatment lagoons 

under current (2013) conditions, Menominee Indian Reservation, Wisconsin. Note that results of this method are 

considered to be qualitative only. Note also that the wastewater infiltration lagoons were redesigned subsequent 

to when this aerial photograph was taken.
32 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

Given the relatively high hydraulic conductivity of the 

aquifer, relatively high groundwater-recharge rate, and presence of perennial streams, the groundwater and surface-water 

systems are assumed to be in close hydrologic connection in 

the study area. As a result, elevations of surface-water features 

are assumed to approximate heads in the underlying groundwater system. An areal two-dimensional groundwater-flow 

model was assumed to be appropriate for this application 

because the aquifer is relatively thin (hundreds of feet thick 

compared to a study area that extends many miles), areally 

extensive, and assumed to have uniform hydraulic conductivity with depth in each hydraulic-conductivity zone. Limitations of the model arise from these assumptions; specifically, 

local three-dimensional flow and transient system response 

expected near wells, infiltration lagoons, and shorelines of 

surface-water features are not represented. In addition, heterogeneities of the groundwater system (for example, variations in hydraulic conductivity and aquifer base elevation) 

are only approximated by using zones. These limitations add 

to the uncertainty associated with estimating areas contributing recharge to wells and mapping wastewater plume extents 

(Franke and others, 1999). However, the Monte Carlo methods 

that were employed for this study take uncertainty associated with parameter values partially into account; that is, 

even though computer modeling always requires simplification of the aquifer system, the use of statistical Monte Carlo 

methods (or other uncertainty methods) provides a means by 

which the sensitivity of simulated results to the model input 

can be tested and visualized (for example, probabilistic areas 

contributing recharge to pumped wells). Other uncertainties 

resulting from the model structure were not evaluated (for 

example, alternative zonations of parameters or alternative 

conceptual models); therefore, some additional uncertainty 

in simulated results is expected beyond that presented in this 

report. Whereas the probabilistic areas contributing recharge 

are reasonable and supported by observed water-level measurements, focused data collection in combination with local 

refinement of the model parameters and boundaries could 

further decrease uncertainty in the simulated areas contributing recharge to wells. In particular, local field investigations 

of riverbed and lakebed resistance could be useful for reducing uncertainty associated with groundwater/surface-water 

interaction and could improve simulated areas contributing 

recharge to wells where surface-water bodies exert strong 

control over the groundwater-flow system. 

Results of Monte Carlo simulations for tracking of 

wastewater from the infiltration lagoons are sensitive to the 

number of mathematical water particles (156, each 16 ft apart) 

used to represent infiltrating wastewater and the level of detail 

in the synthetic grid (16- by 16-ft-square cells). The resolution of the synthetic grid matched the particle-tracking step 

size for the simulated area, as computed by GFLOW, and the 

number of particles represented one particle per synthetic grid 

cell covering the infiltration lagoons. Although these settings are thought to be reasonable (Juckem and Fienen, 2013), 

results from the source-water-tracking simulation (fig. 14) 

are considered qualitative and are to be used with caution 

for informing management decisions. Moreover, this sourcetracking approach is expected to be most useful for identifying 

areas of concern downgradient from infiltrated wastewater. 

A quantitative evaluation of the probability that a well or 

stream captures wastewater could then be performed for the 

newly identified areas of concern by use of the “area contributing recharge” method.

Summary and Conclusions

A regional, one-layer, analytic element groundwaterflow model was developed to simulate the groundwater-flow 

system in the Menominee Indian Reservation, Wisconsin. 

The model was developed by the U.S. Geological Survey 

(USGS), in cooperation with the Menominee Indian Tribe of 

Wisconsin, to contribute to the fundamental understanding 

of the region’s hydrogeology. The objectives of the regional 

model were to improve understanding of the groundwater-flow 

system, including groundwater/surface-water interaction, and 

to develop a tool suitable for evaluating the effects of potential 

water-management programs. In addition, the regional model 

was designed as a framework in which more detail could be 

added and from which three-dimensional inset models could 

be developed. Simulations made with the regional model 

reproduce groundwater levels and stream base flows representative of recent conditions (1970–2013) and illustrate groundwater-flow patterns with maps of (1) the simulated water 

table and groundwater-flow directions, (2) probabilistic areas 

contributing recharge to production wells, and (3) estimation 

of the extent of infiltrated wastewater from treatment lagoons. 

The two-dimensional model used for this study was 

developed with GFLOW, an analytic element groundwaterflow modeling code. The regional aquifer is thin and overlies 

sloping bedrock, which is problematic for analytic element 

simulations. To improve model performance, the aquifer 

base was set below the actual bedrock surface and compensatory hydraulic conductivity and porosity values were used 

to maintain realistic simulation results. The GFLOW model 

was calibrated by means of parameter-estimation techniques. 

During the calibration, hydraulic conductivity and groundwater recharge were adjusted by using the parameter-estimation 

code PEST to improve the match between simulated and 

measured water-level and base-flow targets. Sensitivities of 

calibration targets to parameters were computed by PEST to 

guide calibration, but such sensitivities do not always identify 

parameters that are important for prediction scenarios, such as 

the area contributing recharge to wells. To address this potential inconsistency, Monte Carlo techniques were employed 

for scenario testing. Parameter covariance matrices computed 

by PEST during the calibration process were used to guide 

adjustment of parameters around their calibrated values for 

about 800 simulations for a given scenario. Maps of the probability for areas contributing recharge to pumped wells were 
Summary and Conclusions 33

then produced by evaluating the percentage of simulations for 

which an individual modeled particle of water was captured by 

a well. A map of the probability for the extent of a wastewater plume downgradient of the Neopit wastewater infiltration 

lagoons was produced by evaluating the percentage of simulations for which one or more modeled particles of wastewater 

traversed a cell from a synthetic grid placed across the aquifer 

near the Neopit lagoons.

Major considerations that affect limitations of the model 

and major findings from the hydrogeologic investigations and 

groundwater-flow model simulations include the following:

• Zones of uniform hydraulic conductivity and flat aquifer base elevations were used to characterize major glacial sediment deposits and the regional bedrock slope 

from high elevations in the northwest to low bedrock 

elevations in the southeastern part of the model. Aquifer base elevations, which correspond to the crystalline 

bedrock surface, were purposefully specified below the 

mapped elevations for each zone in order to improve 

model performance and solution stability. Simulated 

flow patterns and groundwater/surface-water interactions were appropriately simulated by deriving 

effective hydraulic conductivities for each zone, which 

compensated for the exaggerated base elevations. 

Similarly, groundwater-flow velocities and mathematical particle traveltimes upgradient of pumped wells 

were appropriately simulated by adjusting porosity 

values near the wells. Nonetheless, porosity values 

were not adjusted throughout the model, and any 

future simulations of groundwater velocity and particle 

traveltimes beyond the areas contributing recharge to 

wells would benefit from local porosity adjustment to 

account for the exaggerated base elevation in the new 

area of interest. 

• Regional groundwater flow generally parallels 

the regional direction of streamflow. West of the 

Wolf River, rivers and groundwater generally flow 

from the northwest toward the southeast. A groundwater divide is present east of the Wolf River, such that 

groundwater flows west toward the Wolf River and 

east toward the Oconto River.

• Areas contributing recharge to pumped wells on the 

Menominee Indian Reservation were delineated by 

tracking modeled particles of water from the water 

table to wells in combination with Monte Carlo 

techniques to produce probability maps for the area 

contributing recharge to each well nest. The simulations show a range in sensitivity of the simulated 

areas contributing recharge to the parameters evaluated through the Monte Carlo analysis; that is, areas 

contributing recharge to some wells show a steep 

gradation from high to low probability of capture by 

the pumped well. Such steep gradients imply relative 

insensitivity to the parameters (hydraulic conductivity, 

recharge, and lakebed and river bed resistance) that 

were evaluated through Monte Carlo analysis. Areas 

contributing recharge to wells that show a diffuse gradient from high to low probability are relatively more 

sensitive to the model parameters. Factors that influence differences in the shape of contributing areas and 

the area over which probabilities transition from high 

to low include (1) the ambient water-table gradient, as 

controlled by the local geology and elevations of local 

surface-water bodies, (2) the magnitude of pumping 

and associated gradients toward the well, (3) the distance between the well and a local groundwater divide 

or mound, and (4) the degree of hydraulic connection 

between the aquifer and nearby surface water. 

• Although the application of Monte Carlo methods 

provided a means to evaluate a range of reasonable 

parameter values and map quantitative probabilities, 

the application cannot evaluate all possible sources 

of uncertainty associated with modeled results—the 

method can test only the properties included in the 

model, such as zones of hydraulic conductivity. To the 

extent that other potentially important complexities 

of the natural system were not included in the model, 

additional uncertainty is inevitable. Further reductions in uncertainty associated with areas contributing 

recharge to wells may be gained by means of local 

field investigations of water levels (monitoring wells), 

hydraulic conductivity (aquifer tests), and groundwater/surface-water interaction that would provide 

information on expected reasonable ranges for model 

parameter values and also identify new complexities 

for incorporation into the model. 

• The likely extent of wastewater in the aquifer downgradient of the Neopit wastewater lagoons was simulated 

by using the groundwater-flow model and Monte Carlo 

techniques. Wastewater infiltrated from the lagoons 

flows predominantly south toward Tourtillotte Creek, 

with low probability that some of the wastewater 

flows beneath Tourtillotte Creek and discharges to the 

nearby West Branch Wolf River. Because the method 

is susceptible to settings specified by the user, results 

for source tracking of wastewater are considered to be 

qualitative only. Although results are considered to be 

qualitative, the nearly linear flow pattern simulated for 

the probable plume extent is expected to be relatively 

less sensitive to user settings than simulated plumes 

that display radial flow patterns. Regardless, it is 

expected that the results may be useful for identifying 

potential areas of concern downgradient of the lagoons, 

which can then be evaluated with quantitative methods 

in the future.
34 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

References Cited

Anderson, M.P., and Woessner, W.W., 1992, Applied groundwater modeling—Simulation of flow and advective transport: San Diego, Academic Press, 381 p.

Attig, J.W., and Ham, N.R., 1999, Quaternary geology 

of northern Oconto County, Wisconsin: Wisconsin 

Geological and Natural History Survey Bulletin 97, plate 1, 

scale 1:100,000.

Bradbury, K.R., and Rothschild, E.R., 1985, A computer technique for estimating the hydraulic conductivity of aquifers 

from specific capacity data: Ground Water, v. 23, no. 2, 

p. 240–246.

Calver, A., 2001, Riverbed permeabilities—Information from 

pooled data: Ground Water, v. 39, no. 4, p. 546–553. 

Cannon, W.F.; Kress, T.H.; Sutphin, D.M; Morey, G.B.; 

Meints, Joyce; and Barber-Delach, Robert, 1997, Digital 

geologic map and mineral deposits of the Lake Superior 

Region, Minnesota, Wisconsin, Michigan: U.S. Geological 

Survey Open-File Report 97–455, scale 1:1,000,000, additional GIS data added in 1999, accessed November 3, 2008, 

at http://pubs.usgs.gov/of/1997/of97-455/.

Dempster, G.R., Jr., comp., 1990 National water information 

system user’s manual; Volume 2, Chapter 3, Automated 

data processing system: USGS Open-File Report 90–116 

[variously paged].

Doherty, J., 2012, PEST—Model-independent parameter 

estimation user manual (5th ed.) and addendum: Brisbane, 

Queensland, Australia, Watermark Numerical Computing.

Farrand, W.R., Mickelson, D.M., Cowan, W.R., and 

Goebel, J.E., 1984, Quaternary geologic map of the Lake 

Superior 4 degrees × 6 degrees quadrangle, United States and 

Canada (edited and integrated by Richmond, G.M., and Fullerton, D.S.): U.S. Geological Survey IMAP 1420(NL–16).

Feinstein, D.T., Fienen, M.N., Kennedy, J.L., Buchwald, C.A., 

and Greenwood, M.M., 2012, Development and application 

of a groundwater/surface-water flow model using MODFLOW-NWT for the Upper Fox River Basin, southeastern 

Wisconsin: U.S. Geological Survey Scientific Investigations 

Report 2012–5108, 124 p.

Franke, O.L., Reilly, T.E., Pollock, D.W., and LaBaugh, J.W., 

1999, Estimating areas contributing recharge to wells—Lessons from previous studies: U.S. Geological Survey Circular 1174, 14 p.

Freeze, R.A., and Cherry, J.A., 1979, Groundwater: Englewood Cliffs, N.J., Prentice-Hall, 604 p. 

Gebert, W.A., Walker, J.F., and Kennedy, J.L., 2011, Estimating 1970–99 average annual recharge in Wisconsin using 

streamflow data: U.S. Geological Survey Open-File Report 

2009–1210, 14 p. plus appendixes.

Hadley, D.W., and Pelham, J.H., comps., 1976, Glacial deposits of Wisconsin: Wisconsin Geological and Natural History 

Survey, Map 10, scale 1:500,000.

Haitjema, H.M., 1995, Analytic element modeling of groundwater: San Diego, Academic Press, 394 p.

Haitjema, H.M., 2005, Modeling lake-groundwater interactions in GFLOW: 13 p., accessed September 20, 2013, at 

http://www.haitjema.com/documents/ModelinglakegroundwaterinteractionsinGFLOW.pdf.

Hill, M.C, 1998, Methods and guidelines for effective model 

calibration: U.S. Geological Survey Water-Resources Investigations Report 98–4005, 90 p.

Hooyer, T.S., and Mode, W.N., 2007, Preliminary Quaternary 

geologic map of the northern Fox River Lowland, Wisconsin: Wisconsin Geological and Natural History Survey 

Open File Report 2007–05, map scale 1:100,000.

Hunt, R.J., Anderson, M.P., and Kelson, V.A., 1998, Improving a complex finite-difference ground water flow model 

through the use of an analytic element screening model: 

Ground Water, v. 36, no. 6, p. 1011–1017.

Hunt, R.J., 2006, Ground water modeling applications using 

the analytic element method: Ground Water, v. 44, no. 1, 

p. 5–15. 

Hunt, R.J., Haitjema, H.M., Krohelski, J.T., and Feinstein, 

D.T., 2003, Simulating ground water-lake interactions—

Approaches and insights: Ground Water, v. 41, no. 2, 

p. 227–237.

Institute of Hydrology, 1980a, Low flow studies—Report No. 1, 

Research report: Wallingford, Oxon, UK, 21 p.

Institute of Hydrology, 1980b, Low flow studies—Report No. 3, 

Catchment characteristic estimation manual: Wallingford, 

Oxon, UK, p. 12–19.

Juckem, P.F., and Fienen, M.N., 2013, Estimating prediction 

uncertainty with Monte Carlo techniques and analytic elements, in Proceedings of the MODFLOW and More 2013 

Conference, Golden, Co.: p. 235–239. 

Juckem, P.F., Fienen, M.N., and Hunt, R.J., 2014, Simulation 

of groundwater flow and interaction of groundwater and 

surface water on the Lac du Flambeau Reservation, Wisconsin: U.S. Geological Survey Scientific Investigations Report 

2014–5020, 34 p., http://dx.doi.org/10.3133/sir20145020.
References Cited 35

Kelson, V.A., Hunt, R.J., and Haitjema, H.M., 2002, Improving a regional model using reduced complexity and parameter estimation: Ground Water, v. 40, no. 2, p. 132–143.

Krohelski, J.T., Kammerer, P.A., and Conlon, T.D., 1994, 

Water resources of the Menominee Indian Reservation of 

Wisconsin: U.S. Geological Survey Water-Resources Investigations Report 93–4053, 54 p.

Martin, Lawrence, 1965, The physical geography of Wisconsin (3d ed.): Madison, Wis., University of Wisconsin Press, 

608 p.

Mickelson, D.M., 1986, Glacial and related deposits of 

Langlade County, Wisconsin: Wisconsin Geological and 

Natural History Survey Information Circular 52, 30 p.

Midwestern Regional Climate Center, 2013, Climatological 

data, accessed Aug., 2013 at http://mrcc.sws.uiuc.edu/climate_midwest/historical/temp/wi/471044_tsum.html.

Milfred, C.J., Olson, G.W., and Hole, F.D., 1967, Soil 

resources and forest ecology of Menominee County, Wisconsin, with chapters by F.P. Baxter and F.G. Goff, W.A. 

Creed, and Forest Stearns: University of WisconsinExtension, Wisconsin Geologic and Natural History Survey 

Bulletin 85, Soil Series No. 60, 203 p.

Mudrey, M.G., Jr., Brown, B.A., and Greenberg, J.K., 

1982, Bedrock geologic map of Wisconsin: Wisconsin 

Geological and Natural History Survey State Map 18, 

scale 1:1,000,000.

Oakes, E.L., and Hamilton, L.J., 1973, Water Resources of 

Wisconsin—Menominee-Oconto-Peshtigo River Basin: 

U.S. Geological Survey Hydrologic Investigations Atlas 

HA–470, sheet 1.

Olcott, P.G., 1968, Water Resources of Wisconsin—Fox-Wolf 

River Basin: U.S. Geological Survey Hydrologic Investigations Atlas HA–321, sheet 1.

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.

Starn, J.J., and Bagtzoglou, A.C., 2012, Programs for calibration-based Monte Carlo simulation of recharge areas: 

Ground Water, v. 50, no. 3, p. 472–476.

U.S. Department of Agriculture, 2008, National Agricultural 

Imagery Program (NAIP) digital orthorectified images 

(DOQ), Wisconsin: Accessed Aug 26, 2013, at 

http://datagateway.nrcs.usda.gov/.

Vogel, J.C., 1967, Investigation of groundwater flow with 

radiocarbon, in Proceedings, Symposium on isotopes in 

hydrology, International Atomic Energy Agency, Vienna, 

November 14–18, 1966: p. 355–369.

Wahl, K.L., and Wahl, T.L., 1995, Determining the flow of 

Comal Springs at New Braunfels, Texas, in Texas Water 

’95, San Antonio, Tex., August 16–17, 1995, Proceedings: 

American Society of Civil Engineers, p. 77–86.

Wisconsin Department of Natural Resources, 2010, Water well 

data files: CD–ROM.

Wisconsin Department of Natural Resources, 1998, 

WISCLAND land cover (WLCGW930): Digital data, 

accessed July, 2013 at http://dnrmaps.wi.gov/imf/imf.

jsp?site=webview.

Wisconsin State Climatology Office, 2013, Northeast Wisconsin average annual precipitation (1895–2013), accessed 

Aug., 2013 at http://www.aos.wisc.edu/~sco/clim-history/

division/graphics/WI-03-precip-annual.gif.

Appendix 1. Data from auger surveys near 

the Villages of Neopit, Zoar, and Keshena

Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin Appendix 1  39﻿ 39

Appendix 1. Data from auger surveys near the Villages of Neopit, Zoar, and Keshena. —Continued

[ft, feet; m, meters; -, not determined; UTM, Universal Transverse Mercator; NAD27, North American Datum of 1927; BLS, below land surface]

Site 

name

Estimated 

ground 

elevation 

(ft)

Depth to 

bedrock 

(ft)

Estimated 

bedrock

 elevation 

(ft)

Depth to 

water 

(ft)

Estimated 

water 

table 

elevation 

(ft)

Saturated 

thickness 

(ft)

 Easting 

(UTM 

NAD27, 

m) 

 Northing 

(UTM 

NAD27, 

m) 

Lithologic description

Village of Neopit

Neopit A 1100 30 1070 13.8 1086.2 16.2 354,105.20 4,981,777.59 Uniform reddish, medium to coarse sand. Abundant clean and white quartz 

grains, common mafic grains.

Neopit B 1095

5 1090 Dry– 0.0 355,046.13 4,981,593.87 Uniform reddish, medium to coarse sand. Abundant clean and white quartz 

grains, common mafic grains.

Neopit C 1095

2 1093 Dry– 0.0 355,791.98 4,981,189.08 –

Neopit D 1081 32 1049 6.3 1074.7 25.7 355,114.85 4,983,370.64 Dark brown, medium-grained quartz sand with some silt, clay, and organic 

content.

Neopit E 1110

8 1102 4.7 1105.3 3.3 353,909.37 4,981,148.54 –

Neopit F 1120 22 1098 4.9 1115.1 17.1 353,397.35 4,981,837.16 

–

Neopit G 1060 14 1046 7.1 1052.9 6.9 357,027.52 4,981,754.78 Reddish, medium sand with minor silty clay.

Neopit H 1100 14 1086 9.3 1090.7 4.7 356,811.09 4,982,903.30 

–

Neopit I 1070 18 1052 8.3 1061.7 9.7 357,607.55 4,982,789.89 

–

Neopit J 1090 14 1076 8.4 1081.6 5.6 356,239.27 4,983,043.99 

–

Neopit K 1100 19 1081 8.0 1092.0 11.0 356,410.57 4,984,102.67 Reddish-brown, medium sand with minor silt and clay.

Neopit L 1070 > 80 below 990 13.0 1057.0 > 67.0 355,136.93 4,982,709.93 Reddish (slightly grayer than typical), medium-grained sand with some silt 

and clay. Sandy clay at bottom of hole.

Neopit M 1080 23 1057 Dry

– 0.0 354,805.07 4,982,893.09 Reddish, medium-grained sand with some silt and clay.

Neopit N 1115 > 60 below 1055 25.0 1090.0 > 35.0 354,767.70 4,983,191.81 Reddish, medium-grained sand with some silt and clay. 

Many cobbles in lower 5 - 10 feet, making augering difficult. 

Drilling aborted at 60 feet.

Neopit O 1080 > 75 below 1005 14.4 1065.6 > 60.6 355,393.94 4,983,644.39 Reddish, medium-grained sand with some silt and clay.

Neopit P 1075 > 75 below 1000 20.8 1054.2 > 54.2 353,233.80 4,982,753.44 

–

NEO 05 1115 108 1007 65.2 1049.8 42.8 355,266.27 4,982,971.73 Moist, reddish-brown sandy and silty clay from 0 to 4 feet. Slightly moist, 

reddish-brown silty sand from 4 to 9 feet. Slightly moist, reddish-brown 

silty sand from 9 to14 feet. Dry, brown sand from 14 to 44 feet. Based 

on penetration rate, driller estimates dry conditions begin at 16 feet. 

Harder drilling at 40 feet. Dry, brown sand from 44 to 94 feet. Change in 

penetration rate at 47 feet interpreted as occasional pebbles. Volume of 

returns decreased below 74 feet. Moist, brown sand from 94 to108 feet. 

Difficult drilling from 94 to 108 feet interpreted as pebbles and cobbles. 

Driller estimates crystalline bedrock at 108 feet. Sediment on auger bit 

and 1 foot of lead auger is coarse sand and sandy clay or clayey sand at 

108 feet. Indication of water on pulled augers at about 60 feet; very wet 

on pulled augers below 70 feet.

NEO 06 1100 51 1049 Dry

– 0.0 355,446.95 4,982,893.76 Moist, reddish-brown clayey soil and sand from 0 to 4 feet. Slightly moist, 

reddish-brown clayey sand from 4 to 19 feet. Nearly dry, reddish-brown 

sand with slight amounts of silt and clay from 19 to 24 feet. Dry, brown 

sand with occasional pebbles from 24 to 51 feet. More difficult drilling 

in places between 44 feet and 51 feet. Auger refusal at 51 feet inter-

preted as crystalline bedrock at 51 feet. Augers dry from 0 to 51 feet. 

NEO 08 1110 74 1036 23.8 1086.2 50.2 354,929.40 4,983,498.11 Dark-brown sand and road gravel from 0 to 4 feet. Slightly moist, reddishbrown silty sand from 4 to 14 feet. Slightly moist, light-reddish-brown 

silty sand with occasional pebbles from 14 to 19 feet. Dry, light-reddishbrown silty sand with occasional pebbles from 19 to 39 feet. Slightly 

moist, very light reddish-brown silty sand with occasional pebbles from 

39 to 44 feet. Moist to wet, very light reddish-brown silty sand with oc-

casional pebbles from 44 to 69 feet. Also coarsens and lightens in color 

with depth and becomes very wet from 44 to 69 feet. Wet, silty finegrained sand from 69 to 74 feet. Driller estimates crystalline bedrock at 

74 feet. Sediment on pulled augers at 74 feet is very wet, sandy silt.

NEO 09 1113 25 1088 Dry

– 0.0 354,760.53 4,983,154.18 Moist, reddish-brown silty sand from 0 to 9 feet. Slightly moist, reddishbrown silty sand with occasional pebbles from 9 to 14 feet. Dry, reddishbrown silty sand with many rounded to angular pebbles (up to 1.5 inches 

in diameter) from 14 to19 feet. Dry, light-brown silty sand with many 

rounded to angular pebbles (up to 1.0 inch in diameter) from 19 to 25 

feet. Driller estimates crystalline bedrock at 25 feet.

NEO 10 1110 > 114 below 996 27.6 1082.4 > 86.4 355,016.14 4,983,240.33 Slightly moist, reddish-brown silty sand from 0 to 9 feet. Slightly moist, 

reddish-brown silty sand from 9 to 19 feet, with some pebbles or cobbles 

within the lower part of this interval. Dry, reddish-brown silty sand from 

19 to 29 feet. Dry, reddish-brown silty sand from 29 to 39 feet. Dry at 29 

feet, becoming moist by 39 feet. Moist, reddish-brown silty sand from 

39 to 49 feet. Moist, brown sandy silt from 49 to 54 feet. Wet by 54 feet. 

Wet, brown sandy silt from 54 to 74 feet. Wet, brown sand from 74 to 

114 feet. Crystalline bedrock was not reached with the available augers. 

Sediment on pulled augers is wet sandy silt.

NEO 11 1120 77 1043 26.6 1093.4 50.4 354,731.65 4,983,787.33 Red sand and soil from 0 to 14 feet. Red, fine- to medium-grained sand 

with some pebbles from 14 to 19 feet. Reddish, fine- to medium-grained 

sand from 19 to 24 feet. Red, fine- to medium-grained sand with rare 

small pebbles from 24 to 29 feet. Red, fine- to medium-grained sand 

from 29 to 34 feet, with increasing coarse fraction and rare small 

pebbles. Red, fine- to coarse-grained sand from 34 to 54 feet with rare 

small pebbles. Moist, red, fine- to coarse-grained sand from 54-59 feet 

with increasing coarse fraction and rare small pebbles. Wet, red fine- to 

coarse-grained sand from 59- 69 feet with rare small pebbles. Moist, red 

fine- to coarse-grained sand from 69-77 feet with rare small pebbles. 

Driller estimates crystalline bedrock at 77 feet. throughout the hole there 

were sections (few feet thick) encountered which caused bumping and 

chattering of drilling (interpreted as cobbles) and tight drilling (high 

torque assumed to be caused by clay). Sediment on bottom of pulled 

augers is wet sand with no substantial clay.

NEO 12 1110 19 1091 Dry

– 0.0 354,845.20 4,983,374.55 Red brown fine- to coarse-grained sand from 0-19 feet. Driller estimates 

crystalline bedrock at 19 feet.

1

 Locations of NEO 20 and 21 could not be determined from GPS due to tree cover. Location is estimated based on aerial photos and perceived location along dirt road approximately midway between 

NEO 10 and 09. 2

 Ground elevation of B1, B2, and B3 were surveyed by Menominee Indian Tribe of Wisconsin.

3

 Depth to water for B4 is an estimated average over period of record from May 1 to July 30, 2002.
40 Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin

Appendix 1. Data from auger surveys near the Villages of Neopit, Zoar, and Keshena.—Continued

[ft, feet; m, meters; -, not determined; UTM, Universal Transverse Mercator; NAD27, North American Datum of 1927; BLS, below land surface]

Site 

name

Estimated 

ground 

elevation 

(ft)

Depth to 

bedrock 

(ft)

Estimated 

bedrock

 elevation 

(ft)

Depth to 

water 

(ft)

Estimated 

water 

table 

elevation 

(ft)

Saturated 

thickness 

(ft)

 Easting 

(UTM 

NAD27, 

m) 

 Northing 

(UTM 

NAD27, 

m) 

Lithologic description

Village of Neopit—Continued

NEO 13 1110 > 119 below 991 22.9 1087.1 > 96.1 354,797.67 4,983,303.19 Red brown sandy soil from 0-9 feet with occasional small pebbles. Red 

brown fine- to medium-grained sand from 9-14 feet with some organics 

and occasional pebbles. Red brown fine- to medium-grained sand from 

14-29 feet with occasional larger pebbles. Moist, red brown fine- to 

medium-grained sand from 29-34 feet with occasional larger pebbles. 

Moist, red-brown, fine- to medium-grained sand from 34 to 69 feet 

with minor silt and occasional larger pebbles. Wet, red-brown, fine- to 

medium-grained sand from 69 to 89 feet with minor silt and occasional 

larger pebbles. No auger returns from 74 to 79 feet. Wet, red-brown, 

fine- to medium-grained sand from 89 to 99 feet with minor silt and 

occasional larger pebbles. Wet, red-brown, fine- to medium-grained sand 

from 99 to 109 feet with minor silt and occasional larger pebbles. Wet, 

red-brown, fine- to medium-grained sand from 109 to119 feet with some 

coarser sand fractions and slight amounts of silt. Crystalline bedrock was 

not reached with the available augers. Sediment on auger bit and pulled 

augers is wet silty sand with no substantial clay.

NEO 201 1100 38 1062 Dry – 0.0 354,903.95 4,983,192.61 –

NEO 211 1100 19 1081 Dry – 0.0 354,903.95 4,983,192.61 –

B12 1117.6 20 1097.6 14.4 1103.2 5.6 356,280.19 4,983,285.36 –

B22 1096.3 8 1088.3 4.8 1091.5 3.2 355,761.16 4,983,491.58 –

B32 1065.1 11 1054.1 6.3 1058.8 4.7 355,917.25 4,983,034.58 –

B43

(Neopit 

MW 1)

1055 45 1010 13.6 1041.4 31.4 355,861.00 4,982,813.56 –

Village of Zoar

Zoar 

Auger A

– 68 – 26.5 – 41.5 Reddish-brown sand with some silt and clay. Difficulty in bringing augers 

to surface interpreted by driller as indication of rock or gravel in sediments.

Zoar 

Auger B

– 43 – 8.5 – 34.5 Sand. Gravel in hole between 30 and 43 feet below ground level.

Zoar 

Auger C

– >76 – 25.8 – >50.2 0-4.5 feet, soil; 4.5-10 feet, uniform fine sand; 10-25 feet, fine to medium 

sand; 25-75 feet, silty sand.

Zoar 

Auger D

– 26 – 5.6 – 20.4 Sand the entire depth.

Village of Keshena

KE 2 910 72 838 40.9 869.1 31.1 370,199 4,972,430 Thin soil, uniform red-brown sand from surface to about 65 feet. Presence 

of some cobbles or gravel between 15 and 20 feet suggested by auger 

behavior. Clayey sand recovered from lead auger at 72 feet.

KE 3 860 >125 below 735 23.9 836.1 >101.1 370,518 4,971,530 Thin, dry soil, then clean, brown sand with few fines. Slightly darker brown 

sand with slightly more silt from 15 to 20 feet BLS. Silty sand with significantly more fines from 20 to 30 feet BLS. Sandy and clayey silt from 

30 to 50 feet BLS. Uniform light brown sandy clayey silt from 50 to 70 

feet BLS. Slightly coarser silty sand from 70 to 75 feet BLS. Uniform 

silty clayey saturated sand from 75 to 125 feet BLS. Augers packed with 

brown clay from 115 to 125 feet BLS.

KE 4 940 >125 below 815 70.0 870.0 >55.0 370,710 4,972,500 Thin soil, red-brown clayey and silty sand from 10 to 30 feet BLS. Transition from abundant pebbles (angular to rounded) to light-brown, clean, 

dry sand with sugary texture from 30 to 35 feet BLS. Increase in silt & 

clay from 100 to 115 feet BLS. Increase in silt & clay from 115 to 120 

feet BLS.

KE 5 940 >100 below 840 87.1 852.9 >12.9 370,707 4,972,230 Thin soil and reddish-brown sand from 0 to 5 feet BLS. Transition to sandy, 

silty red-brown clay from 5 to 10 feet BLS. Sandy, silty red-brown clay 

from 10 to 20 feet BLS. Very little returns from 20 to 40 feet BLS. 

Light-brown, sugary sand from 40 to 100 feet BLS.

KE 6 883 >125 below 758 26.8 856.2 >98.2 370,815 4,971,790 Thin soil with some pebbles and sand from 0 to 5 feet BLS. Uniform lightbrown sand with sugary texture and occasional coarse grains from 5 to 

35 feet BLS. Sand as above but darker in color from 35 to 70 feet BLS. 

Gravelly zone at about 70 feet. Fine to coarse, light-brown sand and 

some silt, with coarse fraction increasing with depth from 70 to 80 feet 

BLS. Poorly sorted, fine to very coarse sand from 80 to 124 feet BLS. 

Abundant angular bedrock grains, considerable silt fraction, mediumbrown color from 124 to 134 feet BLS.

KE 7 845 >130 below 715 19.1 825.9 >110.9 370,850 4,971,520 No soil. Uniform red-brown sand from 0 to 110 feet BLS. Uniform redbrown sand, some cobbles and gravel, with increasing mafic grains and 

silt and clay content with depth from 110 to 130 feet BLS.

KE 8 950 >125 below 825 75.0 875.0 >55.0 Red-brown clayey sand with abundant pebbles from 10 to 15 feet BLS. 

Clean sand (little clay, fewer pebbles) from 15 to 20 feet BLS. Clean, 

dry, brown sand from 20 to 30 feet BLS. Few returns from 30 to 35 feet 

BLS. Uniform clean, light-brown sand from 40 to 125 feet BLS.

KE 9 930 >100 below 830 43.4 886.6 >56.6 370,722 4,972,840 Thin, red, clayey soil from 0 to 5 feet BLS. Light-brown, clean, sugary 

sand from 5 to 50 feet BLS. Transition to darker brown silty sand with 

some clay from 50 to 60 feet BLS. Auger returns diminished and then 

essentially stopped from 70 to 80 feet BLS. Few returns, but occasional 

pebbles up to 1 inch from 80 to 85 feet BLS.

KE 12 930 >125 below 805 76.7 853.3 >48.3 370,263 4,972,000 Thin soil, uniform light-brown sand from 0 to 125 feet BLS.

1

 Locations of NEO 20 and 21 could not be determined from GPS due to tree cover. Location is estimated based on aerial photos and perceived location along dirt road approximately midway between 

NEO 10 and 09.

2

 Ground elevation of B1, B2, and B3 were surveyed by Menominee Indian Tribe of Wisconsin.

3

 Depth to water for B4 is an estimated average over period of record from May 1 to July 30, 2002.

Appendix 1﻿
For more information concerning the research 

in this report contact the

Director, Wisconsin Water Science Center

U.S. Geological Survey

8505 Research Way

Middleton, Wisconsin 53562-3586

http://wi.water.usgs.gov/
Juckem and Dunning—Simulation of the Regional Groundwater-Flow System of the Menominee Indian Reservation, Wisconsin—Scientific Investigations Report 2014–5237

ISSN 2328-0328 (online)

http://dx.doi.org/10.3133/sir20145237