U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigations Report 2014–5167

Prepared in cooperation with the Citizen Potawatomi Nation

Numerical Simulation of Groundwater Flow, Resource 

Optimization, and Potential Effects of Prolonged Drought for 

the Citizen Potawatomi Nation Tribal Jurisdictional Area, 

Central Oklahoma

OKLAHOMA

Citizen Potawatomi Nation

Tribal Jurisdictional Area boundary
Cover:

Front and back cover background, North Canadian River at Harrah, Oklahoma, 2014 (photograph by James 

R. Hanlon, U.S. Geological Survey).

Front cover, Block diagram of Citizen Potawatomi Nation numerical groundwater-flow model and the inset 

model domain for the Citizen Potawatomi Nation Tribal Jurisdictional Area groundwater study, central 

Oklahoma. Perspective is looking northwest.

Back cover, North (left) and south (right) water-supply wells along the North Canadian River near the 

streamflow-gaging station at Harrah, Oklahoma, January 2015 (photograph by William J. Andrews, 

U.S. Geological Survey).
Numerical Simulation of Groundwater 

Flow, Resource Optimization, and Potential 

Effects of Prolonged Drought for the Citizen 

Potawatomi Nation Tribal Jurisdictional 

Area, Central Oklahoma

By Derek W. Ryter, Christopher D. Kunkel, Steven M. Peterson, and 

Jonathan P. Traylor

Prepared in cooperation with the Citizen Potawatomi Nation

Scientific Investigations Report 2014–5167

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/.

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:

Ryter, D.R., Kunkel, C.D., Peterson, S.M., and Traylor, J.P., 2015, Numerical simulation of groundwater flow, resource 

optimization, and potential effects of prolonged drought for the Citizen Potawatomi Nation Tribal Jurisdictional 

Area, central Oklahoma, U.S. Geological Survey Scientific Investigations Report 2014–5167, 27 p., http://dx.doi.

org/10.3133/sir20145167.

ISSN 2328-031X (print)

ISSN 2328-0328 (online)
iii

Acknowledgments

The authors would like to thank the Citizen Potawatomi Nation for support and cooperation 

with the scope and direction of the project and for providing geographic information pertinent 

to the analyses. They would also like to thank Paul Barlow and Edward Banta for valuable 

assistance with the resource optimization. Joseph Nadolski and Robert Frye provided help in the 

preparation of the numerical groundwater-flow model and datasets. Thanks to Paul Barlow and 

Randal Hanson for reviews of the manuscript and the numerical groundwater-flow model.

v

Contents

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

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

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

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

Location and Description of Study Area ..........................................................................................2

Previous Investigations .......................................................................................................................2

Water Use .............................................................................................................................................2

Hydrogeological Framework .......................................................................................................................4

Aquifer Hydraulic Properties .............................................................................................................4

Potentiometric Surface .......................................................................................................................4

Conceptual Flow Model ...............................................................................................................................6

Aquifer Boundaries ..............................................................................................................................6

Stream Base Flow ................................................................................................................................6

Citizen Potawatomi Nation Numerical Groundwater-Flow Model .......................................................7

Model Configuration and Boundary Conditions .............................................................................7

Model Calibration Methods ...............................................................................................................8

Observations ................................................................................................................................9

Initial Trial-and-Error Calibration ..............................................................................................9

Parameters ...................................................................................................................................9

Automated Parameter Estimation ..........................................................................................11

Calibration Results .............................................................................................................................11

Model Sensitivity ................................................................................................................................13

Inset Numerical Groundwater-Flow Model ............................................................................................16

Optimization of Groundwater Withdrawals ............................................................................................16

Water Demands .................................................................................................................................18

Streamflow Depletion by Wells .......................................................................................................18

Optimization Methods .......................................................................................................................18

Optimization Design ..................................................................................................................18

Managed Wells .........................................................................................................................19

Pumping Constraints ................................................................................................................19

Optimization Results ..........................................................................................................................19

Effects of Prolonged Drought ...................................................................................................................19

Model Limitations ........................................................................................................................................23

Summary .......................................................................................................................................................25

References Cited .........................................................................................................................................26
vi

Figures

1. Map showing location of the Citizen Potawatomi Nation Tribal Jurisdictional 

Area, hydrological features, and the Central Oklahoma and alluvial aquifers, 

central Oklahoma .........................................................................................................................3

2. Map showing Citizen Potawatomi Nation numerical groundwater-flow model active 

area with the 1987 potentiometric surface, steady-state head residuals, locations 

of model features, and the inset model domain for the Citizen Potawatomi Nation 

Tribal Jurisdictional Area groundwater study, central Oklahoma .......................................5

3. Map showing recharge multiplier zones, horizontal hydraulic conductivity pilot 

points, and streambed hydraulic conductivity parameters used in the Citizen 

Potawatomi Nation numerical groundwater-flow model, Citizen Potawatomi 

Nation Tribal Jurisdictional Area, central Oklahoma ..........................................................10

4. Graph showing transient measured groundwater heads plotted with simulated 

heads for the Citizen Potawatomi numerical groundwater-flow model, Citizen 

Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma ....................................12

5. Graphs showing simulated and estimated streamflows for North Canadian River 

streamflow-gaging stations A. near Harrah, Oklahoma, and B. at Shawnee, 

Okla., Citizen Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma ............14

6. Graph showing scaled sensitivity of model observations in parameter groups, 

Citizen Potawatomi Nation numerical groundwater-flow model ......................................16

7. Map showing inset model active area, top active layer, and boundary condition 

cells, and locations of economic-development zones included in the optimization 

analysis, Citizen Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma ......17

8. Map showing groundwater Management Process managed wells with 

optimized pumping rates at the geothermal supply area, Citizen Potawatomi 

Nation Tribal Jurisdictional Area, central Oklahoma ..........................................................20

9. Map showing groundwater Management Process managed wells, optimized 

pumping rates, and resulting effect on induced infiltration by stream reach 

number at the Iron Horse Industrial Park, Citizen Potawatomi Nation Tribal 

Jurisdictional Area, central Oklahoma ..................................................................................21

10. Map showing predevelopment base flow and postdevelopment induced 

infiltration for the North Canadian River caused by optimized pumping at the 

Iron Horse Industrial Park, Citizen Potawatomi Nation Tribal Jurisdictional 

Area, central Oklahoma ............................................................................................................22

11. Graph showing observed North Canadian River streamflow at the Shawnee, 

Oklahoma, streamflow-gaging station and simulated flow with percent change 

during and after a hypothetical 10-year drought from 1990 to 2000 ..................................24

Tables

1. Statistical comparison of measured and simulated groundwater heads and 

stream base flows, numerical groundwater-flow model for the Citizen 

Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma ....................................12

2. Statistical summary for calibrated hydraulic parameters...................................................15

3. Numerical groundwater-flow model annual flow budget, Citizen Potawatomi 

Nation Tribal Jurisdictional Area, central Oklahoma ..........................................................15

4. Water budget for a simulation of decrease in recharge during a hypothetical 

10-year drought and effects on groundwater in storage and base flow to streams, 

Citizen Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma .......................23
vii

Conversion Factors

Inch/Pound to International System of Units

Multiply By To obtain

Length

inch (in.) 2.54 centimeter (cm)

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Area

acre 4,047 square meter (m2

)

acre 0.004047 square kilometer (km2

)

square mile 259.0 hectare (ha)

square mile 2.590 square kilometer (km2

) 

Volume

acre-foot 1,233 cubic meter (m3

)

Flow rate

cubic foot per second (ft3

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

cubic foot per second (ft3

/s) 2446.576 cubic meter per day (m3/d)

gallon per minute (gal/min) 5.451 cubic meter per day (m3

/d)

Hydraulic conductivity

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

Datum

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).

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

Abbreviations

BFI A computer program that uses daily streamflow records with a hydrograph recession 

method to separate base flow from runoff; base-flow index

CPN Citizen Potawatomi Nation

EDZ Economic development zones

ET Evapotranspiration

GHB General Head Boundary

GWM Groundwater Management Process for MODFLOW

GWM–VI Groundwater Management Process, Version Independent, for MODFLOW

GTW Geothermal well

IHW Iron Horse well

Kh Horizontal hydraulic conductivity

MNW Multi-node well

OWRB Oklahoma Water Resources Board

R2 Coefficient of determination

SFR2 Streamflow-Routing Package version 2

SWB Soil-Water Balance

Sy Specific yield

USGS U.S. Geological Survey
Numerical Simulation of Groundwater Flow, Resource 

Optimization, and Potential Effects of Prolonged Drought 

for the Citizen Potawatomi Nation Tribal Jurisdictional 

Area, Central Oklahoma

By Derek W. Ryter, Christopher D. Kunkel, Steven M. Peterson, and Jonathan P. Traylor

Abstract

A hydrogeological study including two numerical 

groundwater-flow models was completed for the Citizen 

Potawatomi Nation Tribal Jurisdictional Area of central 

Oklahoma. One numerical groundwater-flow model, 

the Citizen Potawatomi Nation model, encompassed the 

jurisdictional area and was based on the results of a regionalscale hydrogeological study and numerical groundwaterflow model of the Central Oklahoma aquifer, which had 

a geographic extent that included the Citizen Potawatomi 

Nation Tribal Jurisdictional Area. The Citizen Potawatomi 

Nation numerical groundwater-flow model included alluvial 

aquifers not in the original model and improved calibration 

using automated parameter-estimation techniques. The Citizen 

Potawatomi Nation numerical groundwater-flow model was 

used to analyze the groundwater-flow system and the effects 

of drought on the volume of groundwater in storage and 

streamflow in the North Canadian River. A more detailed, 

local-scale inset model was constructed from the Citizen 

Potawatomi Nation model to estimate available groundwater 

resources for two Citizen Potawatomi Nation economic 

development zones near the North Canadian River, the 

geothermal supply area and the Iron Horse Industrial Park. 

Groundwater pumping rates at potential well locations 

were optimized using the most recent version of the U.S. 

Geological Survey Groundwater-Management Process 

for MODFLOW. The objectives of optimization were to 

determine if a total pumping rate of 500 gallons per minute 

could be pumped from 5 wells at the geothermal supply area 

and to maximize discharge from 16 wells at the Iron Horse 

Industrial Park without exceeding specified head drawdown 

constraints at the pumping wells and thus prevent groundwater 

depletion.

The inset model was used to estimate North Canadian 

River streamflow depletion caused by optimized pumping 

at the Iron Horse Industrial Park because water quality was 

a concern, and the river may have degraded water quality 

compared to water in other parts of the alluvial aquifer. The 

fate of streamflow that infiltrates into groundwater because 

of pumping was not directly determined, but it was assumed 

that this water could end up in the well discharge, and was 

considered to be a maximum proportion of well discharge 

derived from the North Canadian River.

The total optimized continuous pumping rate from 

five managed wells at the geothermal supply area was 

638 gallons per minute, which exceeded the target pumping 

rate of 500 gallons per minute. The total continuous pumping 

rate from 16 wells at the Iron Horse Industrial Park was 

1,472 gallons per minute, which induced stream infiltration 

of approximately 4.1 gallons per minute (approximately 

0.3 percent of the total well discharge) from the North 

Canadian River. 

To estimate the effects of drought on water resources in 

the Citizen Potawatomi Nation Tribal Jurisdictional Area, a 

hypothetical 10-year drought during which precipitation would 

decrease by 50 percent was simulated by decreasing model 

groundwater recharge by the same proportion for the period 

1990–2000 of the transient model. The effects of the drought 

were estimated by calculating the change in the volume of 

groundwater storage and groundwater flow to streams at 

the end of the drought period, and the change in simulated 

streamflow in the North Canadian River at the streamflowgaging station at Shawnee, Okla., during and after the drought.

The hypothetical decrease in recharge during the 

simulated drought caused groundwater in storage over the 

entire model in the study area to decrease by 351,500 acrefeet (14,100 acre-feet in the North Canadian River alluvial 

aquifer and 346,400 acre-feet in the Central Oklahoma 

aquifer), or approximately 0.2 percent of the total groundwater 

in storage over the drought period. This small percentage of 

groundwater loss showed that the Central Oklahoma aquifer 

as a bedrock aquifer has relatively low rates of recharge from 

the surface relative to the approximate storage. The budget 

for base flow to the North Canadian River indicated that the 

change in groundwater flow to the North Canadian River 

decreased during the 10-year drought by 386,500 acre-feet, 

or 37 percent. In all other parts of the Citizen Potawatomi 
2 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

Nation Tribal Jurisdictional Area, base flow decreased by 

292,000 acre-feet, or 28 percent. Streamflow in the North 

Canadian River at the streamflow-gaging station at Shawnee, 

Okla., decreased during the hypothetical drought by as much 

as 28 percent, and the mean change in streamflow decreased as 

much as 16 percent. Streamflow at the Shawnee streamflowgaging station did not recover to nondrought conditions until 

about 3 years after the simulated drought ended, during the 

relatively wet year of 2007.

Introduction

Commercial development has grown in recent years 

in the Citizen Potawatomi Nation (CPN) Economic 

Development Zones (EDZs) in the CPN Tribal Jurisdictional 

Area (fig. 1). The EDZs require dependable sources of 

surface water and groundwater—some with specific waterquality requirements—as water demands increase. The study 

described in this report was a cooperative effort between the 

U.S. Geological Survey (USGS) and the CPN to evaluate the 

available groundwater resources for resource management in 

the CPN Tribal Jurisdictional Area by performing a detailed 

analysis of the integrated hydrological system and water use in 

areas of development, including numerical groundwater-flow 

models. 

Part of the numerical groundwater-flow model of the 

Central Oklahoma aquifer described in Mashburn and others 

(2013), referred to in this report as the Central Oklahoma 

aquifer model, was refined and improved to produce a 

subregional-scale numerical groundwater-flow model, referred 

to here as the CPN model. The CPN model was used to 

describe the groundwater-flow system and to estimate the 

effects of future droughts on the water resources in the CPN 

Tribal Jurisdictional Area (fig. 1). A high-resolution numerical 

groundwater-flow model—termed the inset model—was 

constructed in the CPN model active area as a stand-alone 

model for local-scale analysis. The inset model was used to 

optimize use of groundwater resources in two EDZs using 

state-of-the-art applications, including the Groundwater 

Management Process (Banta and Ahlfeld, 2013). 

Purpose and Scope

The purpose of this report is to describe the construction, 

calibration, and analyses performed with two numerical 

groundwater-flow models for the Central Oklahoma aquifer 

and alluvial aquifers in the CPN Tribal Jurisdictional Area 

in central Oklahoma. Analyses include optimization of 

groundwater-resource use at two locations and estimating the 

effects of a prolonged severe drought on the total amount of 

groundwater in storage and streamflow at selected locations. 

The scope of this study includes the geographic extent of 

the CPN Tribal Jurisdictional Area where it overlies the 

Central Oklahoma aquifer (fig. 1). The time period of the 

analysis of the groundwater-flow system includes the period 

1987–2009.

Location and Description of Study Area

The CPN Tribal Jurisdictional Area is in central 

Oklahoma, east of Oklahoma City on the southeastern 

part of the Central Oklahoma aquifer between the North 

Canadian River on the north and the Canadian River on the 

south (fig. 1). The study area was defined as the CPN Tribal 

Jurisdictional Area that overlies the Central Oklahoma aquifer, 

and the numerical groundwater-flow model, as defined by the 

CPN, included an area approximately 20 miles (mi) beyond 

the western boundary of this area and approximately 5 mi 

north of the North Canadian River to minimize external 

boundary effects. The study area is approximately 1,100 

square miles (mi2) (704,000 acres).

Large lakes in the study area include Lake Thunderbird, 

the Shawnee Reservoir 1 and Shawnee Reservoir 2, Wes 

Watkins Reservoir, and Tecumseh Lake (fig. 1). Land use in 

the area is mostly rural with cultivated agricultural practices 

restricted to river valleys, and other areas characterized as 

ranchland or forests (Mashburn and others, 2013).

Previous Investigations

Mashburn and others (2013) describes the hydrogeology, 

surface water, and a regional-scale transient numerical 

groundwater-flow model that included the entire Central 

Oklahoma aquifer. Havens (1989) describes the hydrogeology 

of the North Canadian alluvial and terrace deposits in the 

study area, and Adams and others (1997) is a collection 

of digital geographic data that includes the extents and 

generalized hydraulic properties and groundwater altitudes of 

aquifers in this area.

Water Use

Groundwater demands in the CPN Tribal Jurisdictional 

Area are predominantly from municipal and domestic wells 

(Mashburn and others, 2013). Water use is regulated by the 

Oklahoma Water Resources Board (OWRB) and reported by 

water-use permit holders as total annual pumping. For this 

report, monthly time increments of water pumping amounts 

were used. To convert annual water use from the Central 

Oklahoma aquifer model to monthly time periods, the OWRB 

estimate of the monthly percentage of the annual water 

discharge for different categories of wells pumped by county 

(Oklahoma Water Resources Board, 2012) was used. Other 

municipal demands were met through surface-water storage in 

reservoirs and were not included in this analysis.
Introduction 3

Lake

Tecumseh

Shawnee

Reservoir 2

North

Fork

Salt

Quapaw

Walnut

Canadian River

Little River

Creek

Creek

Creek

Stanley

Draper

Lake

Horseshoe

Lake

Lake

Thunderbird

Wes Watkins

Reservoir

Shawnee

Reservoir 1

North Canadian River

Shawnee

Norman

Moore

Midwest City

Del City

Oklahoma City

The Village 07241520

07241550

07241800

07230000

07230500

COUNTY

COUNTY

CLEVELAND

COUNTY

SEMINOLE

COUNTY

LINCOLN

COUNTY

POTTAWATOMIE

COUNTY

PONTOTOC

COUNTY

MCCLAIN

OKLAHOMA

97°30' 97°15' 97°00' 96°45'

35°30'

35°15'

35°00'

Base modified from U.S. Geological Survey digital data

Universal Transverse Mercator, zone 14

North American Datum of 1983

EXPLANATION

Alluvial aquifer 

 North Canadian alluvial aquifer (modified from Heran and others, 2003)

 Other alluvial deposits

Central Oklahoma aquifer (Mashburn and others, 2013)

Numerical groundwater-flow model active area boundary

Citizen Potawatomi Nation Tribal Jurisdictional Area boundary

Streamflow-gaging station and number 07230000

0 5 10 15 20 MILES

0 5 10 15 20 KILOMETERS

Area enlarged

OKLAHOMA

Figure 1. Location of the Citizen Potawatomi Nation Tribal Jurisdictional Area, hydrological features, and the Central Oklahoma and 

alluvial aquifers, central Oklahoma.
4 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

Hydrogeological Framework

The hydrogeological framework of the Central Oklahoma 

aquifer described in Mashburn and others (2013) was 

modified in the CPN model. Alluvial aquifers are composed of 

Quaternary-age alluvial-valley deposits along major streams 

(fig. 1) and were included as an additional discontinuous 

hydrogeologic unit in the CPN model to better simulate the 

flow between streams and groundwater. The lower boundary 

of the Central Oklahoma aquifer model was defined as the 

transition from overlying freshwater to saline water, and it 

was assumed that there was not substantial flow across this 

boundary (Mashburn and others, 2013). Although there is no 

physical hydrologic boundary and there was not expected to 

be substantial flow in the saline part of the aquifer, the CPN 

model includes the saline groundwater to the base of the 

Wellington Formation. The base of the Wellington Formation 

is a physical boundary where hydraulic conductivity greatly 

decreases (Mashburn and others, 2013). In this report, the 

saline zones of the Wellington Formation are referred to as the 

“saline unit.”

Alluvial deposits along the North Canadian River 

(Havens, 1989) and reaches of Little River and Salt Creek 

(Heran and others, 2003) were included as a discontinuous 

hydrogeologic unit (fig. 1). The OWRB designated the North 

Canadian River alluvium as the North Canadian River alluvial 

and terrace aquifer (Oklahoma Water Resources Board, 2012) 

and is referred to in this report as the North Canadian River 

alluvial aquifer. The alluvial deposits along the Little River 

and Salt Creek were not extensive or thick enough to be 

named as aquifers in this study. Other surficial deposits such 

as upland terrace deposits were not considered to contain 

substantial groundwater resources and were not included in the 

CPN model as a separate unit. 

The base and thickness of alluvial deposits were 

estimated from borehole logs submitted to the OWRB 

(Oklahoma Water Resources Board, 2013). Lithological logs 

that included a depth of the base of alluvial deposits were 

selected and used to construct the bedrock surface. Between 

locations where the base of the alluvial deposits was estimated, 

the base was assumed to be an east-draining erosional surface 

similar to the configuration of modern drainages. The North 

Canadian alluvial aquifer was estimated to have a mean 

thickness of 29 feet (ft) and maximum thickness of 77 ft in the 

study area. The Little River alluvial deposit mean thickness 

was approximately 24 ft with a maximum of 107 ft, and the 

Salt Creek alluvium mean thickness was approximately 27 ft 

with a maximum of 82 ft in the study area. 

Aquifer Hydraulic Properties

Initial hydraulic properties of the freshwater part of the 

Central Oklahoma aquifer were taken directly from Mashburn 

and others (2013). Horizontal hydraulic conductivity (Kh) 

ranged from 0.3 to 6.6 feet per day (ft/d) and locally is highly 

variable because of lenticular beds of fine sand surrounded 

by claystone and shale. The saline unit was assumed to be 

composed of silt and silty fine sand and was assigned an initial 

Kh value of 0.2 ft/d for this lithology from Fetter (1994). 

The values for Kh and vertical anisotropy used in the Central 

Oklahoma aquifer model were used as initial values in the 

CPN model and were adjusted during model calibration. There 

were, however, no available direct measurements of hydraulic 

properties in the Central Oklahoma aquifer, saline unit, or 

North Canadian alluvial aquifer in the study area.

Adams and others (1997) is a digital dataset of the 

hydrogeology of the North Canadian alluvial aquifer, but 

that dataset is limited to the average properties of the entire 

deposits as single values for the entire aquifer. Thus, to 

include an estimate of spatial variation for the CPN model, 

hydraulic properties of textural classes from Fetter (1994) 

were assigned to textural descriptions in lithological borehole 

logs, and the thickness-weighted mean Kh was calculated 

at each log location. The Kh values for the alluvial deposits 

from individual logs were estimated to range from 0.3 to 

120.0 ft/d with a mean of 38.0 ft/d. Mean Kh values from 

each lithological borehole log were used to interpolate Kh 

from various locations across the alluvial aquifer using the 

inverse-distance weighted function in ArcGIS version 10, 

Spatial Analyst (Esri, 2015). As with Kh values of the Central 

Oklahoma aquifer, estimated hydraulic properties of the 

alluvium were changed during numerical model calibration. 

Initial hydraulic properties were used to assign conductance 

values to aquifer boundary conditions.

The Kh of the alluvial material was estimated to be much 

greater than the Kh of the Central Oklahoma aquifer; thus, 

groundwater in the alluvial material is much more responsive 

to stresses such as seasonal recharge, and has a substantial 

hydraulic connection with streams such as the North Canadian 

River. The lesser Kh of the Central Oklahoma aquifer is 

assumed to limit the hydraulic connection between these two 

aquifers.

Potentiometric Surface

The potentiometric surface of the Central Oklahoma 

aquifer in 1986-7 from Mashburn and others (2013) was used 

for this study (fig. 2). The potentiometric surface of the saline 

unit was assumed to be the same as the Central Oklahoma 

aquifer. Detail for the potentiometric surface in the North 

Canadian River alluvial aquifer and Little River and Salt 

Creek alluvial deposits was added using water levels from the 

National Water Information System (U.S. Geological Survey, 

2012) and the altitude of surface-water bodies, such as stream 

channel altitudes derived from digital elevation models (not 

shown on fig. 2). 
Hydrogeological Framework 5

COUNTY

COUNTY

CLEVELAND

COUNTY

SEMINOLE

COUNTY

LINCOLN

COUNTY

POTTAWATOMIE

COUNTY

PONTOTOC

COUNTY

MCCLAIN

OKLAHOMA

97°30' 97°15' 97°00' 96°45'

35°30'

35°15'

35°00'

Base modified from U.S. Geological Survey digital data

Universal Transverse Mercator, zone 14

North American Datum of 1983

Area enlarged

on figure 6

Numerical groundwater-flow model

 Boundary condition, by cell type

 Drain

 General head boundary

 SFR2 stream cell

 Lake cell

 Top active layer

 Alluvial aquifer

 Central Oklahoma aquifer

Economic development zone

Citizen Potawatomi Nation Tribal 

 Jurisdictional Area boundary

Potentiometric contour—Shows altitude

 at which water level would have stood

 in tightly cased wells, 1986–87. Contour 

 interval 50 feet. Datum is NAVD 88 

 (Mashburn and others, 2013)

2

1

2

2

3

EXPLANATION

SFR2, Streamflow-routing package version 2

900

Pumping well

Observation well with layer number, 

 by steady-state head residual, in feet

 –39.0 to –24.3

 –24.2 to –10.1

 –10.0 to –0.2

 –0.1 to 12.2

 12.3 to 44.7

0 5 10 15 20 MILES

0 5 10 15 20 KILOMETERS

Lake

Thunderbird

2

2

2

2

2

1

1

2

1

1

1

2

2

2

2

3

2

2

2

2

2

2

2

2

4

3 2

3

2

2

2

2

2

2

2

2

2

2

2

2

3

2

2

2

2

2

2

3

1

3

2

2

3 3

2

2 2

2

2

2 3

2

2 3

2

2

2

2

2

2

2

2

2

2

2

2

2

2 2

2

2

2

2

2

2

2

2

2

2

2 2

2

2

2

2

2

2

2

1

2

2

3

2 2

2

2 3

2

2

3

2

2

2

2

2

2

2

3

2

2 2 3

2

3

2 2

2

2

3

2

2 2 2

2 2

2

2

2

2

3 1 2

2

2 2

2

1

2

2

1,050

1,000

1,100

1,100

1,050

1,100

1,050

950

1,000

900

950

1,100

1,100

1,050

1,000

950

1,200

1,150

1,100

1,150

1,050

950

900

1,100

1,000

950

1,000

1,050

Figure 2. Citizen Potawatomi Nation numerical groundwater-flow model active area with the 1987 potentiometric surface, steady-state 

head residuals, locations of model features, and the inset model domain for the Citizen Potawatomi Nation Tribal Jurisdictional Area 

groundwater study, central Oklahoma.
6 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

As shown in Mashburn and others (2013) and figure 2, 

groundwater in the Central Oklahoma aquifer flows to the 

east away from recharge areas in the central and western parts 

of the aquifer and discharges to streams across the aquifer. 

Streams are generally considered to be gaining streams that 

receive groundwater discharge along their lengths. The total 

flow to all lakes and streams has not been calculated by 

previous studies or estimated in this study; however, base 

flow for the North Canadian River was estimated and used to 

calibrate the numerical groundwater-flow model.

Conceptual Flow Model

A conceptual flow model is a schematic description of a 

groundwater system including water inflows and outflows that 

take place at aquifer boundaries. Stream base flow and gaining 

and losing stream reaches are also described in this section.

Aquifer Boundaries

Aquifer boundaries include locations, linear features, and 

areas. Location boundaries in this study included pumping 

wells used for agriculture or public water supply and seeps or 

springs. Linear boundaries were stream channels that, in most 

cases in the study area, are sinks where groundwater flows 

to streams, and the margins of the area where groundwater 

moves laterally into or out of the aquifer. Areal boundaries 

include the land surface where recharge from precipitation and 

plant consumption takes place, and lakebeds through which 

groundwater flows to or from the lake.

The boundary condition with the largest area in the 

model is the land surface through which water enters the 

groundwater system from precipitation and leaves the system 

through plant uptake and evapotranspiration (ET). Based on 

the potentiometric surface in Mashburn and others (2013), 

the water table of the Central Oklahoma aquifer was located 

within bedrock and well below the typical deciduous forest 

and grassland root zone depth (Thornthwaite and Mather, 

1957). Thus, plant uptake from groundwater through ET 

was not simulated where the Central Oklahoma aquifer was 

present at the land surface. On alluvial deposits and near 

surface-water bodies where the water table is very shallow, 

recharge was decreased during model calibration to offset 

groundwater discharged to ET. The largest discharges of the 

Central Oklahoma aquifer groundwater occur along streamchannel boundaries as stream base flow (Mashburn and others, 

2013). Stream channels also were determined to be the largest 

discharge boundaries for the North Canadian alluvial aquifer.

Stream Base Flow

The objective of stream base flow analysis was to 

estimate the mean daily base flow during the study period 

(1987–2009) and the monthly stream base flow for the 

period of the transient model (1988-2009) as described in the 

“Model Configuration and Boundary Conditions” section of 

this report for the North Canadian River in the study area. 

Stream base flow, defined in Barlow and Leake (2012) as 

the flow of groundwater into streams, was estimated along 

the North Canadian River in the study area to provide flow 

values for numerical groundwater-flow-model calibration. 

Stream reaches receiving base flow are referred to as gaining, 

and those with infiltration to the aquifer are referred to as 

losing. The North Canadian River reaches studied were 

between the USGS streamflow-gaging stations just upstream 

from the model domain at Britton Road in Oklahoma City, 

Okla. (USGS station 07241520), near Harrah, Okla. (USGS 

station 07241550), and at Shawnee, Okla. (USGS station 

07241800) (fig. 1). Streamflow measurements for all stations 

were accessed from the National Water Information System 

(U.S. Geological Survey, 2012). The period of record for the 

streamflow-gaging station at Harrah was from 1968 to 2009, 

but the period of record for the streamflow-gaging station 

at Shawnee was from 2001 to 2009 and did not include the 

complete period of study (U .S. Geological Survey, 2012). For 

the available periods, base flow was estimated using base-flow 

index (BFI), a computer program that uses daily streamflow 

records with a hydrograph recession method to separate base 

flow from runoff (Wahl and Wahl, 2007). 

During the 84 months for the period from 2001 to 

2009 during which all three streamflow-gaging stations 

recorded streamflow, the North Canadian River between the 

streamflow-gaging stations at Britton and Shawnee was a 

gaining stream during 70 months and a losing stream for 14 

months. From 2001 to 2009, the mean daily North Canadian 

River base flow calculated using BFI at the streamflow-gaging 

station at Britton was 99 cubic feet per second (ft3

/s), at Harrah 

was 185 ft3/s, and at Shawnee was 226 ft3/s. 

To estimate base flow for the reach between the Harrah 

and Shawnee streamflow-gaging stations when streamflow 

observations were not available at the Shawnee station 

(1988–2001) during the model simulation period, the base 

flow upstream from the streamflow-gaging station at Harrah 

was assumed to change by a characteristic amount calculated 

for the period when streamflow data were available at both 

stations. The BFI results indicate that the base flow between 

the streamflow-gaging station at Harrah and the streamflowgaging station at Shawnee was affected by the precipitation 

measured at the National Weather Service cooperative 

observer station in Shawnee, Okla. (USC00348110; National 

Climatic Data Center, 2013). Months were classified based on 

precipitation as being wet (greater than 4.65 inches), average 

(between 4.65 and 1.53 inches), or dry (less than 1.53 inches). 

During wet months, mean base flow at the streamflow-gaging 

station at Shawnee increased 77 ft3/s relative to flow at the 

streamflow-gaging station at Harrah; mean base flow increased 

36 ft3/s during months with average precipitation, and mean 

base flow increased 44 ft3/s during dry months. Streamflow 

for months without streamflow data at the Shawnee 
Citizen Potawatomi Nation Numerical Groundwater-Flow Model 7

streamflow-gaging station was added by classifying each 

month as wet, average, or dry and adding the characteristic 

change in base flow to the measured flow at Harrah.

The BFI program was used to calculate base flows for 

the Little River streamflow-gaging station near Tecumseh, 

Okla. (USGS station 07230500), although these data were not 

used as values for model calibration. Because of effects from 

streamflow releases from Lake Thunderbird, the streamflow 

input to BFI was estimated to be the difference between 

Little River streamflow near Tecumseh and streamflow just 

downstream from the lake measured at the streamflow-gaging 

station at Little River below Lake Thunderbird near Norman, 

Okla. (USGS station 07230000). The long-term mean base 

flow at that streamflow-gaging station on the Little River from 

1968 to 2009 was 2.69 ft3

/s.

Citizen Potawatomi Nation Numerical 

Groundwater-Flow Model

The CPN model was a simulation of the groundwaterflow system using MODFLOW-NWT, a Newton formulation 

for MODFLOW-2005 (Niswonger and others, 2011) that 

included both a steady-state and transient model. The 

potentiometric heads produced by the steady-state model 

were used as the starting heads for the transient model. The 

calibrated transient model was used to estimate the effects of 

drought on water resources of the study area. A transient inset 

model was constructed using calibrated parameters from the 

CPN model for groundwater-pumping optimization.

The numerical groundwater-flow models described in this 

report were based on the assumption that groundwater flows 

through the aquifers according to Darcian flow principles, 

that the variation in the hydraulic properties of the aquifers 

modeled can be resolved with the chosen cell sizes, and that 

the variations in water flux can be simulated with monthly 

stress periods. It was assumed that the water lost to plant ET 

can be simulated by adjusting recharge to the aquifer. It was 

assumed that during 1987 the groundwater-flow system was at 

an approximate equilibrium and could represent a steady-state 

period as used in Mashburn and others (2013). Lastly, the head 

and base-flow observations and estimates were assumed to be 

adequate to calibrate the models.

Model Configuration and Boundary Conditions

The CPN model objectives required more spatial detail 

than the Central Oklahoma aquifer model, which used cell 

dimensions of 3,280 ft on each side. More spatial detail 

was achieved by setting the CPN model cells to 820 ft on 

each side. One of the most important modifications to the 

Central Oklahoma aquifer model was the addition of the 

North Canadian alluvial aquifer and alluvial deposits along 

Little River and Salt Creek as layer 1 (fig. 2), which were not 

included as a separate layer in the Central Oklahoma aquifer 

model.

The CPN model did not require the same level of vertical 

detail included in the Central Oklahoma aquifer model, 

which used 100-ft thick horizontal layers and did not include 

distinct hydrogeologic units. In the CPN model, the Central 

Oklahoma aquifer was configured as two layers (layers 2 

and 3) representing the upper and lower parts of the aquifer. 

The saline unit was assigned to layer 4. Where there was no 

alluvial deposit, the uppermost active model layer was the 

Central Oklahoma aquifer in layer 2 (fig. 2). 

The CPN model time period including both steady-state 

and transient models was from 1987 to 2009, which was 

the same as the transient Central Oklahoma aquifer model 

of Mashburn and others (2013). The steady-state model 

represented mean conditions during 1987, and the transient 

model included the period from 1988 to 2009. The Central 

Oklahoma aquifer model used annual stress periods, which 

were too long to resolve the changes in head and flows in the 

alluvial aquifers used in the CPN model. Thus, the transient 

CPN model used 264 monthly stress periods to increase the 

temporal resolution.

Conceptual flow-model boundaries were simulated in the 

numerical groundwater-flow model by defining model cells 

as having boundary conditions. Boundary conditions included 

head-dependent flow boundaries such as streams, springs, 

and locations where groundwater enters or leaves the model; 

constant-rate flow boundaries such as wells or areas receiving 

recharge; and no-flow boundaries defined as inactive cells or 

the model extent. Flow-boundary conditions were simulated 

with several packages that function with the MODFLOWNWT solver.

Streamflow was routed through the model using 

the Streamflow-Routing Package version 2 (SFR2) for 

MODFLOW (Niswonger and Prudic, 2005) (fig. 2). This 

package simulates water in stream channels and flow between 

streams and groundwater with specified parameters for stream 

channel dimensions, streambed thickness and hydraulic 

conductivity, and channel gradient. For the streams overlying 

the Central Oklahoma aquifer (layer 2), the parameter values 

derived in Mashburn and others (2013) were used as initial 

values. For stream segments overlying the Salt Creek, Little 

River, and North Canadian River alluvial deposits, initial 

streambed conductivity was estimated to be within the range 

of silty, fine sand in Fetter (1994), approximately 13 ft/d. 

Streambed conductivity in four segments of the North 

Canadian River was adjusted during automated parameter 

estimation.

The SFR2 Package was linked with the Lake Package 

for MODFLOW (Merritt and Konikow, 2000), which was 

used to simulate Lake Thunderbird (fig. 2). Stream segments 

simulated in SFR2 can discharge into lakes, and reservoirs can 

release water into stream segments. Streamflow from local 

base-flow fed streams was routed into Lake Thunderbird, and 

releases from Lake Thunderbird recorded by the U.S. Army 
8 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

Corps of Engineers (2013) were routed to the adjacent 

downstream segment of the Little River. The Lake Package 

also simulated groundwater flow into and out of Lake 

Thunderbird.

Several other head-dependent flow boundary conditions 

in the steady-state and transient models were simulated 

using the General Head Boundary (GHB) Package for 

MODFLOW (Harbaugh and others, 2000). A GHB cell has 

a specified groundwater head for each stress period and a 

hydraulic conductance of the interface between the GHB 

and the adjacent aquifer cell, which governs groundwater 

flow between the GHB and adjacent cells. For GHB cells, 

the conductance was set at the initial hydraulic properties of 

the model at that cell and was adjusted only during the initial 

calibration.

Shawnee Reservoirs 1 and 2, Horseshoe Lake, Wes 

Watkins Reservoir, and Tecumseh Lake (fig. 2) did not release 

substantial water to streams and were simulated using GHB 

cells instead of the Lake Package. Approximate lake stages 

were taken from the National Hydrography Dataset (U.S. 

Geological Survey, 2013b) and held constant during the 

transient model period.

Groundwater entered the model from areas to the west 

and north where the edge of the CPN model did not coincide 

with a hydrological boundary in the Central Oklahoma 

aquifer model through GHB cells (fig. 2). In these areas, the 

heads at GHB cells were taken from the groundwater heads 

simulated by the Central Oklahoma aquifer model (Mashburn 

and others, 2013), and the extent of the CPN model was set 

beyond the CPN Tribal Jurisdictional Area to minimize any 

boundary effects from the GHBs. On the western upstream 

model margin, groundwater entering the North Canadian 

River alluvial aquifer was simulated using GHBs. Heads in 

GHBs were set at the estimated initial heads from groundwater 

observations and surface-water bodies. The GHB cells were 

held constant during the transient model period. As with the 

GHBs used to simulate lakes, conductance of the peripheral 

GHBs was approximately the same as the aquifer in which 

they were located.

Because of the increased detail of the CPN model, 

springs and seeps that discharge groundwater along the eastern 

margin of the Central Oklahoma aquifer were simulated with 

the Drain Package for MODFLOW (Harbaugh and others, 

2000) as shown in figure 2. Seeps and springs were not 

simulated in the Central Oklahoma aquifer model (Mashburn 

and others, 2013). Where groundwater in alluvial deposits 

left the model on the eastern margin, drain cells were used 

instead of GHBs. The Drain Package is independent of other 

boundaries and only allows groundwater to leave the model 

based on the relative altitude of the drain and the aquifer 

head and drain conductance. Drain conductance was set 

equal to that of the aquifer at each drain cell, and the drain 

altitude was set to the land-surface altitude from a digitalelevation model of the area (U.S. Geological Survey, 2013a). 

Drain cell properties were only adjusted during initial model 

calibration.

Pumping wells are constant-rate flow boundaries and 

were simulated using the Multi-node Well (MNW) Package 

Version 2 for MODFLOW (Konikow and others, 2009), which 

allows wells to withdraw groundwater from multiple layers. 

Well depths reported to the OWRB were used with model 

layer thicknesses to determine from which layer each well 

pumped water (Oklahoma Water Resources Board, 2013). All 

of the well construction information was taken from Mashburn 

and others (2013) and was applied to the CPN model layers.

The land surface through which recharge flows is a 

constant-rate flow boundary. The Soil-Water-Balance (SWB) 

code of Westenbroek and others (2010) used in the Central 

Oklahoma aquifer model (Mashburn and others, 2013) was 

used in this study with increased spatial detail to estimate 

amounts and distribution of recharge to the aquifer from 

the land surface. This recharge estimate was applied to the 

CPN model using the Recharge Package for MODFLOW 

(Harbaugh and others, 2000). Mean daily recharge for 1987 

from the SWB analysis was input to the steady-state model, 

and mean daily recharge for each month from 1988 to 

2009 was input to each stress period of the transient model. 

Recharge was applied to the uppermost active layer in each 

model cell and scaled during model calibration.

Because SWB estimates recharge as deep percolation 

through the root zone of the landscape, this package does not 

provide an estimate of the seepage through the unsaturated 

zone and into bedrock. A substantial amount of the deep 

percolation can flow laterally through the soil and discharge 

into streams. Thus, where the Central Oklahoma aquifer 

was the top active model layer, recharge was decreased to 

adjust for losses between the land surface and the bedrock. 

During initial calibration, the scaling factor for recharge 

where the Central Oklahoma aquifer was the top active layer 

was determined to be 0.2, which was an estimate of the 

general fraction of the soil deep percolation that may enter 

the Central Oklahoma aquifer as recharge. Because of the 

spatial variability of the Central Oklahoma aquifer and plant 

evapotranspiration, recharge was further scaled locally during 

calibration, as described in the “Model Calibration Methods” 

section of this report.

Model Calibration Methods

Model calibration is the process by which certain model 

inputs are adjusted so that simulated potentiometric heads 

and base flows more closely match measured or estimated 

potentiometric heads and base flows. The difference between 

simulated heads and measured heads or base flows and 

estimated base flows is referred to as the residual. The residual 

is calculated as measured or estimated value minus the 

simulated value. 

The initial aquifer hydraulic properties, boundary-control 

parameters, and recharge were adjusted manually using trial 

and error to minimize head residuals in the steady-state model. 

When the head residuals were manually minimized as much 

as possible, parameters were further adjusted using automated 
Citizen Potawatomi Nation Numerical Groundwater-Flow Model 9

parameter estimation that also calibrated to estimated stream 

base-flow targets. This section describes inputs that were 

adjusted, manual and automated approaches to parameter 

estimation, and comparison of the observations to the final 

calibrated simulation.

Observations

To ensure that the groundwater-flow simulation 

accurately reproduced natural conditions and processes in 

the study area, simulation results were compared to field 

measurements or estimated values at specific locations 

(observations), also referred to as targets or calibration targets. 

For the study described in this report, observations include 

groundwater heads measured in wells retrieved from the 

National Water Information System (U.S. Geological Survey, 

2012) and estimated base flows for the North Canadian River 

as described in the “Stream Base Flow” section of this report. 

Groundwater-head observations were placed in three groups: 

(1) 149 water-level measurements used for the steady-state 

model, (2) 126 Central Oklahoma aquifer water levels 

measured during the transient period, and (3) 5 water-level 

measurements in the North Canadian River alluvial aquifer 

during the transient period. Of the 149 head observations for 

the steady-state model, 9 observations were from wells open 

to layer 1, 121 were from wells open to layer 2, 18 were from 

wells open to layer 3, and 1 was from a well open to layer 4. 

Similar to the steady-state model, head observations for 

the transient period were mostly available from layer 2 with 

a total of 97 head observations, 4 from wells open to layer 1, 

80 from wells open to layer 2, 12 from wells open to layer 

3, and 1 from a well open to layer 4. Many of the wells were 

measured multiple times with 5 of the head measurements 

being from wells open to layer 1, 113 from wells open to 

layer 2, 12 from wells open to layer 3, and 1 from a well 

open to layer 4. There was also a temporal bias with most 

of the measurements made as part of synoptic water-level 

measurements during the steady-state period of 1987 and 

during the last year of the transient model period of 2009. 

Only 35 of the 131 measurements were made between 1987 

and 2009.

Estimated base flows of the North Canadian River 

consisted of 528 observations with 264 values each being used 

for the streamflow-gaging stations at Harrah and Shawnee. 

Estimated flows represented monthly mean streamflow values 

from 1988 through 2009 for the transient model.

Initial Trial-and-Error Calibration

An initial trial-and-error calibration was performed on 

the steady-state and transient models to estimate values and 

ranges for model parameters and help guide and constrain the 

automated parameter estimation. During the initial calibration, 

parameters for Kh, specific yield (Sy; the ratio of unit of water 

produced per unit of drawdown) in the uppermost active layer, 

specific storage in lower layers, and recharge were changed 

manually on a cell-by-cell basis to reduce head-observation 

residuals. Vertical hydraulic conductivity anisotropy (the ratio 

of horizontal to vertical hydraulic conductivity) for layers 

2 and 3 was adjusted as a single value for each layer. The 

initial trial-and-error calibration reduced head residuals, and 

both head residuals and base-flow residuals were expected to 

improve during the automated parameter estimation.

From the initial calibration, the recharge entering the 

Central Oklahoma aquifer was determined to be highly 

variable across the aquifer. Instead of changing recharge 

estimates from SWB directly for each stress period, recharge 

arrays for all stress periods in the steady-state and transient 

models were adjusted using a multiplier array that included 16 

different zones (12 of which are shown on fig. 3) with different 

multiplier values. This array preserved the seasonal recharge 

variability and some of the spatial recharge variability. The 

numerous recharge multiplier zones for part of the CPN model 

area in figure 3 show the spatial variability of the Central 

Oklahoma aquifer in contrast to the larger multiplier zones 

in the North Canadian alluvial aquifer. Multiplier zones were 

set up independently in the Central Oklahoma aquifer and 

the North Canadian alluvial aquifer to accommodate much 

higher recharge rates in the alluvial deposits. Each multiplier 

zone was adjusted to reduce residuals at groundwater-head 

observation locations and improve the calibration. When 

only minor improvement in head residuals was achieved 

by manually changing parameters between model runs, 

the manual calibration was terminated and the automated 

parameter estimation was used.

Parameters

Though a model simulation might consist of thousands 

of input values to represent groundwater flow in an area, 

only some of those inputs are adjusted to calibrate a model; 

these inputs are called parameters. After initial manual trialand-error calibration, automated parameter estimation was 

used to further calibrate the model. For automated parameter 

estimation, the recharge multiplier zones used in the initial 

trial-and-error calibration were used and new parameters 

for Kh and streambed hydraulic conductivity were added. 

Vertical hydraulic conductivity anisotropy values for layers 

2 and 3 were included as parameters. The same parameters 

were used for the steady-state and transient groundwater-flow 

models. For the sensitivity analysis, some parameters were 

also included in parameter groups. Recharge multiplier zone 

parameters in the initial calibration were included in the group 

transrch.
10 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

!(

!( !(

!( !(

!(

!( !( !( !( !( !( !(

!(

!( !(

!( !(

!(

!( !( !(

!(

!( !(

!( !( !(

!(

!( !(

!( !(

!(

!(

!(

!( !(

!(

!( !( !( !( !( !(

!(

!(

!(

!( !(

!(

!( !( !(

!( !( !(

!( !(

!(

!( !( !( !(

!(

!( !( !( !( !( !(

!( !( !( !( !( !( !(

!(

!( !( !( !( !( !( !( !(

!( !( !( !( !( !( !( !(

!( !( !( !( !( !( !( !( !(

!( !( !( !( !( !( !( !( !(

!( !( !( !( !( !( !( !( !(

!( !( !( !( !( !( !( !( !(

!( !( !( !(

!( !( !( !( !( !(

!( !( !( !( !( !( !(

!( !( !( !( !( !( !(

!( !( !( !( !( !( !( !(

!(

!( !( !( !( !( !( !( !(

!(

!( !( !( !( !( !( !( !( !(

!( !( !( !( !( !( !( !( !(

!( !( !( !( !(

!( !( !( !( !( !(

!( !( !( !( !( !(

!( !( !( !( !( !(

!(

!( !( !( !( !( !( !(

!( !( !( !( !( !( !(

10

11

11

11

12

12

12

13 2

2

2

2

2

3

4

4

5

5

5

5

5

5

5

5

6

6

6

6

6

6

6

6

6

7

7

7

7

7

7

7

7

7

7

7

8

8

8

8

8

8

8

8

8

9

9

9

9

3

EXPLANATION

Recharge multiplier zone

 2

 3

 4

 5

 6

 7

SFR2 streambed hydraulic 

 conductivity parameter

 hc-1

 hc-2

 hc-3

North Canadian River alluvial

 aquifer extent

North Canadian River

Horizontal hydraulic conductivity pilot

 point, by layer (parameter group)

 1 (klay1)

 2 (klay2)

 3 (klay3)

 4 (klay4)

!(

!(

!(

!(

8

9

10

11

12

13

SFR2, Streamflow-routing package version 2

Area enlarged

Model active area

97°14' 97°12' 97°10' 97°08' 97°06' 97°04' 97°02'

35°30'

35°28'

35°26'

35°24'

35°22'

Base modified from U.S. Geological Survey digital data

Universal Transverse Mercator, zone 14

North American Datum of 1983 0 1 2 3 4 MILES

0 1 2 3 4 KILOMETERS

Figure 3. Recharge multiplier zones, horizontal hydraulic conductivity pilot points, and streambed hydraulic conductivity parameters 

used in the Citizen Potawatomi Nation numerical groundwater-flow model, Citizen Potawatomi Nation Tribal Jurisdictional Area, central 

Oklahoma.
Citizen Potawatomi Nation Numerical Groundwater-Flow Model 11

Streambed hydraulic conductivity in the SFR2 Package 

was divided into four parameters (hc-1 through hc-4) for four 

reaches of the North Canadian River; parts of hc-1 through 

hc-3 are shown on figure 3. Parameter hc-1 included the North 

Canadian River in the model active area upstream from the 

streamflow-gaging station at Harrah (USGS station 07241520) 

(fig.1). Parameter hc-2 included the reach between hc-1 and 

the river near Shawnee Reservoir 2. Parameter hc-3 included 

the reach of the North Canadian River between hc-2 and just 

upstream from the streamflow-gaging station at Shawnee 

(USGS station 07241800), and parameter hc-4 included the 

reach of river in the model active area downstream from hc-3. 

All stream reaches in each parameter segment were assigned 

the same value for streambed hydraulic conductivity and were 

adjusted by the same amount during each iteration of the 

automated parameter estimation process. The four streambed 

hydraulic conductivity parameters were included in the model 

parameter group sfrhc.

The adjustment of Kh on a cell-by-cell basis in automated 

parameter estimation, as was done in the initial manual 

calibration, was not feasible because doing so would require 

determination of a Kh parameter at every cell. Although 

estimating Kh for large groups of cells or zones is feasible, 

such an approach can result in large changes and unrealistic 

contrasts in Kh at zone boundaries; also, there may not be 

sufficient data to delineate such zones. A preferable approach 

is to use pilot points, which are discrete points located 

within the active model cells, in this case grouped by each 

layer; Kh was estimated for each pilot point individually and 

subsequently interpolated to every active model cell using the 

kriging spatial interpolation algorithm (Doherty, 2010). Pilot 

points were placed where points used in Mashburn and others 

(2013) caused a distinct pattern in the Kh values for each 

layer and were manually placed where they were near inactive 

model cells. A subset of the pilot points symbolized by layer is 

shown in figure 3. This approach produces smooth variations 

in Kh during calibration and thus introduces fewer artifacts 

to the simulation. This procedure used a total of 1,473 pilot 

points that were assigned to four parameter groups by layer: 

176 points for layer 1 in group klay1, 566 for layer 2 in group

klay2, 448 for layer 3 in group klay3, and 283 for layer 4 in 

group klay4.

Automated Parameter Estimation

Automated parameter estimation was undertaken 

using the statistical techniques provided in the PEST suite 

of software (Doherty, 2010) after the initial trial-and-error 

calibration. PEST is a model-independent tool that at its 

simplest level, automatically adjusts parameters within 

predefined limits, reruns the model, checks the calibration, 

and continues until the calibration cannot be improved further. 

Complex algorithms integral to PEST allow a model to be 

improved efficiently and quickly, allowing a large number 

of parameters to be used, and enabling a more accurate 

representation of the intrinsic variability of natural hydrologic 

systems.

In addition to the efficiencies gained through automation 

and powerful statistical approaches, PEST also facilitates 

Tikhonov regularization, wherein a penalty is applied to the 

objective function when parameters deviate from their initial 

values (Doherty, 2010). The result of that regularization 

is a more stable parameter estimation process that allows 

for introduction of variability when necessary to improve 

calibration but preserves initial values that may be based 

upon the best data and reasoning available to the investigator. 

For this study, Tikhonov regularization was applied to 

Kh parameters of aquifer layers 1–4, with initial values 

determined during the trial-and-error calibration. Parameters 

for hydraulic properties were limited by end-member values 

estimated for aquifer properties in Mashburn and others 

(2013). This limiting of parameters prevented the automated 

process from choosing values that are not realistic for aquifers 

in this study to achieve a best match with observations.

Calibration Results

Figure 2 shows the distribution of steady-state head 

residuals for all observation wells and the associated model 

layer. The mean residual for groundwater heads in the steadystate model was -1.1 ft (table 1), indicating that on average, 

simulated groundwater heads were 1.1 ft higher than measured 

groundwater heads. For transient groundwater heads in the 

Central Oklahoma aquifer, the mean residual was -1.7 ft, 

indicating that on average, simulated groundwater heads 

were 1.7 ft higher than measured groundwater heads. The 

mean residual was -4.1 ft for transient groundwater heads 

in the North Canadian River alluvial aquifer (table 1). The 

root-mean-squared groundwater head residuals for steadystate, transient Central Oklahoma aquifer, and transient North 

Canadian River alluvial aquifer had a maximum value of 

about 3 percent of the total groundwater-level relief (450 ft) 

in the study area, indicating that errors or inaccuracies were 

only a small part of the total model response (Anderson and 

Woessner, 1992). Figure 4 shows the correlation between 

simulated and measured head observations. Scatter about 

the regression line was caused by local variations in the 

Central Oklahoma aquifer not included in the model, but 

there is a reasonable general correlation with a coefficient 

of determination (R2

) value of 0.98. There was no apparent 

characteristic of the simulated head accuracy associated with 

higher or lower altitudes.
12 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

Table 1. Statistical comparison of measured and simulated groundwater heads and stream base flows, numerical groundwater-flow 

model for the Citizen Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma.

[Residual is calculated as the measured or estimated value minus the simulated value; thus, a negative residual for water levels or streamflows indicates that 

simulated values are larger than estimated or measured values]

Observation group

Steady-state 

groundwater head

(feet)

Transient groundwater head 

(feet) Transient 

streamflows 

(cubic feet per second) Central Oklahoma 

aquifer

North Canadian River

alluvial aquifer

Observation count 149 126 5 528

Minimum residual -39.0 -45.7 -12.1 -2,853

Mean residual -1.1 -1.7 -4.1 74

Maximum residual 44.7 28.8 7.3 9,122

Mean absolute residual 9.4 6.3 7.0 204

Root-mean-squared residual 13.0 9.0 7.6 576

Measured head, in feet

850 900 950 1,000 1,050 1,100 1,150 1,200 1,250

Simulated head, in feet

850

900

950

1,000

1,050

1,100

1,150

1,200

1,250

Simulated head = 0.9833 (measured head) + 19.575

Coefficient of determination (R2

) = 0.98

Linear regression (transient 

 groundwater head)

Transient groundwater head

EXPLANATION

Figure 4. Transient measured groundwater heads plotted with simulated heads for the Citizen Potawatomi numerical groundwaterflow model, Citizen Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma.
Citizen Potawatomi Nation Numerical Groundwater-Flow Model 13

The CPN model provided reasonable matches between 

simulated stream base flow and estimated stream base flow 

at the Harrah streamflow-gaging station (fig. 5A) and the 

Shawnee streamflow-gaging station (fig. 5B). The mean 

estimated stream base flow used in the transient model was 

604 ft3

/s. The mean residual for transient stream base flow 

was 74 ft3/s (table 1), indicating that on average, simulated 

stream base flow was 74 ft3/s, or about 12 percent, less than 

the mean estimated stream base flow. This difference was 

caused, at least in part, by runoff peaks coming into the model 

from upstream and no routing of runoff within the model 

area. The numerical flow model only affects the streamflow 

routed through the model by changes in base flow, which is 

a relatively small fraction of streamflow during peak-flow 

periods.

A statistical summary of the final calibrated values 

for Kh, Sy, and vertical anisotropy estimated during model 

calibration is listed in table 2, and the numerical groundwaterflow model annual flow budget and total flow for 1987–2009 

are listed in table 3. The annual model flow budget includes 

the net flow—inflow minus outflow—for principal model 

boundaries, with positive flow indicating groundwater 

entering the model and negative flow leaving the model. The 

largest inflow was from recharge, which was nearly equal to 

the largest discharge, stream base flow (table 3). Negative 

volumes for change in storage represent adding groundwater 

to storage during wet years, and positive values represent 

water taken from storage during drier years. 

Model Sensitivity

The PEST suite of software adjusts all parameters 

individually and determines the relative influence of 

parameters and parameter groups on the head and flow model 

observations. Figure 6 presents the relative magnitude of 

changes in residuals caused by a 1-percent change in the value 

of each parameter estimated during the initial trial-and-error 

calibration. 

Changes in parameters for group sfrhc caused the 

largest change of residuals because of the direct influence 

on the residuals for streamflows in the North Canadian 

River, which had the largest contribution to the parameter 

estimation objective function. However, sensitivity to sfrhc

was only approximately twice that of the least sensitive 

parameter group klay1. Parameter group klay1 was the group 

that represented the North Canadian River alluvial aquifer 

hydraulic conductivity. This parameter group has more effect 

on head than on streamflow, but there were only nine head 

observations in layer 1 in the steady-state model and five in 

the transient model. Thus, changes in klay1 may have caused 

a greater effect on residuals if there had been more head 

observations.

Simulated heads were sensitive to changes in parameter 

groups klay2–4 and transrch. Because layer 2 was the 

uppermost active layer in the active model area that did not 

include alluvial deposits, it received recharge, and the majority 

of head observations in layer 2 were affected by transrch as 

well as klay2. Base flow was affected by transrch in areas 

where recharge was applied to the North Canadian alluvial 

aquifer.

The parameters to which model observations were most 

sensitive included streambed hydraulic conductivity, Kh, and 

recharge. Streambed hydraulic conductivity was limited to the 

North Canadian River overlying the North Canadian alluvial 

aquifer. Thus, model observations in the Central Oklahoma 

aquifer were only responding to changes in klay2–4 and 

transrch.
14 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

1988198919901991199219931994199519961997199819992000200120022003200420052006200720082009

19881989199019911992199319941995

19961997199819992000

20012002200320042005200620072008

2009

Date

10

100

1,000

10,000

10

100

1,000

10,000

100,000

Streamflow, in cubic feet per second

B. North Canadian River at Shawnee, Oklahoma

A. North Canadian River near Harrah, Oklahoma

EXPLANATION

Simulated

Estimated

Figure 5. Simulated and estimated streamflows for North Canadian River streamflow-gaging stations A. near Harrah, Oklahoma, and 

B. at Shawnee, Okla., Citizen Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma.
Citizen Potawatomi Nation Numerical Groundwater-Flow Model 15

Table 2. Statistical summary for calibrated hydraulic parameters.

Layer Mean Minimum Maximum Standard 

deviation

Horizontal hydraulic conductivity 

(feet per day)

1 43.42 0.04 315.53 151.39

2 11.90 0.03 82.02 14.61

3 1.79 0.03 16.40 1.20

4 1.18 0.03 16.40 1.37

Specific yield 1 0.19 0.11 0.20 0.02

2 0.12 0.07 0.20 0.01

3 0.11 0.07 0.16 0.01

4 0.10 0.10 0.10 0.00

Vertical anisotropy (vertical/

horizontal hydraulic 

conductivity)

1 19 4 20 67

2 5 1 10 4

3 516 100 960 395

4 4 1 8 4

Table 3. Numerical groundwater-flow model annual flow budget, Citizen Potawatomi Nation Tribal Jurisdictional Area, central 

Oklahoma.

[1987*, is from the steady-state model; Positive flow is to the model and negative flow is from the model; all flows are net flow (balance of inflow and outflow) 

rounded to 100 acre-feet] 

Annual model flow in acre-feet

Year Recharge Change in 

storage

General head 

boundaries 

and drains

Stream 

base flow Lake Wells

1987* 170,000 0 4,000 -172,000 0 -2,000

1988 334,700 -60,900 3,500 -249,900 -24,700 -2,800

1989 306,600 -4,900 2,600 -285,300 -15,800 -3,100

1990 480,400 -114,700 1,400 -348,600 -15,400 -3,000

1991 452,600 -92,400 1,400 -344,100 -14,100 -3,400

1992 522,400 -97,300 -300 -406,500 -15,200 -3,000

1993 474,300 -33,000 -1,700 -420,500 -15,800 -3,300

1994 378,300 -3,700 -700 -356,600 -14,400 -2,900

1995 437,100 -30,400 -1,400 -385,900 -16,200 -3,300

1996 436,900 -46,200 -1,000 -373,200 -13,000 -3,500

1997 452,500 -44,500 -1,500 -390,800 -12,200 -3,400

1998 390,400 10,100 -1,600 -385,600 -10,600 -2,600

1999 331,700 34,100 -1,000 -348,000 -13,100 -3,700

2000 386,900 -19,900 -700 -348,900 -14,300 -3,100

2001 374,800 4,600 -1,000 -358,800 -15,800 -3,800

2002 381,000 -12,600 -700 -348,800 -15,600 -3,400

2003 183,900 111,800 200 -276,700 -16,100 -3,200

2004 423,600 -65,700 0 -337,700 -16,900 -3,300

2005 236,200 89,200 -100 -305,700 -16,300 -3,300

2006 204,000 60,200 1,300 -245,700 -16,400 -3,400

2007 550,900 -131,800 -500 -391,800 -22,400 -4,400

2008 298,100 53,400 -500 -329,200 -18,400 -3,500

2009 381,200 -25,700 -100 -332,000 -19,900 -3,600

Total 8,588,500 -420,300 1,600 -7,742,300 -352,600 -75,100
16 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

Inset Numerical Groundwater-Flow 

Model

A numerical groundwater-flow model was required to 

optimize groundwater pumping at two EDZs, the geothermal 

supply area and the Iron Horse Industrial Park just south of 

the city of Shawnee, Okla. (fig. 7). Because of the small area 

and increased detail required to analyze water resources in 

the EDZs, a detailed transient inset model with smaller cell 

sizes was constructed using the calibrated model inputs and 

layers from the CPN model. The inset model was constructed 

completely within the CPN model active area, on the North 

Canadian River near the northern margin (fig. 2) spanning 

6.5 mi east to west and 5 mi north to south (fig. 7). The 

analysis performed using the inset model predominantly 

affected the North Canadian River alluvial aquifer, and the 

extent of the inset model is well beyond the extent of the 

alluvium to the north and south, and several stream miles 

upstream and downstream to eliminate boundary effects on the 

analysis. 

The inset model also included heads calculated in the 

CPN model, and these heads were used in GHBs along the 

margin of the inset model. Like the CPN model, the inset 

model used the SFR2 Package to simulate routed streamflow, 

the MNW Package for wells, and the Recharge Package. 

The inset model provided a more efficient analysis of the 

optimized pumping and allowed simulation of local pumping 

effects on groundwater heads and North Canadian River 

streamflow depletion.

To provide more detail, the inset model cells were 

reduced to 164 ft on all sides from the 820-ft cells used 

in the CPN model. Properties such as recharge and Kh 

were transferred from the larger CPN model cells to the 

corresponding smaller cells. Pumping rates for groundwaterwithdrawal wells in the inset model area were taken from the 

CPN model without further refinement. Simulated streamflow 

in the North Canadian River and tributaries in the CPN model 

were used as inflow to those streams that originated outside 

the inset model and crossed into the inset model active area 

(fig. 7). The inset model was not calibrated separately, and 

there were no head observations available within the inset 

model in the North Canadian River alluvial aquifer during the 

transient simulation period.

The starting head for all layers of the inset model was set 

equal to that of the CPN model. The external boundary of the 

inset model was composed of GHB cells with the head set to 

the starting head of the CPN model at each cell location. The 

head in GHB cells remained constant for the duration of the 

transient model simulation.

Optimization of Groundwater 

Withdrawals

Two EDZs along the North Canadian River were selected 

by the CPN for groundwater development: the geothermal 

supply area and the Iron Horse Industrial Park (fig. 7). Both 

areas are along the North Canadian River and are separated by 

about 2 mi. Both areas are underlain by the North Canadian 

River alluvial aquifer, and the groundwater resources in those 

areas are pumped from that aquifer. The areas are sufficiently 

distant from each other, and pumping at one site is assumed 

to have negligible effect on heads or water availability at the 

other.

Parameter group

klay1 klay2 klay3 klay4 sfrhc transrch

Mean of dimensionless scaled sensitivity

4.5

4.0

3.5

3.0

2.5

2.0

1.5

1.0

0.5

0

Hydraulic conductivity for layer 1 (klay1)

EXPLANATION

Hydraulic conductivity for layer 2 (klay2)

Hydraulic conductivity for layer 3 (klay3)

Hydraulic conductivity for layer 4 (klay4)

Streambed hydraulic conductivity (sfrhc)

Transient recharge multiplier (transrch)

Figure 6. Scaled sensitivity of model observations in parameter groups, Citizen Potawatomi Nation numerical groundwater-flow model.
Optimization of Groundwater Withdrawals 17

EXPLANATION

Economic development zone

 boundary

Pumping well

Observation well

Numerical groundwater-flow model

 Top active layer

 Alluvial aquifer

 Central Oklahoma aquifer

 Boundary condition, by cell type

 General head boundary cell

 SFR2 stream cell

SFR2, Streamflow-routing package version 2

Base modified from U.S. Geological Survey digital data

Universal Transverse Mercator, zone 14

North American Datum of 1983

270

HARDESTY ROAD

E WALNUT STREET

0 0.5 1.0 1.5 2.0 MILES

0 0.5 1.0 1.5 2.0 KILOMETERS

North Canadian River

Iron Horse

Industrial

Park 

Geothermal Geothermal supply area

96°58' 96°57' 96°56' 96°55' 96°54' 96°53' 96°52'

35°21'

35°20'

35°19'

35°18'

35°17'

35°16'

Shawnee

Figure 7. Inset model active area, top active layer, and boundary condition cells, and locations of economic-development zones 

included in the optimization analysis, Citizen Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma.
18 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

Water Demands

The geothermal supply area required a continuous 

flow of 500 gallons per minute (gal/min) from one or more 

wells located within the geothermal supply area, just south 

of the North Canadian River (fig. 7). Plans for expansion of 

the Iron Horse Industrial Park have not set a specific water 

requirement, so the objective for that EDZ was to determine 

the total available groundwater at the site. Water quality 

was also an issue of concern for the Iron Horse Industrial 

Park, and stream water in this area has been found to contain 

more dissolved solids and chlorides, nitrogen, and organic 

compounds than the groundwater in alluvial deposits (Becker, 

2014); thus, the amount of stream water entering the Iron 

Horse Industrial Park well field was a concern.

Streamflow Depletion by Wells

Streamflow depletion by a well takes place when a 

well reduces streamflow either by capturing groundwater 

that would otherwise discharge to the stream as base flow 

or by inducing streamflow to flow into the adjoining aquifer 

through induced infiltration, or both (Barlow and Leake, 

2012). At a gaining reach, the aquifer head is higher than 

the stream stage, but if the aquifer head is drawn down 

below the stream stage by pumping, the head gradient is 

reversed, resulting in induced infiltration of streamflow to 

groundwater.

The SFR2 Package was used to determine the rates 

of base flow and streamflow depletion through induced 

infiltration of North Canadian River streamflow that 

could enter wells. The SFR2 Package determines rates of 

groundwater flow to or from each stream reach; a reach is 

defined as the length of stream channel in each model cell. 

The SFR2 results showed that before development of the 

Iron Horse Industrial Park well field, all adjacent North 

Canadian River stream reaches were gaining base flow during 

model simulations. If pumping caused gaining reaches to 

become losing reaches, the volume of induced infiltration 

of streamflow to the aquifer during well pumping was 

calculated using SFR2. Because the North Canadian River is 

a hydrologic groundwater boundary in the North Canadian 

alluvial aquifer and there were no other known groundwater 

withdrawals in the area, all water lost from the stream was 

assumed to enter the pumping wells. However, because the 

fate of induced infiltration was not confirmed in this study, this 

calculated flow is considered to be a maximum value.

As shown in figure 5, streamflow in the North Canadian 

River is highly variable. Flow in the reach adjacent to the Iron 

Horse Industrial Park was tracked with the SFR2 Package. 

The median streamflow in the North Canadian River adjacent 

to the Iron Horse Industrial Park was 302 ft3

/s during the 

transient model period.

Optimization Methods

Pumping optimization was performed by coupling the 

inset numerical groundwater-flow model with the most recent 

version of the USGS Groundwater-Management Process 

(GWM) for MODFLOW (Ahlfeld and others, 2005), which 

is referred to as GWM-VI, version 1.0.1 (Banta and Ahlfeld, 

2013). Optimization determined whether the 500 gal/min 

required for the geothermal supply could be pumped from 

the geothermal supply area, and the maximum amount of 

water that could be pumped from wells in the Iron Horse 

Industrial Park. All of the pumping rates were constrained by 

limits on head drawdown to protect the alluvial aquifer from 

dewatering.

The GWM process uses managed wells, which are 

planned or existing wells, and determines pumping rates 

during optimization as flow variables. The optimization 

process determines the response at each constraint to pumping 

at each managed well independently and then uses an 

objective function method to determine the optimal pumping 

rate at each well that satisfies all constraints. The objective 

function at the geothermal supply area included a criterion 

that the total pumping from managed wells had to be at least 

500 gal/min. At the Iron Horse Industrial Park, the maximum 

pumping rate that satisfied all head drawdown constraints was 

determined.

Optimization Design

Because of the presence of the North Canadian River 

adjacent to the well fields and the relatively thin alluvial 

aquifer, the response of heads to pumping was not entirely 

linear. To correct for this nonlinearity, the sequential linear 

programming algorithm in GWM was used to linearize 

nonlinear processes and calculate the response at constraints 

(Ahlfeld and others, 2005).

The mean thickness of the alluvial aquifer is 

approximately 43 ft in the parts of the EDZs where 

optimization was performed and the drawdown in managed 

wells could cause loss of well yield. The MNW package 

calculates loss of yield and adjusts the flow rate to simulate the 

net flow from the well. To prevent loss of yield, preliminary 

model runs were used to determine the pumping rate at which 

loss of yield takes place at each well without regard to other 

constraints. This flow rate was used in GWM as the maximum 

allowable pumping rate for each managed well.

Pumping rates at all managed wells were considered to 

be the optimal continuous flow rates with all wells pumping 

simultaneously. The drawdown and effects on streamflow 

in the North Canadian River were evaluated over the entire 

transient model period of 21 years; thus, the pumping rates 

were conservative and will be less than short-term pumping 

rates that could be achieved by cycling pumping rates from 

several wells.
Effects of Prolonged Drought 19

Managed Wells

Five proposed managed wells were placed in the 

geothermal supply area (GTW1–5; fig. 8). A total of 16 

managed wells were placed in the northern part of the Iron 

Horse Industrial Park (IHW01–16; fig. 9). A minimum 

well spacing of 650 ft was chosen based on the estimated 

drawdown cone at each pumping well observed during 

preliminary model runs. Similar to recently developed wells in 

the local area, all of the planned wells were open only to the 

North Canadian alluvial aquifer (layer 1) and did not penetrate 

the Central Oklahoma aquifer.

Pumping Constraints

The principal concern at the geothermal supply area 

and the Iron Horse Industrial Park was that head drawdown 

would deplete groundwater in the alluvial aquifer and could 

be unsustainable. Head constraints were placed on all wells 

to limit the head drawdown. The maximum head drawdown 

at each of the geothermal supply managed wells was set at 

16.4 ft, which was estimated to be the maximum drawdown 

that would not cause loss of well yield. The head constraints 

at the Iron Horse Industrial Park well field allowed the head 

to drop as much as one-half of the predevelopment saturated 

thickness (starting head relative to the base of layer 1) at each 

well.

Streamflow depletion was expected to take place under 

pumping conditions at both EDZs. However, preliminary 

runs of the inset model showed that streamflow depletion and 

induced infiltration were a very small fraction of the total 

streamflow in the North Canadian River and total pumping. 

Thus, the streamflow depletion was not used as a pumping 

constraint, but induced infiltration was calculated to estimate 

the total streamflow that could end up in the pumped water at 

the Iron Horse Industrial Park well field.

Optimization Results

Optimized pumping rates for wells at the geothermal 

supply area ranged from 98 to 201 gal/min (fig. 8). The 

combined optimized pumping rate was approximately 

638 gal/min, which exceeded the goal of 500 gal/min pumping 

for the five managed wells. Managed wells located in areas 

where the aquifer is relatively thin or little water is taken from 

the North Canadian River had lower optimized pumping rates. 

The total pumping rate might be increased by increasing well 

spacing or adding wells outside of the well field, or both, 

although these changes would require the wells to be located 

outside the geothermal EDZ boundary.

Optimized pumping rates for the 16 managed wells at 

the Iron Horse Industrial Park ranged from 62 to 141 gal/min 

(fig. 9). The total estimated optimized pumping rate at the Iron 

Horse Industrial Park was estimated to be 1,472 gal/min. Wells 

closer to the North Canadian River and those to the southwest 

generally produced more water as shown in figure 9, in part 

because the stream provided a hydrologic boundary that 

reduced head drawdown under pumping conditions, and the 

aquifer is slightly thicker to the southwest. 

Under pumping conditions, the base flow to the North 

Canadian River adjacent to the Iron Horse well field was 

decreased, and groundwater head declined to the point 

where streamflow began to infiltrate to groundwater. Flow 

into and out of each stream reach defined in the inset model 

was calculated by the SFR2 package, and rates of induced 

infiltration that resulted from the optimal pumping scheme 

at the Iron Horse Industrial Park are shown for each stream 

reach numbered upstream to downstream on figure 9. The total 

rate of induced infiltration from the North Canadian River for 

the optimal pumping scheme was approximately 4.1 gal/min, 

which was approximately 0.3 percent of the total pumping of 

1,472 gal/min.

A graph of the distribution of predevelopment base flow 

and induced infiltration from each inset model reach adjacent 

to the Iron Horse Industrial Park is shown in figure 10. The 

reaches are shown in map view on figure 9. This figure shows 

the change in base flow and amount of induced infiltration 

in each short section of the river that may be captured by 

the Iron Horse wells. The green line in figure 10 shows 

the predevelopment base flow for each stream reach. After 

development, the yellow line shows that all reaches have 

induced infiltration caused by pumping. The greatest loss 

of water occurred in reaches 27 and 29 (fig. 10). The wells 

closest to the reaches that lost the most water were IHW01, 

03, 05, and 06 (fig. 9).

Effects of Prolonged Drought

The effects of drought on a groundwater system can be 

manifested in changes in inflow from decreased precipitationrelated recharge, which can result in decreased groundwater 

discharge to streams and wells. Decreases in recharge also 

can cause a deficit that is balanced by loss of groundwater 

from aquifer storage. To estimate the effects of drought on 

the streamflow and groundwater storage, the CPN model was 

altered to simulate estimated drought conditions.

The CPN Tribal Jurisdictional Area in central 

Oklahoma has had several severe hydrologic droughts since 

recordkeeping began in 1925, including the Dust Bowl 

drought (1929–41) and the droughts of 1952–6, 1961–72, 

and 1976–81 (Shivers and Andrews, 2013). The duration of 

these drought periods ranged from 4 years in the 1950s to 

the 13-yearlong Dust Bowl. The period of this study does 

not include any of these historical droughts, so hypothetical 

drought conditions were estimated based on hydrologic 

characteristics of recorded droughts.
20 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

Geothermal supply area

35°18'45"

35°18'30"

35°18'15"

96°56'45" 96°56'30" 96°56'15" 96°56'00" 96°55'45"

Note: Total estimated rate of optimized pumping was 638 gallons per minute

GTW1

GTW2

GTW3

GTW4 GTW5

201

118

116

98 103

Economic development zone—Geothermal supply area

 (location shown on fig. 6)

Numerical groundwater-flow model cell for managed

 well at the geothermal supply area, by optimized

 pumping rate, in gallons per minute—Well number

 and pumping rate are shown next to cell

 99

 100 to 104

 105 to 119

 120 to 202

EXPLANATION

Base modified from U.S. Geological Survey digital data

Universal Transverse Mercator, zone 14

North American Datum of 1983

North Canadian River

177

HARDESTY ROAD

HARDESTY DRIVE

RANGELINE ROAD

0 500 1,000 1,500 2,000 FEET

0 100 200 300 400 METERS

Figure 8. Groundwater Management Process managed wells with optimized pumping rates at the geothermal supply area, Citizen 

Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma.
Effects of Prolonged Drought 21

Iron Horse Industrial Park

96°53'45" 96°53'30" 96°53'15" 96°53'00" 96°52'45"

35°19'15"

35°19'00"

35°18'45"

35°18'30"

IHW06

IHW05

IHW10

IHW04

IHW02

IHW03

IHW01

IHW09

IHW16

IHW07 IHW08

IHW11 IHW13 IHW14

IHW12

IHW15

110

99

79

73

109

136

62

73

71 81

103 100 82

141

78

75

17

33

27

49

43

29

24

35

19 20

23

48

44

18

40

46

41

31

34

8

21

32

39 38

30

37

28

10

13

45

11

47

14

25 26

16 15

12

22

9

42

36

Note: Total estimated rate of optimized pumping was 1,472 gallons per minute.

 Total rate of induced infiltration from the North Canadian River for optimal 

 pumping scheme was 4.3 gallons per minute (0.3 percent of pumping).

Base modified from U.S. Geological Survey digital data

Universal Transverse Mercator, zone 14

North American Datum of 1983

0 500 1,000 1,500 2,000 FEET

0 100 200 300 400 METERS

EXPLANATION

Economic development zone—Iron Horse Industrial Park

 (location shown on fig. 6)

Numerical groundwater-flow model cell for managed well

 at the Iron Horse Industrial Park, by optimized pumping 

 rate, in gallons per minute—Well number and pumping 

 rate are shown next to cell

 62

 63 to 75

 76 to 82

 83 to 110

 111 to 141

Stream reach and stream reach number, by induced infiltration from North

 Canadian River for optimal pumping scheme, in cubic feet per second 

0 to 0.006

0.007 to 0.014

0.015 to 0.021

0.022 to 0.032

0.033 to 0.048

40

22

15

8

25

Figure 9. Groundwater Management Process managed wells, optimized pumping rates, and resulting effect on induced infiltration by 

stream reach number at the Iron Horse Industrial Park, Citizen Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma.
22 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

Years in the historical record that had the largest negative 

departure from the median precipitation were 1953 and 2011; 

both were nearly a 50-percent decrease in annual precipitation 

(Shivers and Andrews, 2013). Because groundwater recharge 

is tied to precipitation and in many areas estimated to be a 

percentage of precipitation, drought was approximated by 

decreasing recharge over the entire CPN model for each 

stress period by 50 percent for a 10-year period from 1990 to 

2000. The results are approximate because changes in other 

flows such as runoff and streamflow entering the model and 

processes such as plant consumption were not changed during 

the hypothetical drought. 

Changes to the volume of groundwater in storage in the 

CPN Tribal Jurisdictional Area were estimated by comparing 

the model run under calibrated recharge conditions to the 

model run with hypothetical drought recharge conditions using 

the “ZONEBUDGET” program (Harbaugh, 1990). The total 

water in storage was estimated using the saturated thickness in 

the model and the Sy of the aquifer in each cell. The Central 

Oklahoma aquifer (layers 2-4) had a mean Sy of 0.11 and 

the North Canadian alluvial aquifer (layer 1) had a mean Sy 

of 0.19 after model calibration (table 2). The total calibrated 

model groundwater in storage at the end of the hypothetical 

drought period (2000) without the hypothetical drought 

conditions was approximately 148 million acre-feet (acre-ft) in 

the Central Oklahoma aquifer (table 4), and 163,100 acre-ft in 

the North Canadian River alluvial aquifer.

During the hypothetical drought, model recharge to 

the North Canadian River alluvial aquifer decreased by 

441,000 acre-ft and to the Central Oklahoma aquifer by 

696,000 acre-ft (table 4). This decrease in recharge caused 

a decrease of 351,500 acre-ft of groundwater in storage 

(14,100 acre-ft or 8.6 percent in the North Canadian alluvial 

aquifer and 346,400 acre-ft or 0.2 percent in the Central 

Oklahoma aquifer) at the end of the hypothetical drought. The 

total change was approximately 0.2 percent of all groundwater 

in storage at the end of the drought period. This volume of 

groundwater storage loss shows that the rate of recharge to the 

Central Oklahoma aquifer over a 10-year period is much less 

than the total water in storage. 

North Canadian River reach number

7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

Base flow and induced infiltration, in cubic feet per second

–0.15

–0.10

–0.05

0

0.05

0.10

EXPLANATION

Predevelopment base flow

Postdevelopment induced infiltration

North Canadian River

Figure 10. Predevelopment base flow and postdevelopment induced infiltration for the North Canadian River caused by optimized 

pumping at the Iron Horse Industrial Park, Citizen Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma.
Effects of Prolonged Drought 23

The budget for base flow in the model during the 10-year 

simulated drought decreased by 678,500 acre-ft (33 percent); 

base flow to the North Canadian River decreased by 386,500 

acre-ft (37 percent), and for all other parts of the CPN Tribal 

Jurisdictional Area base flow decreased by 292,000 acre-ft 

(28 percent). Approximately 11 percent of the decrease in 

recharge to the North Canadian alluvial aquifer, and 66 percent 

of the change in recharge to the Central Oklahoma aquifer 

were offset by change in storage and base flow to streams. 

The balance of the change in recharge was due to decreases in 

discharge to GHBs, drains, and lakes.

The CPN model also was used to simulate changes to 

streamflow in the North Canadian River routed through the 

model during and after the drought period. Immediately after 

the start of the hypothetical drought, simulated streamflow 

in the North Canadian River at the Shawnee streamflowgaging station decreased because of decreased base flow and 

continued to decrease until the end of the simulated drought in 

2000 (fig. 11). Routed streamflow slowly recovered through 

2008 as recharge restored the amount of groundwater in 

storage. The largest estimated decrease in mean streamflow 

was 28 percent during September of 1994, and the percent 

decrease approximated by a 50-month moving average was 

approximately 16 percent in 1994. Streamflow did not recover 

to near nondrought streamflow conditions until the wet year 

of 2007, which was 7 years after the end of the hypothetical 

drought. The delay in streamflow recovery was caused by the 

aquifer recharge replenishing groundwater storage that had 

been depleted during the drought.

Model Limitations

All models are based on limited amounts of data and 

thus are necessarily simplifications of actual systems. Model 

limitations are a consequence of uncertainty in three basic 

aspects of the model, including inadequacies, inaccuracies, 

or simplifications in observations used in the model; the 

representation of geologic complexity in the hydrogeologic 

framework; and the representation of the groundwater 

system in the model. It is important to understand how these 

characteristics limit the use of the model.

Layer 2 in the model had many more head observations 

than other layers, and there were many more observations in 

1987 and 2009 than in the intervening years. Base flow in the 

North Canadian River was estimated throughout the transient 

model period but there were few head observations in layer 1. 

These biases in observation introduce uncertainty in the results 

of model simulations. 

Although the CPN and inset models were refined 

from the numerical groundwater-flow model of the Central 

Oklahoma aquifer in Mashburn and others (2013), they do 

not include any more detailed hydraulic or hydrogeologic 

information, and local analyses can be improved by additional 

detailed hydrogeological characterization. The Central 

Oklahoma aquifer and alluvial aquifer hydraulic parameters 

were highly generalized and could not include the substantial 

local variations in hydraulic properties. Thus, the results of 

the simulations presented in this report are approximate mean 

values of local conditions. Furthermore, the Central Oklahoma 

aquifer is not composed of geologic layers that can be defined 

individually or grouped as hydrogeologic units, which 

introduces uncertainty in the model because model layers do 

not have unique hydraulic properties.

Recharge was applied only to layer 2 (Central Oklahoma 

aquifer) and to layer 1 (alluvial units). The Central Oklahoma 

aquifer was calibrated to head observations using Kh for all 

layers and recharge. Kh and recharge are correlated, which 

limits the ability to uniquely estimate parameter values 

because the same result could be achieved with multiple 

combinations of the parameter values. The North Canadian 

alluvial aquifer was calibrated to head observations and baseflow estimates using recharge, Kh, and streambed hydraulic 

conductivity parameters. However, the small number of head 

observations and having data from only two streamflowgaging stations in the model active area introduces uncertainty, 

and again, different combinations of parameter values could 

provide the same calibrated model result.

Table 4. Water budget for a simulation of decrease in recharge during a hypothetical 10-year drought and effects on groundwater in 

storage and base flow to streams, Citizen Potawatomi Nation Tribal Jurisdictional Area, central Oklahoma.

[Volumes in acre-feet] 

Aquifer

Change in 

recharge during 

drought

Total 

groundwater 

in storage in 

2000 with no 

drought

Net change in 

groundwater 

storage during 

drought

Percent 

change in 

total aquifer 

storage

Net change 

in base flow 

during drought

Percent 

change in 

base flow

North Canadian alluvial aquifer -441,000 163,100 -14,100 -8.6 -386,500 -37

Central Oklahoma aquifer -696,000 147,910,000 -346,400 -0.2 -292,000 -28

Total -1,137,000 148,073,000 -351,500 -0.2 -678,500 -33
24 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

Representation of the groundwater-flow system in the 

transient model was difficult because of the lack of stresses in 

the aquifers such as pumping and a dearth of observations of 

the system response to stresses. This lack of data introduces 

uncertainty to estimates of aquifer storage parameters and was 

seen in the low sensitivity of model observations to storage 

parameters.

Analyses of available groundwater at EDZs in this study 

are approximate and are limited by the size of model cells 

and limited available subsurface hydrogeological data. The 

inset model used calibrated hydraulic values estimated by 

the CPN model and did not include local variation in the Iron 

Horse Industrial Park or the geothermal supply area. The 

GWM analysis could be improved by obtaining additional 

subsurface data such as direct measurements of the alluvium 

thickness, aquifer properties, and water levels from the local 

area. The continuous pumping rates presented in this report are 

approximate, and additional subsurface and aquifer test data 

could be acquired to verify available water.

The simulation of a hypothetical 10-year drought is 

limited to local changes in recharge and there are no data 

available for how the system responded to previous droughts. 

The analysis presented in this report did not include changes 

in runoff within the model area, flow entering the model in 

the North Canadian River, or pumping and diversions that 

may take place under drought conditions. Although recharge 

is estimated to be the greatest inflow to the aquifers, a more 

detailed analysis of drought-related flow changes in the North 

Canadian River and changes in water consumption by humans 

and plants may reduce the uncertainty of drought effects.

Date

1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008

Streamflow, in cubic feet per second

0

500

4,000

3,500

2,500

3,000

2,000

1,500

1,000

Percent change in streamflow

–60

–50

–40

–30

–20

–10

0

Drought begins

Drought ends

EXPLANATION

Streamflow during drought

Percent change

Normal streamflow

Percent change 50-month moving average

North Canadian River at Shawnee, Oklahoma

Figure 11. Observed North Canadian River streamflow at the Shawnee, Oklahoma, streamflow-gaging station and simulated flow with 

percent change during and after a hypothetical 10-year drought from 1990 to 2000.
Summary 25

Summary

A numerical groundwater-flow model was constructed 

using hydrogeologic information from the Central Oklahoma 

aquifer model that covered the Citizen Potawatomi Nation 

(CPN) Tribal Jurisdictional Area. The new model, termed 

the CPN model, estimated available groundwater resources 

and the effects of groundwater development and prolonged 

drought. The CPN model used a cell size of 820 feet (ft) on 

each side, which was substantially smaller than the cells in 

the Central Oklahoma aquifer model, which were 3,280 ft 

on each side. The CPN model also included alluvial deposits 

as hydrogeological units not differentiated in the Central 

Oklahoma aquifer model. Because of increased resolution and 

simulation of detailed flow on the eastern model margin, the 

CPN model also included drain cells to allow groundwater 

to discharge to the east. The CPN model simplified Central 

Oklahoma aquifer model layers to include only two layers 

in the Central Oklahoma aquifer—the upper and lower 

parts—and added an additional layer representing the part of 

the Wellington Formation that contains saline groundwater, 

referred to as the “saline unit.” The CPN model was calibrated 

using trial-and-error and automated parameter-estimation 

techniques to refine the calibration of the subregional Central 

Oklahoma aquifer model. 

To optimize the use of water resources at two economic 

development zones (EDZs) just south of the city of Shawnee, 

Okla., the U.S. Geological Survey Groundwater-Management 

Process (GWM) for MODFLOW, referred to as the GWM-VI 

version, was used with an inset model constructed from 

a part of the CPN model. Model cells in the inset model 

were reduced to 164 ft on all sides to analyze local-scale 

groundwater flow and interactions between pumping wells 

and the North Canadian River. Model layers and hydraulic 

properties were transferred directly from the calibrated CPN 

numerical groundwater-flow model. Groundwater withdrawal 

amounts were optimized using GWM at the geothermal supply 

area and the Iron Horse Industrial Park. The target combined 

pumping rate for all wells at the geothermal supply area was 

500 gallons per minute (gal/min). The Iron Horse Industrial 

Park required an estimate of the maximum groundwater 

development and the amount of North Canadian River water 

that may enter the wells. There were 5 managed wells used at 

the geothermal supply area and 16 managed wells used at the 

Iron Horse Industrial Park. 

Because the estimated mean thickness of the North 

Canadian River alluvial aquifer was approximately 43 ft, 

streamflow depletion was less of a concern than local well 

drawdown, and pumping was constrained by head drawdown 

controls at managed wells at both EDZ areas. Induced 

infiltration of streamflow was calculated to determine the 

potential for lower-quality river water entering the water 

pumped for the Iron Horse Industrial Park. At the geothermal 

supply area, maximum allowed drawdown was set at 16.4 ft, 

which was estimated to be the maximum drawdown that 

would not cause loss of well yield. The head constraints at the 

Iron Horse Industrial Park well field allowed the head to drop 

as much as one-half of the predevelopment saturated thickness 

at each well.

Five managed wells at the geothermal supply area 

produced a maximum combined simulated flow of 638 gallons 

per minute (gal/min) of continuous pumping, which exceeded 

the goal of 500 gal/min. At the Iron Horse Industrial Park, 

the 16 managed wells had a combined withdrawal rate of 

1,472 gal/min. Only 4.1 gal/min (0.3 percent) of the pumped 

water was estimated to come from induced infiltration of the 

North Canadian River streamflow

The effects of a hypothetical 10-year drought on 

groundwater in the CPN Tribal Jurisdictional Area and 

streamflow at the Shawnee, Okla., streamflow-gaging station 

were estimated using the CPN numerical groundwater-flow 

model. The drought was assumed to reduce precipitation and, 

in turn, recharge by 50 percent for the period 1990–2000. The 

model was run with recharge scaled down during the drought 

period and the resulting changes in total groundwater in 

storage and simulated streamflow at the Shawnee streamflowgaging station were estimated.

The decrease in recharge during the hypothetical 

drought caused groundwater in storage over the entire model 

to decrease by 351,500 acre-feet (acre-ft) (14,100 acre-ft 

or 8.6 percent in the North Canadian River alluvial aquifer 

and 346,400 acre-ft or 0.2 percent in the Central Oklahoma 

aquifer), or approximately 0.2 percent of all groundwater 

in storage at the end of the drought period. This volume of 

groundwater loss showed that the Central Oklahoma aquifer 

is a bedrock aquifer that has relatively low rates of recharge 

from the land surface. The budget for base flow to the North 

Canadian River estimated that the change in groundwater flow 

to the North Canadian River decreased during the 10-year 

drought by 386,500 acre-ft (37 percent). In all other parts of 

the CPN Tribal Jurisdictional Area, base flow decreased by 

292,000 acre-ft (28 percent). The total change in base flow 

to streams in the model was a decrease of 678,500 acre-ft 

(33 percent).

Streamflow in the North Canadian River was also 

analyzed during the hypothetical drought. The decrease 

of recharge to the North Canadian River alluvial aquifer 

was 441,000 acre-ft, which caused streamflow in the North 

Canadian River at the streamflow-gaging station at Shawnee 

to decrease by as much as 28 percent, and a 50-month moving 

average of simulated streamflow decreased by 16 percent. 

Simulated streamflow returned to nondrought streamflow rates 

approximately 7 years after the end of the drought.
26 Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN

References Cited

Adams, G.P., Runkle, D., Rea, A., and Becker, C.J., 1997, 

Digital datasets that describe aquifer characteristics of the 

alluvial and terrace deposits along the North Canadian 

River from Oklahoma City to Eufala Lake in East-Central 

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

96–448. [Also available at: http://pubs.usgs.gov/of/1996/

ofr96-448/].

Ahlfeld, D.P., Barlow, P.M., and Mulligan, A.E., 2005, 

GWM—A ground-water management process for the 

U.S. Geological Survey modular ground-water model 

(MODFLOW-2000): U.S. Geological Survey Open-File 

Report 2005-1072, 124 p. [Also available at http://pubs.

usgs.gov/of/2005/1072/pdf/ofr2005_1072.pdf.]

Anderson, M. P. and W. W. Woessner, 1992, Applied 

groundwater modeling: simulation of flow and advective 

transport; San Diego, London, Academic Press Inc., 381 p.

Banta, E.R., and Ahlfeld, D.P., 2013, GWM–VI—

Groundwater management with parallel processing for 

multiple MODFLOW versions: U.S. Geological Survey 

Techniques and Methods, book 6, chap. A48, 33 p. [Also 

available at http://dx.doi.org/10.3133/tm6a48/.]

Barlow, P.M., and Leake, S.A., 2012, Streamflow depletion 

by wells—Understanding and managing the effects of 

groundwater pumping on streamflow: U.S. Geological 

Survey Circular 1376, 84 p. [Also available at http://pubs.

usgs.gov/circ/1376/.] 

Becker, C.J., 2014, Stream-water and groundwater quality 

in and near the Citizen Potawatomi Nation Tribal 

Jurisdictional Area, Pottawatomie County, Oklahoma, 

2012–13: U.S. Geological Survey Scientific Investigations 

Report 2014–5178, 102 p., [Also available at http://dx.doi.

org/10.3133/sir20145178.]

Doherty, J., 2010, PEST, model independent parameter 

estimation user manual (5th ed.): Brisbane, Australia, 

Watermark Numerical Computing, accessed October 2, 

2008, at http://pesthomepage.org/.

Esri, 2015, An overview of the Spatial Analyst Toolbox, 

accessed January 28, 2015, at http://resources.arcgis.com/

en/help/main/10.2/index.html#/An_overview_of_the_

Spatial_Analyst_toolbox/009z00000003000000/.

Fetter, C.W., 1994, Applied hydrogeology (3d ed.), Upper 

Saddle River, New Jersey, Prentice-Hall, 691 p.

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

subregional water budgets using results from the U.S. 

Geological Survey modular three-dimensional ground-water 

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

90-392, 46 p. [Also available at http://water.usgs.gov/nrp/

gwsoftware/zonebud3/ofr90392.pdf.]

Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, 

M.G., 2000, MODFLOW-2000, the U.S. Geological 

Survey modular ground-water model -- User guide to 

modularization concepts and the Ground-Water Flow 

Process: U.S. Geological Survey Open-File Report 2000-92, 

121 p.

Havens, J.S., 1989, Geohydrology of the alluvial and terrace 

deposits of the North Canadian River from Oklahoma City 

to Eufaula Lake, central Oklahoma: U.S. Geological Survey 

Water-Resources Investigations Report 88-4234, 32 p., 12 

plates.

Heran, W.D., Green, G.N., and Stoeser, D.B., 2003, A digital 

geologic map database for the state of Oklahoma: U.S. 

Geological Survey Open-File Report 03-247, 10 p.

Konikow, L.F., Hornberger, G.Z., Halford, K.J., and Hanson, 

R.T., 2009, Revised multi-node well (MNW2) package for 

MODFLOW ground-water flow model: U.S. Geological 

Survey Techniques and Methods 6–A30, 67 p.

Mashburn, S.L., Ryter, D.W., Neel, C.R., Smith, S.J., and 

Magers, J.S., 2013, Hydrogeology and simulation of 

groundwater flow in the Central Oklahoma (GarberWellington) Aquifer, Oklahoma, 1987 to 2009, and 

simulation of available water in storage, 2010–2059: U.S. 

Geological Survey Scientific Investigations Report 2013–

5219, 92 p. [Also available at http://dx.doi.org/10.3133/

sir20135219.]

Merritt, M.L., and Konikow, L.F., 2000, Documentation of 

a computer program to simulate lake-aquifer interaction 

using the MODFLOW Ground-Water Flow Model and the 

MOC3D Solute-Transport Model: U.S. Geological Survey 

Water-Resources Investigations Report 00-4167, 146 p.

National Climatic Data Center, 2013, Locate weather 

observation station record: Asheville, N.C., National 

Climatic Data Center, data accessed June 15, 2013, at http://

www.ncdc.noaa.gov/oa/climate/stationlocator.html.

Niswonger, R.G. and Prudic, D.E., 2005, Documentation 

of the Streamflow-Routing (SFR2) Package to include 

unsaturated flow beneath streams--A modification to SFR1: 

U.S. Geological Survey Techniques and Methods, Book 6, 

Chap. A13, 47 p.

Niswonger, R.G., Panday, S., and Ibaraki, M., 2011, 

MODFLOW-NWT, A Newton formulation for 

MODFLOW-2005: U.S. Geological Survey Techniques and 

Methods 6–A37, 44 p.

Oklahoma Water Resources Board, 2012, 2012 Oklahoma 

Comprehensive Water Plan Executive Report: Oklahoma 

Water Resources Board, 159 p.

Oklahoma Water Resources Board, 2013, Water Well Record 

Search, accessed January 1, 2013, at http://www.owrb.

ok.gov/wd/search/search.php?type=wl.
References Cited 27

Shivers, M.J., and Andrews, W.J., 2013, Hydrologic drought 

of water year 2011 compared to four major drought periods 

of the 20th century in Oklahoma: U.S. Geological Survey 

Scientific Investigations Report 2013–5018, 52 p. [Also 

available at http://pubs.usgs.gov/sir/2013/5018/.]

Thornthwaite, C.W., and Mather, J.R., 1957, Instructions and 

tables for computing potential evapotranspiration and the 

water balance: Centerton, N.J., Laboratory of Climatology, 

Publications in Climatology, v. 10, no. 3, p. 185–311.

U.S. Army Corps of Engineers, 2013, Reservoir release, 

storage, evaporation, and precipitation data for lakes, 

accessed January, 2013 at: http://www.swt-wc.usace.army.

mil/stations.htm. 

U.S. Geological Survey, 2013a, National Elevation Dataset: 

accessed January, 2013 at http://ned.usgs.gov/index.html.

U.S. Geological Survey, 2013b, National Hydrography Dataset: 

accessed January, 2013 at http://nhd.usgs.gov/data.html.

U.S. Geological Survey, 2012, National Water Information 

System data available on the World Wide Web (USGS Water 

Data for the Nation), accessed December 2, 2013, at http://

waterdata.usgs.gov/nwis/.

Wahl, T.L., and Wahl, K.L., 2007, BFI—A computer program 

for determining an index to base flow, version 4.15: Bureau 

of Reclamation software release, accessed July 23, 2008, 

[Also available at http://www.usbr.gov/pmts/hydraulics_lab/

twahl/bfi/.]

Westenbroek, S.M., Kelson, V.A., Dripps, W.R., Hunt, R.J., 

and Bradbury,K.R., 2010, SWB—A modified ThornthwaiteMather Soil-Water-Balance code for estimating groundwater 

recharge: U.S. Geological Survey Techniques and Methods 

6-A31, 60 p.

Publishing support provided by

Lafayette Publishing Service Center
Ryter and others—Numerical Simulation of Groundwater Flow, Resource Optimization, and Potential Effects of Prolonged Drought for the CPN—SIR 2014–5167

ISSN 2328-031X (print)

ISSN 2328-0328 (online)

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

I SBN 978- 1- 4113-3932- 3