Prepared in cooperation with the

BUREAU OF RECLAMATION

Scientific Investigations Report 2005–5277

U.S. Department of the Interior

U.S. Geological Survey

Numerical Ground-Water Change Model of the C Aquifer

and Effects of Ground-Water Withdrawals on Stream

Depletion in Selected Reaches of Clear Creek,

Chevelon Creek, and the Little Colorado River, 

Northeastern Arizona

Numerical Ground-Water Change Model of 

the C Aquifer and Effects of Ground-Water 

Withdrawals on Stream Depletion in Selected 

Reaches of Clear Creek, Chevelon Creek, and 

the Little Colorado River, Northeastern Arizona

by S.A. Leake, J.P. Hoffmann and Jesse E. Dickinson

Scientific Investigations Report 2005–5277

Sacramento, California

2005

In cooperation with the Bureau of Reclamation
U.S. Department of the Interior

Gale A. Norton, Secretary

U.S. Geological Survey

P. Patrick Leahy, Acting Director

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

For Sale by U.S. Geological Survey, Information Services

Box 25286, Denver Federal Center

Denver, CO 80225-0286

For more information about the USGS and its products:

Telephone: 1-888-ASK-USGS

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

Any use of trade, product, or firm names in this publication is for descriptive purposes only and does not 

imply endosement by the U.S. Government.

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

owners to reproduce any copyrighted materials contained within this report.

Leake, S.A., Hoffmann, and Dickinson, J.E., 2005, Numerical ground-water change model of the 

C aquifer and effects of ground-water withdrawals on stream depletion in selected reaches of Clear Creek, 

Chevelon Creek, and the Little Colorado River, northeastern Arizona: U.S. Geological Survey Scientific 

Investigations Report 2005-5277, 29 p.
Contents

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

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

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

Previous Investigations ........................................................................................................................................... 3

Acknowledgments.................................................................................................................................................... 3

Physiographic and Hydrogeologic Setting................................................................................................................... 3

Boundary of the C Aquifer....................................................................................................................................... 4

Geologic Structure ................................................................................................................................................... 4

Geologic Units of the C Aquifer.............................................................................................................................. 4

Aquifer Properties .................................................................................................................................................... 6

Ground-Water Flow and Discharge....................................................................................................................... 6

Simulation Approach and Methods............................................................................................................................... 8

Summary of Numerical Ground-Water Change Model .............................................................................................. 9

Model Discretization and Boundary...................................................................................................................... 9

Aquifer Thickness................................................................................................................................................... 11

Surface-Water Features........................................................................................................................................ 11

Storage Properties ................................................................................................................................................. 13

Hydraulic Conductivity........................................................................................................................................... 13

Transmissivity.......................................................................................................................................................... 13

Withdrawal Scenarios ........................................................................................................................................... 14

Model Results.................................................................................................................................................................. 17

Monte Carlo Analysis of Model Results...................................................................................................................... 19

Limitations of the Model ................................................................................................................................................ 26

Summary and Conclusions............................................................................................................................................ 27

References....................................................................................................................................................................... 28

iii
Figures

Figure 1. Location of modeled area of the C aquifer, northeastern Arizona ................................................ 2

Figure 2. Generalized stratigraphic section of rock units in and adjacent to the C aquifer in the

study area, northeastern Arizona .................................................................................... 5

Figure 3. Estimated hydraulic conductivity values for the C aquifer, northeastern Arizona, compiled 

by Southwest Ground-water Consultants (2003), hydraulic conductivity zones, and contours of 

drawdown to conversion to unconfined conditions used in the change model ................................. 7

Figure 4. Lateral model boundaries, perennial reaches, confined areas, and initial transmissivity 

used in simulations of effects of withdrawals in the Leupp, Arizona, area .................................... 10

Figure 5. Thickness of the C aquifer, northeastern Arizona, as derived from a hydrogeologic 

framework model ................................................................................................... 12

Figure 6. Vertical section of a hypothetical partially confined aquifer. A, Actual geometry. 

B, Transformed geometry for use in a change model ............................................................ 14

Figure 7. Locations of existing and conceptual wells used in model simulation scenarios A and B 

for the C aquifer, northeastern Arizona ........................................................................... 15

Figure 8. Withdrawal scenarios for model simulations of changes in the C aquifer, 

northeastern Arizona ................................................................................................ 16

Figure 9. Total depletion for all stream reaches for model simulation scenarios A and B for the 

C aquifer, northeastern Arizona .................................................................................... 18

Figure 10. Depletion for lower Clear and Chevelon Creeks for model simulation scenarios A and B 

for the C aquifer, northeastern Arizona ........................................................................... 18

Figure 11. Computed maximum drawdown in the 101-year model simulation period for scenarios 

A and B for the C aquifer, northeastern Arizona .................................................................. 20

Figure 12. Locations of withdrawal centers, stream reaches, and wells with observations in a 

change model for the C aquifer, northeastern Arizona, that simulates historical conditions ............... 21

Figure 13. Observed and computed head change for three wells completed in the C aquifer, 

northeastern Arizona ................................................................................................ 24

Figure14. Results of Monte Carlo analyses of computed stream depletions in the study area, 

northeastern Arizona ................................................................................................ 26

iv

Tables

Table1. Locations of withdrawals used in model simulation scenarios A and B for the C aquifer, 

northeastern Arizona ................................................................................................ 17

Table 2. Locations and rates of withdrawals in the C aquifer used in the Western Navajo Hopi Study 

change model to simulate historical conditions .................................................................. 22

Table 3. Parameters used in the Monte Carlo analysis of computed stream depletions in the 

study area, northeastern Arizona ................................................................................. 25
Abbrevations

NAVD 88, Vertical coordinate information is referenced to the North American Vertical Datum of 1988. 

NAD 27, Horizontal coordinate information is referenced to the North American Datum of 1927.

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

Transmissivity, The standard unit for transmissivity is cubic foot per day per square foot times foot of 

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

(ft2/d), is used for convenience. 

acre-ft, acre foot

acre-ft/yr, acre foot per year

ft, foot

ft/day, foot per day

ft2/d, foot squared per day

ft3/s, cubic foot per second

ft-1, per foot

gal/day, gallon per day

in./yr, inch per year

m, meter

mi, mile

mi2, square mile

Acronyms

EDNA, Elevation Derivatives for National Applications

GIS, Geographic Information System

GWF, Ground-Water Flow (procedures of MODFLOW-2000)

KRB, riverbed conductance parameter

LCR, Little Colorado River

MODFLOW-2000, Modular Three-Dimensional Finite Difference Ground-Water Flow Model

OBS, Observation (procedures of MODFLOW-2000)

PES, Parameter Estimation (procedures of MODFLOW-2000)

RIVGRID, River Package data set program

SEN, Sensitivity (procedures of MODFLOW-2000)

USGS, U.S. Geological Survey

 v
This page intentionally left blank.
Abstract

The base flow in parts of Chevelon and Clear Creeks and 

of the Little Colorado River near Blue Spring in northeastern 

Arizona is sustained by discharge from the C aquifer, and in 

some reaches supports threatened and endangered fish species. C aquifer water is proposed as a replacement supply to 

relieve pumping from the N aquifer—the current source of water for a coal slurry pipeline used to transport coal mined from 

Black Mesa to Laughlin, Nevada. Locations of the proposed 

withdrawals are in the area of Leupp, Arizona, about 25 miles 

from a perennial reach of lower Clear Creek. A simulation tool 

was needed to determine possible effects of the proposed withdrawals from the C aquifer, particularly the effects of depletion of streamflow in Clear Creek, Chevelon Creek, and the 

Little Colorado River in the area of Blue Spring. A numerical 

ground-water change model was developed for this purpose. 

The model uses the U.S. Geological Survey finite-difference 

model code MODFLOW-2000 and data sets representing key 

features of the C aquifer to simulate change in the system that 

would result from withdrawing water at proposed locations. 

Aquifer thickness was estimated from a hydrogeologic framework model, and values of aquifer properties such as hydraulic 

conductivity and specific yield were estimated from aquifer-test 

data. Two scenarios with differing withdrawal rates were run 

for a 101-year period that included 51 years of withdrawals followed by 50 years of no withdrawals. About 6 percent of the 

ultimate volume of depletion occurs in the 101-year period for 

either scenario. The maximum streamflow depletion rate for all 

reaches in the scenario with the greatest withdrawal rates was 

computed to be about 0.6 cubic foot per second. The depletion 

rate was highest in lower Clear Creek, the reach that is closest 

to the well field. A model that simulates historical conditions 

was used to help select the most reasonable parameter sets for a 

Monte Carlo analysis of computed stream depletions.

Introduction

The C aquifer underlies much of northeastern Arizona, 

and parts of northwestern New Mexico (fig. 1), and is an important source of water for many users. It is named for the primary 

water-bearing rock unit of the aquifer, the Coconino Sandstone. 

The aquifer is the most productive aquifer in the Little Colorado 

River Basin. Industrial users, individual homeowners, agriculture, and the municipalities of Flagstaff, Holbrook, Leupp, Show 

Low, St. Johns, Joseph City, Snowflake, Taylor, and Winslow, 

Arizona, depend upon water drawn from wells completed in the 

aquifer. Because the aquifer lies beneath the Navajo and Hopi 

Indian Reservations, it represents an important present and future source of water to support these native peoples and their 

economies.

Water from the C aquifer is proposed as a replacement supply for mine-related withdrawals from the N aquifer—the current source of water for a coal slurry pipeline used to transport 

coal mined from Black Mesa to the Mojave Generating Station, 

an electrical power plant in Laughlin, Nevada. Water from the N 

aquifer generally is of good quality. Replacement of the supply 

by poorer quality from the C aquifer would reduce industrial 

use of the N aquifer water, a current supply for tribal communities. The proposed area of withdrawals for replacement water 

from the C aquifer is near Leupp, Arizona. A pipline would be 

required to transport the water to the coal mine on Black Mesa. 

Along with withdrawals for the mine, withdrawals also are proposed for Hopi and Navajo water supply.

The base flow of parts of some streams, including 

Chevelon and Clear Creeks, and of the Little Colorado River is 

sustained by discharge from the C aquifer, and in some cases, 

these streams support threatened and endangered fish species. 

The C aquifer discharges into the Verde River and Salt River 

Basins to the south, and to underlying formations and ultimately to the lower Little Colorado River to the northwest (fig. 1). 

Blue Spring, about 13 river miles (mi) upstream from the mouth 

of the Little Colorado River, is the largest of many springs that 

discharge in the lower reach of the river.

Numerical Ground-Water Change Model of the C Aquifer and 

Effects of Ground-Water Withdrawals on Stream Depletion in 

Selected Reaches of Clear Creek, Chevelon Creek, and the 

Little Colorado River, Northeastern Arizona

by S.A. Leake, J.P. Hoffmann, and Jesse E. Dickinson
Figure 1. Location of modeled area of the C aquifer, northeastern Arizona.

2

Colorado

River

Lit tle

Clear Creek

Page

Kayenta

Window Rock

Show Low

Snowflake

Taylor

Joseph City

Holbrook

Winslow

Leupp

Cameron

St. Johns

Blue Spring

BLACK

MESA

STRUCTURAL BASIN

MOGOLLON RIM

DEFIANCE UPLIFT

C o c o n i n o

A p a c h e

N a v a j o

G i l a

M c K i n l e y

C i b o l a

Greenlee

Salt Lake City

Mexican Hat

Flagstaff

Phoenix

Santa Fe

Albuquerque

A R I Z O N A

N E W M E X I C O

U T A H

Active

Model

Area

Model Boundary Origin

Model Boundary

33o

34o

36o

35o

108o 109o

112o 111o 110o

0

0 25 MILES

25 KILOMETERS

Flagstaff

Base from U.S. Geological Survey

digital data, 1:100,000, 1980

Lambert Conformal Conic projection

Standard parallels 2930' and 4530',

central meridian –11130'

U T A H

A R I Z O N A

C O L O R A D O

N E W M E X I C O

NAVAJO INDIAN

RESERVATION

RESERVATION

INDIAN

HOPI

113o

Color ado

River

EXPLANATION

PERENNIAL STREAM REACH

EXISTING WELL USED IN SCENARIOS

 TO COMPUTE STREAM DEPLETION

Chevelon

Creek

C a t r o n
Purpose and Scope

In 2004, the U.S. Geological Survey (USGS), in cooperation with the Bureau of Reclamation, was requested to participate on a Federal team to evaluate potential effects of future 

ground-water withdrawal scenarios from the C aquifer on base 

flow in the Chevelon and Clear Creek drainages, and groundwater discharge to the Little Colorado River in the area of Blue 

Spring. USGS participation on the Federal team involved assisting with test drilling and associated aquifer testing in the 

Leupp, Arizona, area (Hoffmann and others, in press) and with 

development of a numerical ground-water model of the C aquifer, discussed in this report.

The purpose of this report is to describe results of an 

analysis of pumping effects on base flow using a numerical 

ground-water model. The model was used to evaluate the effect 

of pumping on base flows of Chevelon and Clear Creeks and on 

ground-water discharge to the Little Colorado River near Blue 

Spring area for specified withdrawals and periods. Evaluations 

were carried out using two scenarios of withdrawals provided 

by a technical advisory group (Craig Sommers, Economist, 

ERO Resources Corporation, Denver, Colorado, written commun., 2005). Results are limited to calculated stream depletions 

and drawdown values. 

Previous Investigations

Parts of the C aquifer have been investigated to evaluate water availability for local water supply. Peter Mock Groundwater 

Consulting, Inc. (2003) constructed a ground-water flow model 

of a 2,400-mi2 rectangular part of the C aquifer encompassing 

Clear and Chevelon Creeks south and east of Leupp, Arizona. 

The ground-water flow model was used to estimate the changes 

in ground-water levels near the proposed well fields, at other 

existing and proposed water demand centers, and in base flow 

of Clear and Chevelon Creeks. Boundary conditions for this local-scale model were derived using a regional analytical model 

(Southwest Ground-water Consultants, Inc., 2003). The analytical model used the Theis nonequilibrium equation with the 

Jacob water-table correction, with a transmissivity of 25,000 

(gal/d)/ft (3,300 ft2/d) and a specific yield of 0.07.

Two recent studies have attempted to characterize the 

aquifer on a regional scale or over its full extent. In the area 

of Flagstaff, Arizona, the C aquifer was investigated by Bills 

and others (2000). This study included extensive structural 

mapping, geophysical surveying and mapping, determination 

of hydraulic properties, and chemical and isotopic analysis of 

ground water. A ground-water budget and a generalized characterization of the geohydrology of the entire C aquifer were 

developed by Hart and others (2002). Several other reports provide information on water resources in the study area including 

the C aquifer. Cooley and others (1969) describe the hydrogeology of the Navajo and Hopi Reservations and McGavock and 

others (1986), Mann and Nemecek (1983), and Mann (1976) 

describe water resources in southern parts of Coconino, Navajo, 

and Apache Counties, respectively.

Acknowledgments

Bill Greenslade and John Ward provided data sets that 

were helpful in providing information used to constrain aquifer thickness. Craig Sommers, Jason John, and Kevin Black 

were instrumental in providing information related to pumping scenarios. Don Bills provided guidance on the geology and 

hydrology of the study area, Alice Konieczki assembled much 

of the data needed from USGS sources, and Marilyn Flynn carried out Geographic Information System (GIS) analyses that 

were essential to constructing the model and presenting results. 

Claire Tiedeman provided guidance in strategies for calibration 

and Monte Carlo analyses. The modeling subcommittee members of the Federal team provided constructive input on model 

development.

Physiographic and 

Hydrogeologic Setting

The Little Colorado River Basin is in the southern part 

of the Colorado Plateau physiographic province (Fenneman, 

1946). The Little Colorado River drains an area of about 27,000 

mi2 of northeastern Arizona and northwestern New Mexico (fig. 

1). Elevations in the basin range from about 2,700 ft near the 

mouth of the Little Colorado River to 12,633 ft near Flagstaff. 

Average annual precipitation ranges from less than 6 in. in the 

lower areas to more than 30 in. in the highest elevations. Total 

annual average precipitation within the basin is estimated to 

be approximately 17,000,000 acre-ft. Potential evapotranspiration is estimated to range from less than 30 in./yr to more 

than 60 in./yr (Hart and others, 2002).

3
Boundary of the C Aquifer

The southern and southwestern boundary of the C aquifer extends south of the Little Colorado River Basin into the 

Verde River and Salt River Basins where rocks of the C aquifer 

are exposed along the Mogollon Rim (fig. 1). The northeastern 

boundary of the aquifer extends into New Mexico where the 

aquifer thins substantially (Hart and others, 2002). Few waterlevel and geologic data are available; therefore, the boundary 

is uncertain in these locations. The western boundary of the C 

aquifer is west of the Little Colorado River where the aquifer is 

unsaturated. Rocks of the C aquifer crop out in the southern and 

southwestern areas along the Mogollon Rim, in the northeastern area near the Defiance Uplift, and in the western area near 

the Little Colorado River (fig. 1).

Geologic Structure

The dominant structural feature within the C-aquifer area 

is the Black Mesa Basin, located in the north central part of the 

study area. The C aquifer generally dips toward the center of the 

basin from its outcrops along the margins of the study area. In 

the center of the Black Mesa Basin the aquifer is buried beneath 

nearly 5,000 ft of overlying geologic units (Southwest Groundwater Consultants, Inc., 2003; fig. 2) The Defiance Uplift is 

near the eastern edge of the C aquifer. Monoclines on either 

side of the uplift plunge to the west and east in the Black Mesa 

and San Juan Basins. Several regional-scale normal faults exist 

in the area and generally have two principal strike directions: 

north-northeast and north-northwest (Bills and others, 2000).

Geologic Units of the C Aquifer

The C aquifer comprises a sequence of sedimentary rock 

units between the top of the Kaibab Formation, which may 

underlie the Shinarump Member of the Chinle Formation or 

Moenkopi Formation, and the middle to upper part of the Supai 

Group (fig. 2; Cooley and others, 1969). The primary aquifer 

unit is the Coconino Sandstone and its lateral equivalents (the 

DeChelly Sandstone and Glorieta Sandstone in New Mexico); 

however, the overlying Kaibab and Toroweap Formations and 

the underlying Schnebly Hill Formation and Upper and Middle 

Supai Formations of the Supai Group locally can be significant 

water producing units. In the western part of the C aquifer near 

Flagstaff, units above the Supai Group are dry; the water table 

of the C aquifer there is in the upper part of the Supai Group 

(Bills and others, 2000).

The Kaibab Formation is primarily a limestone with a 

thickness of as much as 650 ft. It is absent in the eastern and 

northern part of Little Colorado River Basin. The Toroweap 

Formation underlies the Kaibab Formation and comprises carbonate sandstone, red beds, silty sandstone, siltstone, limestone, 

and thin layers of gypsum. Similar to the Kaibab Formation, 

the Toroweap is absent in the eastern and northern parts of the 

Little Colorado River Basin.

The Coconino Sandstone is a white to tan to light brown, 

crossbedded, eolian, fine-grained quartz sandstone. This unit 

generally is 300 to 1,100 ft thick (Bills and others, 2000). It can 

be extensively fractured where faulted. These fractured zones 

are likely to be areas of high permeability.

The Schnebly Hill Formation comprises a sequence of reddish-brown to reddish-orange very fine to silty sandstone, mudstone, limestone, and dolomite (Blakey, 1990). The formation 

ranges in thickness from a few feet to as much as 800 ft and 

generally underlies the Coconino Sandstone; however, the contact between the two units can be gradational and intertonguing. 

Where saturated, the Schnebly Hill Formation is in hydraulic 

connection with the Coconino Sandstone.

The Supai Group was divided into the Upper, Middle, and 

Lower Supai Formations by Blakey (1990). The Upper Supai 

Formation is a complex series of horizontally bedded reddish to 

brown sedimentary units that are mostly fine-grained sandstone, 

siltstone, and mudstone. The Middle Supai Formation is a grayish-orange, calcareous, very fine grained sandstone to siltstone. 

The Lower Supai Formation is a red to purple sandstone and 

siltstone and gray limestone and dolomite. The Supai Group is 

typically 600 to 2,400 ft thick and is exposed in the southern 

and southwestern part of the study area. Where saturated, the 

Upper Supai Formation is the most permeable unit in the group. 

The Lower Supai Formation typically forms a confining unit 

that separates the C aquifer from underlying saturated units.

4
Figure 2. Generalized stratigraphic section of rock units in and adjacent to the C aquifer in the study area, northeastern Arizona.

5

YOUNGER

VOLCANIC ROCKS

AND BASALTS

OLDER BASALTS

MOENKOPI FORMATION

KAIBAB FORMATION

COCONINO SANDSTONE

UPPER SUPAI FORMATION

SCHNEBLY HILL FORMATION

MIDDLE SUPAI FORMATION

LOWER SUPAI FORMATION

TEMPLE BUTTE (MARTIN) FORMATION

REDWALL LIMESTONE

MUAV LIMESTONE

SANDSTONE

GRANITE AND SCHIST

TAPEATS

BRIGHT ANGEL SHALE

TOROWEAP FORMATION

FORMATION

CHINLE

YOUNGER

BASALTS

EVAPORITE

DEPOSITS

UNCONSOLIDATED

SEDIMENTS

)5',)(:%7)5

bearing zone

Solution

channels

in caves

C aquifer

Supai

Group

Confining layer

)(:%//8%9

aquifer

The Moenkopi and Chinle Formations

form a confining layer where present

0–5,000

10’s–100’s

0–150

100–650

0–150

300–1,100

0–800

600–2,400

0–400

0–200

50–300

Modified from Bills and others (2000)

THICKNESS, IN FEET
Aquifer Properties

Aquifer tests results reported for these sedimentary units 

in the vicinity of Flagstaff indicate a significant variability in 

rock fracturing leading to a range of transmissivities, hydraulic 

conductivities, and storage coefficients (Bills and others, 2000; 

Hoffmann and others, in press). Transmissivity and hydraulic 

conductivity values for the Coconino Sandstone generally were 

higher than the values for the Kaibab Formation, the Schnebly 

Hill Formation, and the Supai Group. Regardless of lithology, 

wells with the highest transmissivity and hydraulic conductivity values were coincident with the greatest degree of fracturing 

and the lowest values were coincident with areas of limited or 

no visible fracturing. Transmissivity values range from 10 to 

4,700 ft2/d and hydraulic conductivity values range from 0.019 

to 6.88 ft/d for the Flagstaff area. In addition to the values published by Bills and others (2000), Southwest Ground-water 

Consultants, Inc. (2003) compiled aquifer-parameter estimates 

of the C aquifer from multiple sources. These sources indicate 

a large variation in transmissivity estimates (1 to 13,000 ft2/d). 

The large variance is attributed to factors such as degree of fracturing and differences in the penetration depths of wells used 

for aquifer tests. Many of the wells were designed to meet small 

demands for water and therefore were drilled to a depth where 

the required yield could be obtained. In many cases these wells 

penetrate only a part of the total aquifer thickness, resulting in 

an underestimate of transmissivity. In addition to transmissivity, Southwest Ground-water Consultants, Inc. (2003) compiled 

estimates of hydraulic conductivity values from tests in the 

C aquifer (fig. 3).

In confined parts of the aquifer, the storage 

coefficient [dimensionless] is related to the compressibility 

of the rock matrix and water. In terms used by ground-water 

hydrologists, the storage coefficient, S, can be computed as

 

 where b [length] is aquifer thickness, and Ss [length-1] is total 

specific storage, which is the sum of Ssw and Ssk, the components 

of specific storage due to compressibility of water and the aquifer skeleton, respectively. Values of Ss

 can be estimated from 

well-designed aquifer tests and also from properties of water 

and rocks. 

In unconfined parts of the area, the most important storage 

property is specific yield [dimensionless], which is a measure 

of the relative amount of water draining or filling pore spaces 

when the height of the water table changes. According to Bills 

and others (2000), the average of eight specific yield estimates 

derived from aquifer tests in the Flagstaff area was 0.077.

Several investigators (Mann, 1976; Daniel, 1981; Mann 

and Nemecek, 1983) have delineated areas where the C aquifer 

is confined and unconfined (fig. 3). For analyses of effects of 

withdrawal of water from the C aquifer, specific-yield values 

must be applied in areas where the aquifer is unconfined and 

storage coefficient values must be applied where the aquifer is 

confined. According to data compiled by Southwest Groundwater Consultants, Inc. (2003), the Coconino Sandstone responds to pumping as an unconfined or semiconfined aquifer, 

with storage properties of 0.07 for specific yield and 1 × 10-4 for 

the confined storage coefficient; however, there are few values 

of confined storage coefficient reported.

Ground-Water Flow and Discharge

Ground water flows from areas of recharge to areas of discharge. In the C aquifer, recharge occurs primarily in the southern part of the study area in the highlands north the Mogollon 

Rim, including the San Francisco volcanic field in the Flagstaff 

area, and on the Defiance Uplift in the eastern part of the area. 

From the recharge areas along the Mogollon Rim, some water flows to drainage basins south of the rim and some water 

flows northward, within the Little Colorado River drainage, 

ultimately toward the Blue Spring area. Ground water generally is unconfined in the southern and western parts of the 

aquifer (fig. 3) where the water table is below the elevation of 

the bottom of the Moenkopi Formation (fig. 2). The Moenkopi 

Formation comprises red to dark-red to reddish-brown siltstone, 

silty sandstone, fine- to very fine-grained sandstone, mudstone, 

and gypsum and is a confining layer in the north and east parts 

of the aquifer.

The primary discharge for the C aquifer in the Little 

Colorado River Basin is through a series of springs in the lower 

reaches of the Little Colorado River (Hart and others, 2002). 

The main spring is Blue Spring, which issues from solution 

channels in the Redwall-Muav aquifer in the lower reach of the 

Little Colorado River, about 13 mi upstream of the mouth of 

the river. The springs sometimes are collectively referred to as 

“Blue Springs.” Discharge from Blue Spring and other springs 

in the area is estimated to be approximately 164,000 acre-ft/

yr (Hart and others, 2002). All this discharge is thought to be 

downward leakage from the C aquifer because the RedwallMuav aquifer is deeply buried and does not crop out within the 

Little Colorado River Basin. An interpretation of the flow system is that water moves downward through the Supai Group 

through a series of normal faults (Cooley and others, 1969; Hart 

and others, 2002).

S b S , s = u

6
Figure 3. Estimated hydraulic conductivity values for the C aquifer, northeastern Arizona, compiled by Southwest Ground-water Consultants 

(2003), hydraulic conductivity zones, and contours of drawdown to conversion to unconfined conditions used in the change model.

7

C

olorado

River

109o

110o

111o

112o

0

0 25 MILES

25 KILOMETERS

Flagstaff

Base from U.S. Geological Survey

digital data, 1:100,000, 1980

Lambert Conformal Conic projection

Standard parallels 2930' and 4530',

central meridian –11130'

N E W M E X I C O

A R I Z O N A

EXPLANATION

1

1

1

3 4

4

2

2

2

2

36o

35o

34o

00

3 00

10 0

C AQUIFER CONFINED

C AQUIFER DRY

HYDRAULIC CONDUCTIVITY ZONE

 USED IN THE CHANGE MODEL

Less than 1

1–5

5–10

Greater than 10

ESTIMATED HYDRAULIC CONDUCTIVITY, IN FEET PER DAY:DRAWDOWN TO CONVERSION CONTOUR,

 IN FEET

Modified from Mann, 1976

Modified from Mann and Nemecek, 1983

Modified from Daniel, 1981

Estimated

20 0

0

30 0

100

1

0

0

100

200
The Redwall-Muav aquifer comprises the Redwall and 

Muav Limestones, the Martin Formation, and other Devonian 

rock units, where present. The Redwall Limestone is a lightgray to gray aphanitic to coarsely crystalline bedded limestone. 

Some beds within the formation are fractured, and the limestone 

typically contains solution cavities and caverns, some of which 

have collapsed. The Muav Limestone is a dolomitic limestone. 

Thickness of the Redwall-Muav aquifer ranges from about 200 

to more than 500 ft in the study area.

Ground water from the C aquifer also has historically 

discharged along or adjacent to various reaches of the Little 

Colorado River and its tributaries. The two tributary reaches 

of greatest interest for this study are the lower reaches of Clear 

and Chevelon Creeks. USGS streamflow-gaging stations were 

operated from 1906 to 1982 at Clear Creek near Winslow and 

from 1906 to 1972 at Chevelon Creek near Winslow. Few recent base-flow measurements exist along these creeks; however, measurements made in the summer of 2005 indicate baseflow values of about 5.6 and 2.7 ft3

/s (a combined total of about 

6,000 acre-ft/yr) in the lower reaches of Clear and Chevelon 

Creeks, respectively.

Simulation Approach and Methods

Ground-water models can be valuable tools for understanding the possible responses of aquifers to stresses such 

as withdrawal of water by wells. One approach to simulating 

the effects of withdrawal is to construct a numerical groundwater flow model with a program such as MODFLOW-2000 

(Harbaugh and others, 2000). In this approach, a ground-water 

flow equation with head as the dependent variable is solved. 

Flow into and out of the aquifer system is simulated, and the 

computed distribution of head can be compared with observations of head in wells. A simpler approach is to simulate 

changes in head in an aquifer. This approach is most common 

with analytical models, such as the Theis (1935) equation, in 

which change in head or drawdown is the dependent variable. 

If the aquifer-system response is linear, then the principle of 

superposition applies and simulation of changes can be referred 

to as “superposition modeling.” A problem with application of 

analytical superposition approaches (such as the Theis equation) for large-scale systems is that complex aquifer geometry 

and heterogeneous aquifer thickness and hydraulic properties 

cannot be simulated. For these situations, however, a superposition approach can be used with a numerical model program 

such as MODFLOW-2000. For details on the principle of superposition as used in ground-water hydrology, see Reilly and 

others (1987). For systems with nonlinearities such as 

transmissivity that changes as a function of head, the principle 

of superposition does not strictly apply.

Using a code such as MODFLOW-2000, a model can be 

set up to simulate changes in head in the aquifer and changes 

in flow to or from boundaries, while accounting for nonlinearities such as variable transmissivity. For the approach, referred 

to as “change modeling” here, the initial water surface (head) 

is set to a uniform elevation, referenced to an arbitrary datum. 

The head at all head-dependent boundary features is set to the 

same elevation and no recharge or discharge is specified so that 

no flow into or out of the aquifer or within the aquifer occurs in 

the initial state. Imposing a stress such as withdrawal of water 

by a well has the effect of propagating head changes into the 

surrounding model domain. In time, the effects can cause head 

changes in model cells containing surface-water boundaries. In 

response to that change in head, change in flow to or from the 

boundary is calculated.

In theory, any uniform initial head value could be used; 

however, numerical precision of calculations made with the 

model will be greatest when an initial elevation of 0.0 ft is used. 

This value also should be used for head in surface-water boundaries. To correctly simulate the aquifer thickness, the bottom 

of the aquifer must be set to an elevation equal to the initial 

head minus aquifer thickness. With an initial head of 0.0 ft, 

an estimated aquifer thickness of 520 ft would result in a bottom elevation of −520 ft. If at sometime during the simulation 

a head decline of 40 ft was computed in the change model for 

that location, the elevation of head in the model would be −40 ft 

and the saturated thickness would be updated automatically in 

MODFLOW to 480 ft and used in calculations of conductance 

between model cells.

The purpose for developing the change model described 

here was to estimate likely rates of capture of flow in selected 

streams as a result of proposed withdrawals of ground water in 

the area of Leupp, Arizona. Because of the short time frame 

provided for development of a simulation tool, a generalized simulation approach was used. Some key aspects of the 

approach are as follows:

1. The horizontal domain of the model is represented with 

a grid of cells that covers most of the C aquifer, using distant 

lateral boundaries so that computed effects in the area of interest will not be affected by uncertainties in the assignment of 

types and locations of boundaries.

2. The vertical domain of the model includes a single layer 

to represent the saturated part of the C aquifer, primarily the 

Coconino Sandstone. In the western part of the modeled area, 

the C aquifer is dry and the single layer represents the RedwallMuav aquifer that receives and transmits water flowing out of 

the C aquifer. No attempt is made to simulate the continuum 

of flow from the C aquifer to the underlying rock units with a 

configuration other than a single model layer.

8
3. The system is treated as a regional, nearly homogeneous 

porous media. No attempt is made to account for the effects of 

fractures and small-scale heterogeneities in aquifer properties.

4. Key perennial stream reaches are assumed to be 

hydraulically connected to the C aquifer.

5. Attempts to calibrate the change model are limited. 

Instead, reasonable average aquifer properties are used.

6. The model considers the effects of the proposed withdrawals near Leupp. Effects from other pumping in the regional 

aquifer are not considered. Scenarios for withdrawal provided 

by a technical team include 51 years of withdrawals at various 

rates and locations within a well field, followed by 50 years of 

no withdrawals.

Further details of the simulations are included in the 

following sections.

Summary of Numerical Ground-Water 

Change Model

The numerical code used to simulate changes in water 

levels in the C aquifer and ground water flow at the head-dependant boundaries is the Modular Three-Dimensional Finite 

Difference Ground-Water Flow Model (MODFLOW-2000) 

developed by the USGS (McDonald and Harbaugh, 1988; 

McDonald and Harbaugh, 1996; Harbaugh and others, 2000). 

Some of the Redwall-Muav aquifer in the western part of the 

domain also is included to form a single continuous transmissive system that extends to the major discharge area at and near 

Blue Spring.

Model Discretization and Boundary

The active model domain encompasses an area of about 

32,150 mi2, which includes much of the C aquifer beneath the 

Colorado Plateau and parts of the Redwall-Muav aquifer near 

Blue Spring (fig. 3). The model uses a uniform grid spacing of 

0.5 mi with 358 rows and 538 columns of finite difference cells. 

The model uses one layer; therefore, the model has a total of 

192,604 cells. There are 128,609 active cells, which represent 

areas where the C aquifer is saturated and part of the RedwallMuav aquifer near Blue Spring. The corner of the grid at row 

1, column 1 is at Universal Transverse Mercator (UTM) zone 

12 easting 518,610.25 meters, northing 4,155,142.5 meters. 

The grid is rotated 45 degrees clockwise so that columns 

increment in a southeasterly direction and rows increment in 

a southwesterly direction. This rotation angle corresponds 

approximately with direction of regional geologic structure.

The lateral boundary of the active area of the model shown 

on figure 4 is represented as a no-flow boundary. Individual 

segments of the boundary generally coincide with features as 

follows:

Boundary Segment Approximate location of feature

A to B Physical extent of the C aquifer along the Mogollon 

Rim

B to C Transition along the Mogollon Rim to a groundwater divide

C to D Ground-water divide

D to E Mesa Butte and Gray Mountain Faults

E to F Colorado River

F to G Line from the Colorado River to ground-water divide

near the Defiance Uplift

G to A Boundary where ground-water conditions are poorly

defined

The boundary segment E to F along the Colorado River is 

likely to have some hydraulic connection with the ground-water 

system yet is represented as a no-flow boundary. This representation of a boundary type, along with treatment of divides as 

immoveable no-flow boundaries, must be done with care in a 

model. In this case, the boundaries are at such great distances 

from the proposed withdrawals near Leupp that incorrect locations or boundary types will not affect computed results of capture in stream reaches of interest within the time frames used 

for the simulations.

The vertical domain of the model was simulated with 

one model layer. The system simulated includes the C aquifer and underlying aquifers in the western part of the area 

where the C aquifer is dry. For parts of the study area where 

the Coconino Sandstone has a substantial saturated thickness, 

most of the transmissivity of the C aquifer can be attributed 

to the Coconino Sandstone. Aquifer tests done in the C aquifer in the vicinity of Leupp, Arizona (Hoffmann and others, in 

press), indicate that the hydraulic conductivity of the Coconino 

Sandstone can be three orders of magnitude greater than that of 

the most permeable formation of the Supai Group, the Upper 

Supai Formation.

9
Figure 4. Lateral model boundaries, perennial reaches, confined areas, and initial transmissivity used in simulations of effects of withdrawals in 

the Leupp, Arizona, area.

10

Clear Creek

109o

110o

111o

112o

Flagstaff

Base from U.S. Geological Survey

digital data, 1:100,000, 1980

Lambert Conformal Conic projection

Standard parallels 2930' and 4530',

central meridian –11130'

N E W M E X I C O

A R I Z O N A

Chevelon

Creek

36o

35o

34o

Little

Colorado

River

0

0 25 MILES

25 KILOMETERS

EXPLANATION

Leupp

C AQUIFER TRANSMISSIVITY, IN

THOUSANDS OF FEET PER DAY

DRY 1 2 3 4 7 10

C AQUIFER CONFINED

NUMBERED PERENNIAL REACH

POINT ALONG MODEL BOUNDARY USED TO

 REFERENCE SEGMENTS FOR DISCUSSION

 OF BOUNDARY PLACEMENT

EXISTING WELL USED IN WITHDRAWAL

 SCENARIO

3–Upper

 Clear

 Creek 5–Upper

 Chevelon

 Creek

4–Lower

 Chevelon

 Creek

A

A

B

C

D

G

E

F

2–Lower Clear Creek

1–Little Colorado River

 below Blue Spring
Aquifer Thickness

A geologic framework model was used to estimate aquifer thickness throughout the area. The framework model represents three-dimensional relationships between the C aquifer, 

the surrounding confining units, and the water table. In the 

confined area, aquifer thickness is defined as the difference in 

elevation between the top of the Kaibab Formation or Coconino 

Sandstone and the top of the Supai Group. In the unconfined 

area, aquifer thickness is defined as the difference in elevation 

between the water table and the top of the Supai Group. The 

framework model was constructed using Dynamic Graphics, 

Inc. EarthVision® software. The grid consists of 358 rows and 

538 columns and is rotated at 45 degrees from the horizontal, 

which is the same discretization and orientation of the numerical ground-water model. Thickness is calculated for coordinates 

that coincide with the cell centers of the numerical model.

The framework model was constrained on the basis of 

well logs from 795 wells throughout the study area and from 

elevations of stratigraphic units exposed on the Mogollon Rim 

and on the south rim of the Grand Canyon (fig. 2). Of these 

wells, 559 have well logs that identify elevations for the top 

of the C aquifer; 386 well logs have elevations for the top of 

the Supai Group. Geologic sections (Southwest Ground-water 

Consultants, Inc., 2003) and structural contours of the top of the 

Coconino Sandstone (Mann, 1976; Mann and Nemecek, 1983) 

also were used. Because cross-section and contour data were 

mostly interpolated or estimated on the basis of prior geologic 

judgment, these data were given less weight in the calculations 

of the top and bottom of the C aquifer. The piezometric surface 

was digitized from maps in Hart and others (2002). Geologic 

faults were not included in the construction of the framework 

model because of a lack of data for these features. This simplification results in a continuous layer in which steep variations 

in thickness can occur where faults might exist. The framework 

model calculated some areas with thicknesses of less than 300 

ft and pinch outs of the C aquifer. The thin areas and pinch 

outs occurred because the framework model fits smooth planes 

using minimum-tension gridding to represent the top and bottom of the unit—calculations for areas having sparse data could 

result in thinning because one or both of the surfaces was either 

too low or too high in elevation. A thickness of 300 ft was assumed for these areas to assure horizontal continuity of the C 

aquifer. The framework model was assumed to be correct in 

areas having large thicknesses. The thickness of the C aquifer 

is generally unknown in the New Mexico part of the study area, 

but the productive thickness of the C aquifer is reported to be 

about 300 ft (Cooley and others, 1969).

The thickness of the C aquifer in the framework model 

varies between 300 and 2,089 ft (fig. 5); the average thickness 

of the aquifer is 390 ft. The C aquifer is dry in areas where 

faulting or folding has caused the Coconino Sandstone to be 

above the water table. Unless dry, for the purpose of this investigation, the minimum thickness applied to the aquifer is 300 ft. 

To the west, where the C aquifer is dry owing to drainage into 

the Redwall-Muav aquifer, aquifer thickness is generally held 

constant at about 300 ft (fig. 5). The thickness in the extreme 

western part of the model is typically greater than 700 ft, whereas much of the eastern part generally has a thickness less than 

500 ft. The two thickest areas include 1,600 ft or more of saturated sediments in the western and northern parts of the study 

area. A thickness of 300 ft was used for the middle of the basin, 

northwest of Holbrook.

Surface-Water Features

For simulations of effects of possible streamflow capture 

by proposed ground-water withdrawals near Leupp, the following five stream reaches were represented in the model: (1) 

upper Clear Creek, (2) lower Clear Creek, (3) upper Chevelon 

Creek, (4) lower Chevelon Creek, and (5) the Little Colorado 

River below Blue Spring (fig. 4). Upper and lower limits of 

perennial reaches 1 and 3 were taken from Brown and others 

(1981). Upper and lower limits of perennial reaches 2 and 4 

were taken from USGS field investigations made during the 

summer of 2005 (D.J. Bills, U.S. Geological Survey, written commun., 2005). Reach 5 was taken to be the trace of the 

Little Colorado River from Blue Spring to the mouth, where 

much of the discharge of water moving through the C aquifer 

ultimately occurs. Coordinates for traces of all segments were 

derived from the USGS Elevation Derivatives for National 

Applications (EDNA) database (U.S. Geological Survey, 2005). 

All five reaches were represented in the model using the River 

Package of MODFLOW-2000 and a river head elevation of 0.0 

ft, the same elevation as starting head in the aquifer. With this 

approach, flow between a surface-water feature and an aquifer 

is proportional to the head difference between the surface water 

and ground water in the underlying model cell. The coefficient 

of proportionality is the riverbed conductance. In the case of 

the change model constructed for this study, the river package computes change in flow in response to change in head in 

model cells underlying stream segments represented with the 

River Package. The change model does not address whether the 

change is a decrease in existing ground-water discharge to a 

stream, or an increase in existing leakage from a stream to the 

aquifer. In either case, however, a computed change in response 

to ground-water withdrawal would result in a net decrease in 

stream base flow. For this study, the change model was not configured to limit the amount of capture that would be available 

from a given stream reach.

11
Figure 5. Thickness of the C aquifer, northeastern Arizona, as derived from a hydrogeologic framework model.

12

C

olorado

River

109o

110o

111o

112o

Base from U.S. Geological Survey

digital data, 1:100,000, 1980

Lambert Conformal Conic projection

Standard parallels 2930' and 4530',

central meridian –11130'

N E W M E X I C O

A R I Z O N A

36o

35o

34o

0

0 25 MILES

25 KILOMETERS

EXPLANATION

Page

Window Rock

St. Johns

Kayenta

800

600

4

4 0 0

400

60 0

80 0

1,0 0 0

,01 00

800

600

400

800

60

400

0

600

800

800

1,000

2 1, 00

00

Leupp

Holbrook

Cameron

C AQUIFER THICKNESS, IN FEET

300 400 600 800 1,000 1,200 1,400 1,600 1,800

TRACE OF SECTION FOR WHICH ELEVATION

 DATA ARE AVAILABLE

WELL HAVING ELEVATION DATA:

 Top of C aquifer

 Bottom of C aquifer LINE OF EQUAL THICKNESS OF

 THE C AQUIFER—Approximately

 located. Interval 100 and 200 feet

600

Flagstaff

Winslow

Show Low
The riverbed conductance parameter required for the River 

Package data set is not known and cannot be estimated by calibration of the change model because observations of changes 

in ground-water flow to or from the streams are not available. 

For simulations of a 51-year withdrawal scenario, reducing the 

value of riverbed conductance has the effect of reducing the 

maximum rate of computed captured streamflow and extending the period in time over which capture occurs. Increasing 

riverbed conductance has the converse effect. A point exists, 

however, for which further increases in riverbed conductance 

will result in no further increases in maximum capture rate. The 

approach taken here is conservative with respect to the potential for capturing streamflow. The conservative approach was 

accomplished by using relatively large riverbed conductance 

values so that simulations would yield possible maximum rates 

in maximum capture over time.

The River Package data set for all five reaches was constructed with program RIVGRID (Leake and Claar, 1999). This 

program intersects a stream trace with individual model grids 

and computes the length of the stream in each model cell traversed. For all stream reaches in the C aquifer change model, 

initial riverbed conductance was computed by using a vertical 

hydraulic conductivity of 1 ft/d and a riverbed thickness of 20 ft. 

The assumed widths used for calculating riverbed-conductance 

values were 10 ft for reaches 1 through 4 and 20 ft for reach 5. 

The elevation for converting to steady river leakage was set to 

−30 ft so that capture would not be limited by drawdown values that might be computed under streams in this model. With 

these values, RIVGRID computed initial conductance values 

that averaged about 1,000 ft2/d for reaches 1 through 4 and 

about 2,000 ft2/d for reach 5. For final simulations, a riverbed 

conductance parameter (KRB) was specified in the Sensitivity 

Process of MODFLOW-2000 as a multiplier for all riverbedconductance values in all five reaches. Final simulations used 

a value of KRB of 100, resulting in average conductance values of about 1×105

 ft2/d for reaches 1 through 4 and about 

2×105 ft2/d for reach 5. The higher conductance values lead to 

computation of a maximum capture rate.

Storage Properties

The change model was constructed using the LayerProperty Flow Package of MODFLOW-2000 with the option of 

“convertible” layer type. With this option, changes in head can 

cause parts of the model to switch between confined and unconfined conditions. This option requires input of specific-yield 

and specific-storage values. MODFLOW-2000 uses specific 

yield as the storage property if head in a cell is lower than the 

top of an aquifer, and the product of specific storage and aquifer thickness if head is higher than the top of the aquifer. For 

the change model in which only drawdown (downward head 

changes) will be simulated, the top elevation in the unconfined 

area can be arbitrarily specified above the initial head of 0.0 ft. 

In an area that is initially confined, the top elevation must be the 

elevation at which the system would become unconfined when 

drawdown occurs. The strategy for computing top and 

bottom elevation arrays is outlined in figure 6. Mann (1976) and 

Mann and Nemecek (1983) provide contours of the pressure 

head above the top of confined parts of the C aquifer in southern Navajo and Apache Counties, Arizona, respectively. These 

pressure heads were used to compute the “drawdown to conversion” for southern Navajo and Apache Counties. The maximum 

drawdown to conversion, Smax on figure 6, was set at 300 ft. The 

shape of the surface of drawdown to conversion was extended 

outside of southern Apache and Navajo Counties to approximate the surface within the entire confined area (fig. 3).

For final simulations, a value of 0.06 was used for specific yield. This was selected on the basis of aquifer tests in 

the Leupp area (Hoffmann and others, in press.) A value of 

2×10-6 ft-1 was used for specific storage on the basis of 

properties of water and estimated skeletal specific storage for 

sandstone.

Hydraulic Conductivity

Hydraulic conductivity in the C aquifer is highly variable. 

Values from individual wells are subject to local influences 

such as well construction and small-scale geologic features that 

cannot be represented in a regional flow model. On the basis 

of hydraulic conductivity estimates compiled by Southwest 

Ground-water Consultants, Inc. (2003), the four zones for hydraulic conductivity shown in figure 3 were initially selected. 

Preliminary sensitivity tests using the change model showed 

that predicted stream depletion was sensitive to the value of hydraulic conductivity in zone 1 but was much less sensitive to 

the values in other zones. This result is to be expected because 

the proposed withdrawals near Leupp as well as all but one of 

the stream reaches simulated are within zone 1. The value of 

hydraulic conductivity for zone 1 was set at 6 ft/d on the basis 

of specific capacity data compiled by Southwest Ground-water Consultants, Inc. (2003) and hydraulic conductivity values 

determined for the Coconino Sandstone and the Schnebly Hill 

Formation by Hoffmann and others (in press). The value of hydraulic conductivity for zones 2, 3, and 4 was set at 5 ft/d, an estimate of the average value of the property for the C aquifer. For 

the Monte Carlo analysis, discussed in a later section, values of 

hydraulic conductivity were allowed to vary in all four zones.

Transmissivity

MODFLOW-2000 computes transmissivity as the product 

of hydraulic conductivity and saturated thickness for each grid 

cell in the active flow domain. For purposes of displaying this 

fundamental aquifer property, the distribution of initial transmissivity values was computed from initial thickness and hydraulic 

conductivity for each model cell location and mapped for the 

model area (fig. 4). This distribution is largely a reflection of 

the variations in the saturated thickness. The minimum value is 

1,500 ft2/d, the product of the minimum thickness (300 ft) and 

the minimum hydraulic conductivity (5 ft/d).

13
Withdrawal Scenarios

Two withdrawal scenarios were simulated with the numerical change model. Both scenarios include withdrawals for 

the 51-year period 2010–2060 and no withdrawals for the 50-

year period 2061–2110. For the period 2010–2060, multiple 

MODFLOW-2000 stress periods with 1-year time steps were 

used to represent variations in rates of withdrawals over time. 

A single MODFLOW-2000 stress period with fifty 1-year time 

steps was used to represent the recovery from 2061 to 2110. 

The locations of four existing wells and as many as 14 “conceptual” wells were used for both scenarios (table 1, fig. 7). 

Total withdrawals were subdivided into two use categories: 

(1) project withdrawals that include mine slurry and mine domestic use, and (2) withdrawals for Hopi and Navajo Tribal use. 

For simulations of scenarios A and B, project withdrawals for 

mine slurry and mine domestic use were assumed to take place 

at three test wells near Leupp, Arizona, drilled in 2005, and 

withdrawals for tribal use were distributed among the existing 

Sunshine well and the conceptual well locations (table 1, fig. 7). 

Scenario A simulates a maximum withdrawal rate of about 

6,500 acre-ft/yr between 2010 and 2060 (fig. 8A). Scenario B 

simulates as much as 11,500 acre-ft/yr (fig. 8B). Scenario A represents the minimum projected demands for the area, whereas 

scenario B represents the maximum projected demands for the 

area. For scenario A, 6,000 acre-ft/yr is withdrawn for industrial use between 2010 and 2026; the remainder (500 acre-ft/yr) 

is withdrawn for tribal use. Between 2026 and 2028, industrial 

use is reduced to a total of 505 acre-ft/yr while the tribal use 

is ramped up toward 6,500 acre-ft/yr (fig. 8A). After 2029, all 

pumpage is used for tribal purposes. Industrial pumpage is extracted from three wells (PW-1, PW-2, and PW-3) that are simulated in three model cells (table 1). Tribal pumpage is extracted 

from seven wells (NPW-1 through NPW-7) that are simulated 

in seven model cells (table 1).

Figure 6. Vertical section of a hypothetical partially confined aquifer. A, Actual geometry. B, Transformed geometry for use in a change model.

14

Top of aquifer

Bottom of aquifer

Smax

Smax

Bottom of aquifer

Initial potentiometric surface

Smax

B. Transformed geometry

NOTE: Initial aquifer thickness at all points, represented here with values at selected

 points b1, b2,....b10, is preserved in the transformed system. Values of drawdown

 to conversion to water table conditions, represented here with values at selected

 points S1, S2,....Smax, also is preserved in the transformed system up to the

 value of Smax. Use of the transformed geometry in a change model will be valid

 as long as computed drawdown does not exceed Smax.

A. Actual geometry

Smax

Smax S2

S1

S2 S1

b6

b5

b6

b5

b4

b4

b3

b3

b2

b1

b2

b1

Initial water table

b10

b10

Smax

b9

Smax

b8

b7 b8

b7

Smax

b9

Initial water table Initial potentiometric surface

Top of aquifer

Confining unit

Aquifer

Aquifer
Figure 7. Locations of existing and conceptual wells used in model simulation scenarios A and B for the C aquifer, northeastern Arizona.

15

EXPLANATION

C

C

l

Che e onl v

ear

reek

Creek

CONCEPTUAL WELL

EXISTING WELL AND NAME

0

0 5 MILES

5 KILOMETERS

Base from U.S. Geological Survey

digital data, 1:100,000, 1980

Lambert Conformal Conic projection

Standard parallels 2930' and 4530',

central meridian –11130'

111o 45’ 110o30’

Leupp

L

C

River

o orado

l

ittle

Winslow

34o45’

111o15’

35o30’

35o

15’

PW-1

PW-2

PW-3

PW-3

Sunshine

PERENNIAL STREAM

 REACHES
Figure 8. Withdrawal scenarios for model simulations of changes in 

the C aquifer, northeastern Arizona.

16 WITHDRAWAL RATE, IN ACRE-FEET PER YEAR

12,000

10,000

8,000

6,000

4,000

2,000

0

12,000

10,000

8,000

6,000

4,000

2,000

0

2000 2010 2020 2030 2040 2050 2060 2070

B. Scenario B

A. Scenario A

Project

withdrawals Tribal withdrawals

Project

withdrawals Tribal withdrawals
For scenario B, industrial withdrawal rates are the same as 

those in scenario A. The increased withdrawal rates in scenario 

B, relative to scenario A, result from increased withdrawals for 

tribal purposes. For scenario B, tribal withdrawals are assumed 

to occur in as many as 15 wells represented in 15 model cells 

(table 1).

The total volume withdrawn in the 51-year period is about 

331,000 acre-ft for scenario A and about 464,000 acre-ft for 

scenario B. Because withdrawal locations are far from perennial stream reaches, depletion can occur long after withdrawals 

cease. This lack of immediate response occurs because recovery from shutting off withdrawals takes time to reach distant 

parts of the outward propagating cone of depression from the 

51-year period of withdrawals. To simulate the continued depletion, each scenario therefore includes a 50-year period of no 

withdrawals.

Model Results

The model can make calculations of stream depletion and 

drawdown and display numbers with seven to eight significant 

digits of precision. For this discussion and the following section 

on the Monte Carlo analysis, most values of computed depletion of streamflow are less than 1 ft3

/s. Numbers for computed 

depletion are presented to the nearest 0.01 ft3/s; however, this 

level of precision does not imply predictive accuracy to that 

level. Also, changes in base flow in the range of 0.01–0.1 ft3/s 

would be hard to detect in small streams without special monitoring equipment such as precision flumes. Furthermore, natural variations in base flow likely would be this amount or larger. 

This level, however, allows for comparison of relative amounts 

of computed depletion among different reaches, some of which 

have computed depletions in the range of hundredths of a cubic foot per second. Graphs constructed from model-calculated 

depletion rates use the full precision of numbers recorded in 

model runs. Drawdown values are presented here for 1- and 

10-ft contour intervals.

Table 1. Locations of withdrawals used in model simulation scenarios A and B for the C aquifer, 

northeastern Arizona.

Well name UTM easting 

(meters)

UTM northing 

(meters)

Model 

row

Model 

column

Active in 

scenario A?

Active in 

scenario B?

Existing wells

PW-1 490360 3892073 257 207 yes yes

PW-2 496465 3895430 248 209 yes yes

PW-3 505487 3891199 244 221 yes yes

Sunshine 497291 3885968 256 218 no yes

Conceptual wells

NPW-1 492833 3891689 255 209 yes yes

NPW-3 495797 3891654 252 212 yes yes

NPW-6 499153 3891581 249 215 yes yes

NPW-7 502460 3891581 246 218 yes yes

NPW-2 493310 3893513 253 208 yes yes

NPW-5 499732 3894482 246 213 yes yes

NPW-4 497778 3893592 249 212 yes yes

NPW-8 494834 3895033 250 208 no yes

NPW-9 495593 3893475 251 210 no yes

NPW-10 497032 3894553 248 211 no yes

NPW-11 497471 3891556 251 214 no yes

NPW-12 499629 3893075 247 214 no yes

HPW-1 491677 3888360 259 211 no yes

HPW-2 493875 3887001 258 214 no yes

17
The change model simulates stream depletions resulting 

from withdrawal scenarios A and B. For this analysis, the only 

withdrawals simulated were at rates and locations described in 

the previous section. For the 101-year simulation period, the 

original lengths of the five perennial reaches were assumed 

to remain unchanged by the simulated withdrawals or by any 

other stresses. This has the effect of representing a conservative or “worst-case” scenario of total stream depletion because 

actual stream depletion from the proposed withdrawals would 

be reduced if perennial reaches were shortened or disappeared 

during the simulation period. The model was configured to accumulate computed depletions separately for each of the five 

reaches at 5-year intervals from the start of the simulations and 

at additional times of 1, 51, and 101 years from the start of the 

simulations.

The total stream depletion from all five reaches at the end 

of the 101-year simulation period is about 20,000 acre-ft for 

scenario A and about 26,000 acre-ft for scenario B, about 6 percent of the total volume withdrawn for each scenario. About 

94 percent of the depletion therefore would occur after the 101-

year period, assuming that the major discharge mechanisms are 

represented in the model and that the withdrawals do not induce 

additional recharge to the aquifer.

The time distribution of the rate of total depletion from all 

five stream reaches simulated indicates a gradual increase in 

stream capture to a maximum rate of depletion that occurs near 

the end of the 101-year simulation period (fig. 9). Computed 

depletion for scenario A was 0.31 ft3/s at the end of the 51-year 

withdrawal period and increased to a maximum of about 0.45 

ft3/s around 95 years after the start of the simulation. Computed 

depletion for scenario B was 0.37 ft3/s at the end of the 51-year 

withdrawal period and increased to a maximum of 0.63 ft3/s 

around 95 years after the start of the simulation.

Of the total depletion shown in figure 9 for scenarios A 

and B during the 101-year period, the relative amount occurring in each stream reach varies only slightly. The combined 

rates of depletion from upper Clear Creek, upper Chevelon 

Creek, and the Little Colorado River below Blue Spring increase from zero to about 3 percent of the total rate of depletion 

over the 101-year simulation period. The maximum computed 

depletion rate for any of these three reaches for scenario B is 

0.015 ft3/s. Because of the small values of computed depletion 

in these three reaches, only depletion in lower Clear Creek and 

lower Chevelon Creek will be discussed in detail in this report.

Computed rates of depletion for lower Clear Creek (fig. 10

are about 83 to 84 percent of the total depletion rates for both 

scenarios. Most of the depletion occurs in lower Clear Creek 

because the creek is the closest stream reach to the withdrawal 

locations. For scenario A, computed depletion in lower Clear 

Creek was 0.26 ft3/s at the end of the 51-year withdrawal period and increased to a maximum of about 0.39 ft3

/s around 85 

years after the start of the simulation. For scenario B, computed 

Figure 9. Total depletion for all stream reaches for model simulation 

scenarios A and B for the C aquifer, northeastern Arizona.

Figure 10. Depletion for lower Clear and Chevelon Creeks for model 

simulation scenarios A and B for the C aquifer, northeastern Arizona.

18

0.4

TOTAL STREAM DEPLETION, IN CUBIC FEET PER SECOND

0.3

0.2

0.1

0.7

0.6

0.5

0

TIME, IN YEARS SINCE START OF WITHDRAWALS

0 10 20 30 40 50 60 70 80 90 100 110

Scenario

B

Scenario A

0.4

TOTAL STREAM DEPLETION, IN CUBIC FEET PER SECOND

0.3

0.2

0.1

0.6

0.5

0

TIME, IN YEARS SINCE START OF WITHDRAWALS

Lower Chevelon Creek

Lower Clear Creek

0 10 20 30 40 50 60 70 80 90 100 110

Scenario

B

Scenario A

Scenario A

Scenario B
depletion was 0.31 ft3

/s at the end of the 51-year withdrawal 

period and increased to a maximum of about 0.53 ft3/s around 

90 to 95 years after the start of the simulation.

Computed rates of depletion for lower Chevelon Creek 

(fig. 10) are 14–17 percent of the total rate. Computed depletion 

for scenario A in lower Chevelon Creek reaches 0.05 ft3/s at the 

end of the 51-year withdrawal period and continues to increase 

to a maximum of about 0.06 ft3/s near the end of the simulation. 

Computed depletion for scenario B reaches 0.06 ft3/s at the end 

of the 51-year withdrawal period and continues to increase to a 

maximum of about 0.09 ft3/s near the end of the simulation.

Another use of the change model is to compute drawdown 

(negative change in head) in the aquifer. From the model results, 

time series of drawdown can be displayed for any model cell or 

the distribution of drawdown in the active model area can be 

displayed for different times in the 101-year simulation period. 

The model was configured to save head changes at 5-year intervals starting with 1 year and ending at 101 years. The model includes increases and decreases in withdrawals, and head 

change therefore could be negative (decreasing head) or positive (increasing head). For purposes of displaying drawdown 

in this report, all head changes were analyzed to determine the 

maximum drawdown at each active model cell for the 101-year 

period. Contours of 1 ft and 10 ft of maximum drawdown (fig. 

11) reflect some of the large-scale features represented in the 

model. The contours bulge noticeably into the confined area, 

reflecting the lower storage factor in that area. The 1-ft contour 

is noticeably affected by the presence of lower Clear Creek, 

a source of water that limits the propagation of the computed 

cone of depression in this area.

Monte Carlo Analysis of Model Results

An indication of possible ranges of computed future stream 

depletions can be determined with a Monte Carlo analysis of 

scenarios A and B. The general procedure for these analyses is 

to (1) generate N sets of model parameters by perturbing each 

parameter with random numbers, (2) run the model with each 

parameter set, resulting in N sets of predictions, and (3) compute the statistical properties of the predictions. Ideally, the N 

sets of parameters should be generated by using standard deviations of parameters and the covariances between parameters 

determined by using an automatic calibration procedure.

Such an effort was attempted for the C aquifer change model using the Ground-Water Flow (GWF), Parameter Estimation 

(PES), Sensitivity (SEN), and Observation (OBS) Processes of 

MODFLOW-2000 (Harbaugh and others, 2000; Hill and others, 2000). For this purpose, a model was constructed to simulate changes in ground-water levels for the period 1961–2001, 

using five stress periods covering time periods 1961–1970, 

1971–1980, 1981–1990, 1991–2000, and 2001, respectively. 

Withdrawal rates and locations for the first four periods for 

various withdrawal centers (fig. 12, table 2) were taken from the 

Western Hopi Navajo analytical model of effects of withdrawals in the C aquifer (Bill Greenslade, Southwest Ground-water 

Consultants, Inc., written commun., 2005.) The final stress period, year 2001, was added because of the many water levels 

that were measured during that year and that were potentially 

useful for calibration of the model. Withdrawal locations and 

rates were assumed to be the same as those for the previous period, 1991–2000. For the Western Hopi Navajo withdrawal data 

set used here, exact locations of wells withdrawing water were 

not known. The Western Hopi Navajo Study assumed locations 

for two or more wells in each center that were generally within 

a 2- to 3-mi radius of the locations of the withdrawal centers 

shown on figure 12.

The OBS Process in the MODFLOW-2000 provides a 

means to enter observations for use in constructing an objective 

function. The PES Process uses the objective function, the sum 

of the squares of weighted differences between computed values and observations (weighted residuals), to select an optimum 

set of parameters. For a change model, potential observations 

could include changes in head and changes in some flow quantity (such as ground-water discharge to a stream) that occur in 

response to stresses such as withdrawals by wells. Flow data 

for the C aquifer were insufficient for use in the change model; 

however, ground-water databases of the USGS and Arizona 

Department of Water Resources contain more than 3,000 water-level measurements in wells in the C aquifer. Some of these 

form time series measurements for individual wells that could 

be used to compute head change for parts of the period of historical withdrawals from 1961–2001. The distribution of potentially useful time series of water-level change includes well 

locations adjacent to the Little Colorado River east of Winslow, 

locations in the Flagstaff area, and some scattered locations in 

the unconfined area south of the Little Colorado River.

The change model was configured to simulate historical 

changes using the same parameter zones and aquifer-thickness distribution described earlier for the change model used to 

compute stream depletion. The model for simulating historical 

change, however, includes additional stream reaches, shown on 

figure 12, that may have been connected to the C aquifer during 

the period 1961–2001. The added streams were simulated with 

the River Package in the same manner described earlier for the 

five streams in the model used to calculate future stream depletions. In addition to the streams, two ground-water discharge 

areas, McDonald and Obed Meadows (fig. 12) were added to 

the model. Areas were delineated from USGS topographic 

maps and registered with the model grid to compute meadow 

area in each of 12 model cells intersected. For each of those 

cells, conductance was calculated as the product of that area, 

an assumed vertical hydraulic conductivity of 1 ft/d, divided by 

an assumed bed thickness of 20 ft. The resulting bed conductance ranged from 1.4×103

 to 9.7×104 ft2/d with an average of 

1.2×104 ft2/d.

19
Figure 11. Computed maximum drawdown in the 101-year model simulation period for scenarios A and B for the C aquifer, northeastern Arizona.

20

Clear

109o

110o

111o

112o

Flagstaff

Base from U.S. Geological Survey

digital data, 1:100,000, 1980

Lambert Conformal Conic projection

Standard parallels 2930' and 4530',

central meridian –11130'

N E W M E X I C O

A R I Z O N A

Chevelon

Creek

36o

35o

34o

Little

Colorado

River

0

0 25 MILES

25 KILOMETERS

EXPLANATION

1

NUMBERED PERENNIAL REACH

CONCEPTUAL WELL

EXISTING WELL USED IN WITHDRAWAL

 SCENARIO

3–Upper

 Clear

 Creek 5–Upper

 Chevelon

 Creek

4–Lower

 Chevelon

 Creek

2–Lower Clear Creek

1–Little Colorado

River below

Blue Spring

C AQUIFER CONFINED

LINE OF COMPUTED MAXIMUM

 DRAWDOWN IN THE 101-YEAR

 SIMULATION PERIOD—Values

 are in feet:

 Scenario A

1 Scenario B

Creek

01

1

1

1

1

01
Figure 12. Locations of withdrawal centers, stream reaches, and wells with observations in a change model for the C aquifer, northeastern 

Arizona, that simulates historical conditions.

21

Clear

109o

110o

111o

112o

Flagstaff

Winslow

APS

Joseph City

Holbrook

Hayhollow

SRP Concho

SRP Patterson

St. Johns St. Johns

Tepco

Springerville

Showlow

Eager

Pinetop-Lakeside

Woodruff

Abitibi

Snowflake

Taylor

Base from U.S. Geological Survey

digital data, 1:100,000, 1980

Lambert Conformal Conic projection

Standard parallels 2930' and 4530',

central meridian –11130'

N E W M E X I C O

A R I Z O N A

Chevelon

Creek

36o

35o

34o

(A-13-20)13ddd

(A-18-14)13abd

(A-18-15)28aad

Little

Colorado

River

0

0 25 MILES

25 KILOMETERS

EXPLANATION

NUMBERED PERENNIAL REACH USED FOR

 PROJECTIONS

NUMBERED STREAM REACH OR MEADOW

 USED IN THE CALIBRATION MODEL

WITHDRAWAL CENTER AND NAME

WELL WITH OBSERVATIONS USED IN OBJECTIVE

 FUNCTION OF MODEL USED TO SIMULATE

 HISTORICAL PERIOD AND NAME

3–Upper

 Clear

 Creek 5–Upper

 Chevelon

 Creek

4–Lower

 Chevelon

 Creek

2–Lower

 Clear

 Creek

6–Little Colorado River

between Chevelon Creek

and Clear Creek

9–Little Colorado River above the meadows

 10–Little Colorado River

 above the meadows

 13–Little Colorado River

 above Concho Creek

 14–Upper Little

 Colorado River

 11–Silver Creek

 12–Concho Creek

 7–McDonald Meadows

 8–Obed Meadows

1–Little Colorado

River below

Blue Spring

(A-18-14)13abd

Creek

Winslow
Table 2. Locations and rates of withdrawals from the C aquifer used in the Western Navajo Hopi Study change model to simulate 

historical conditions.

[Data source: Bill Greenslade, Southwest Ground-water Consultants, Inc. written commun., 2005]

Withdrawal rate, in acre-feet per year

Withdrawal center name Water use UTM easting 

(meters)

UTM northing 

(meters) 1961–1970 1971–1980 1981–1990 1991–2000

MANY FARMS Municipal and 

industrial

624036 4023766 333 333 466 746

CHINLE Municipal and 

industrial

630197 4001806 1133 1133 1587 2538

GANADO Municipal and 

industrial

632017 3942729 701 701 982 1571

EAGER Municipal and 

industrial

657570 3775603 650 650 990 1415

FLAGSTAFF Municipal and 

industrial

435464 3892934 480 1800 4934 7235

HOLBROOK Municipal and 

industrial

576866 3862326 1469 1469 1504 1539

PINETOP-LAKESIDE Municipal and 

industrial

595824 3778254 380 380 453 674

SHOW LOW Municipal and 

industrial

589344 3790577 1287 1287 1961 2801

SNOWFLAKE Municipal and 

industrial

584571 3819265 474 474 564 839

SPRINGERVILLE Municipal and 

industrial

658082 3778074 322 322 490 700

ST JOHNS Municipal and 

industrial

650498 3819268 100 100 210 281

TAYLOR Municipal and 

industrial

583426 3813899 350 350 417 621

WINSLOW Municipal and 

industrial

517320 3868734 1725 1725 1878 2494

TEPCO Industrial 660918 3794561 0 0 2484 8854

APS Industrial 563685 3858068 3500 3182 7754 13335

SRP CONCHO Industrial 629633 3821469 0 439 4705 5216

SRP PATTERSON Industrial 659160 3824565 0 370 3174 3461

ABITIBI Industrial 561051 3818273 8446 11430 15687 17415

JOSEPH CITY Agriculture 561013 3869422 3163 7762 7815 8950

WOODRUFF Agriculture 587813 3849446 253 621 625 716

HAY HOLLOW Agriculture 606106 3833548 506 1242 1250 1432

SNOWFLAKE Agriculture 584069 3816804 5568 13660 13754 15752

ST JOHNS Agriculture 627191 3815033 3163 7762 7815 8950

22
Examples of computed and observed water-level time series for three locations south of the Little Colorado River [A18-15)28aad, (A-18-14)13abd, and (A-13-20)13ddd] (fig. 12) 

show varying degrees of goodness of fit (fig. 13) using the parameters described earlier for the streamflow-depletion analysis. The first of these time series is for well (A-18-15)28aad, 

about 3 to 4 mi west of the start of the perennial part of lower 

Clear Creek. The time series starts in 1969, and the water level 

declines slightly more than 1 ft by the end of 2001. Observed 

head changes from 1982 to 1986 were not used because they 

were above zero and no process in the model could account 

for positive head changes. In the early 1980s, climate conditions were wetter than normal, and the positive values of head 

change in the well during and after this period are evidence that 

the aquifer responded to high runoff in lower Clear Creek and 

(or) other recharge. The model calculated about 5.7 ft of head 

decline at this location, significantly more than observed. The 

second time series, 1968–2001, is for well (A-18-14)13abd 

near the Winslow withdrawal center. Some of the individual 

observed values may have been influenced by recent and (or) 

nearby withdrawals, but the model reasonably matches the 

overall trend. The third time series, 1963–1985, is for well 

(A-13-20)13ddd near Snowflake. Observed declines are more 

than 40 ft over this period; the model underestimates the total 

decline at this well.

The entire procedure for the Constrained Monte Carlo 

analysis is as follows:

1. Generate N sets of parameters by perturbing values of 

horizontal-hydraulic conductivity for four zones, specific yield 

in the unconfined area, specific storage in the confined area, and 

riverbed conductance for all reaches.

2. For each of the N sets of parameters, run the model 

of historical conditions and compute the value of an objective function on the basis of the observed and computed headchange values shown in figure 13. Rank the N parameter sets 

from lowest to highest values of the objective function. Select 

the top M parameter sets on the basis of the lower objective 

function values.

3. Run the streamflow-depletion model simulating scenarios A and B with each of the top M parameter sets and compute the mean and standard deviation of the results of computed 

streamflow depletion in lower Clear and Chevelon Creeks.

The names of the seven model parameters allowed to vary 

in the Monte Carlo analysis are given in table 3. To generate the 

new parameters outlined in step 1 above, all but specific yield 

were log-transformed. For the log-transformed parameters, the 

equation for computing a perturbed value of the nth parameter, 

pn′, from the base value of the parameter, pn, is

 

where α1 is the i

th number in a series of normally distributed 

random numbers with a mean of 0.0 and a standard deviation of 

1.0 and Vn is the standard deviation of the nth parameter, in log10

space. The equation used for generating perturbed values of the 

untransformed parameter, specific yield (Sy), is

where Sy′ is the perturbed specific yield, and σSy is the standard 

deviation of specific yield. In generating the new parameters by 

using the above equations, values outside the upper and lower 

bounds given in table 3 were rejected and a new value was computed by using the next random number in the series. Values of 

and the upper and lower bounds in table 3 were selected on the 

basis of hydrologic judgment. 

For the first step, 1,829 parameter sets were generated and 

used in the model that simulates the historical period of withdrawals as outlined in step 2. The objective function for each 

of these sets was computed using the head-change observations 

from the three time series in figure 13. In computing the objective function, the first point in each series is the head value from 

which the head changes were calculated, and thus this point 

was not included in the objective function. The standard deviation of measurement error (input value STATdd in Hill and 

others, 2000) was set to 1.0 for the head-change observations 

for wells (A-18-14)13abd and (A-13-20)13ddd and was set to 

0.5 for the observations for well (A-18-15)28aad. This has the 

effect of increasing the weights for observations in the latter series, which is needed because the head-change values there are 

small relative to the changes in the two other series and would 

otherwise not contribute much to the objective function. The resulting values of the objective function ranged from 1.063×103

to 7.35×104. For comparison, the objective function for a simulation using the base parameter values is 2.405×103

. Many 

parameter sets resulted in lower objective-function values for 

these observations than the value from the base run using the 

parameter sets selected for the simulation of scenarios A and 

B. The use of the model to simulate historical changes at three 

sites, however, is not considered to be a rigorous calibration for 

selecting parameters. The process is instead a way of selecting 

the most reasonable parameter sets for use in the Monte Carlo 

analysis and of eliminating the worst sets in matching change 

observations at the three sites selected.

For step 3, the sets of parameters with objective functions 

less than 1.5 times the objective function from the base run 

were used. This resulted in 900 parameter sets with objective 

functions ranging from 1.063×103 to 3.604×103. For each run of 

the streamflow-depletion model, calculated depletions in lower 

Clear and lower Chevelon Creeks were saved at 5-year intervals over the 101-year simulation period with additional values 

saved at the end of the first time step, 1 year; the end of the period of withdrawals, 51 years; and at the end of the simulation, 

101 years. The saved values provide a data set for calculating 

basic statistics from the results for the entire simulation period. 

For each time interval saved, the mean and standard deviation 

of the 900 results were computed.

23

S , c

y S

y

DiVSy = +

p , c

n 10 10 log pn Di

Vn + © 1 § ·

=
Figure 13. Observed and computed head change for three wells 

completed in the C aquifer, northeastern Arizona.

24 HEAD CHANGE, IN FEET

0

1960 1965 1970 1975 1980 1985 2005

EXPLANATION

C. (A-13-20)13ddd

B. (A-18-14)13abd

A, (A-18-15)28aad

OBSERVED COMPUTED

1990 1995 2000

–5

–10

0

–5

–10

–15

–20

0

–5

–10

–15

–20

–25

–30

–35

–40

–45
The results for scenarios A and B in figure 14 show the 

upper and lower curves computed by adding and subtracting the 

standard deviation from the mean for each time interval saved. 

In all cases, the mean depletion from the 900 Monte Carlo runs 

is higher than the results for the model using the parameters 

selected for the depletion calculations. The first reason for 

the higher mean is that the model does not include a way for 

negative streamflow depletion to occur, and therefore the lower 

bound for streamflow depletion is zero. There is no corresponding bound for positive streamflow depletion, however, and the 

mean values from the Monte Carlo results therefore are skewed 

higher than the results using the standard parameters. The second reason is that the elimination of parameter sets with higher 

objective functions resulted in parameter distributions with 

median values that differ from the values of the corresponding 

parameters in the forward run. For example, the median of the 

900 HK1 values used in the Monte Carlo analysis is 9.3 ft/d, 

and the value of HK1 in the forward run is 6 ft/d. The higher 

median HK1 value leads to a higher mean depletion rate. In 

computing mean minus one standard deviation, the values less 

than zero in the lower curves on figure 14 were truncated at 

zero. The upper curves (mean plus one standard deviation) are 

substantially larger than the depletions simulated in the forward 

run, but the formulation of the Monte Carlo analyses does set 

upper bounds for reasonable parameter values. These curves 

can be viewed as a general measure of uncertainty in results 

using these limits and assumed standard deviations of parameters. Analyses of uncertainty could be improved with additional 

work in estimating parameters and determining parameter correlations, as well as in data collection to improve the model. The 

curves of estimated depletion in figure 10 should be used as the 

best estimates of computed depletion using the change model 

documented in this report.

Table 3. Parameters used in the Monte Carlo analysis of computed stream depletions in the study area, northeastern Arizona

[Standard deviations are in transformed space for HK1, HK2, HK3, HK4, SS1, and KRB and are in native space for SY1. Upper and lower bounds 

are in native space. — , unitless}

Parameter 

name Parameter description

Untransformed base 

value

Units Transformation type

Standard 

deviation, 

σ

Upper bound Lower bound

HK1 Horizontal hydraulic conductivity in zone 1 6 Per foot log10 0.5 60 0.6

HK2 Horizontal hydraulic conductivity in zone 2 5 Per foot log10 0.5 50 0.5

HK3 Horizontal hydraulic conductivity in zone 3 5 Per foot log10 0.5 50 0.5

HK4 Horizontal hydraulic conductivity in zone 4 5 Per foot log10 0.5 50 0.5

SY1 Specific yield in unconfined area 0.06 — none 0.02 0.12 0.01

SS1 Specific storage in confined area 2×10-6 Per foot log10 1 2×10-5 1×10-7

KRB Multiplier for riverbed conductance values for 

stream reaches used for projections

100 — log10 2 1000 0.1

25
Limitations of the Model

The ground-water change model of the C aquifer described 

in this report was designed specifically to compute the possible 

effects of ground-water withdrawals in an unconfined part of 

the aquifer near Leupp, Arizona, about 25 mi from the nearest 

connected surface-water feature, lower Clear Creek. Treatment 

of the aquifer as a porous medium with generalized aquifer 

properties is reasonable for this scale of simulation. This model 

should not be used for purposes such as evaluation of possible 

drawdown in and around well fields because local conditions 

such as flow in fractures and heterogeneities not represented in 

the model may be important at that scale.

The model also was not designed to evaluate the effects 

of existing withdrawals throughout the C aquifer on streams 

of interest including lower Clear and Chevelon Creeks. That 

purpose would require a calibrated flow model, rather than a 

change model. A related caution is that the model should not be 

used to evaluate the effects of withdrawals in areas other than 

near Leupp. The perimeter boundaries are distant from this area 

so that possible errors in placement or types of these boundaries 

will not affect the calculation of depletion in stream reaches of 

interest. This is not true, however, for an area such as Flagstaff, 

which is near a model boundary. Also, in the attempt to calibrate 

the change model, many observations near the confined–unconfined boundary could not be matched. More work would need 

to be done before the model could be used with confidence to 

evaluate the effects of withdrawals in those areas.

The two scenarios included simulations of 51 years of 

withdrawals followed by 50 years of no withdrawals. A key assumption was that the lengths of the perennial reaches remain 

the same over this period. This assumption is conservative in 

that the change model calculates the maximum total depletions for lower Clear and Chevelon Creeks. If combinations 

of the proposed withdrawals, other C-aquifer withdrawals, and 

extended drought were to cause the surface-water reaches to 

shrink or dry up, then depletion in lower Clear and Chevelon 

Creeks would be less than was indicated by the model. By the 

principle of conservation of mass, however, the total volume 

of outflow captured by the withdrawal scenarios would be the 

same—only the timing and locations of the captured outflow 

would vary.

Figure 14. Results of Monte Carlo analyses of computed stream depletions in the study area, northeastern Arizona.

Mean, upper and lower curves are mean results, mean plus one standard deviation, and mean minus one standard deviation, respectively. Forward-run curves are results from standard set of parameters, also shown on figure 10. A, Lower Clear Creek, scenario A. B, Lower Clear Creek, 

scenario B. C, Lower Chevelon Creek, scenario A. D, Lower Chevelon Creek, scenario B.

26

TIME, IN YEARS SINCE START OF SIMULATION TIME, IN YEARS SINCE START OF SIMULATION

0.4

0.3

0.2

0.1

0

Upper

C. Lower Chevelon Creek,

 scenario A

Mean

Lower

Forward run

STREAM DEPLETION, IN CUBIC FEET PER SECOND

0 20 40 60 80 100 120

D. Lower Chevelon Creek,

 scenario B

0.4

0.3

0.2

0.1

0

Upper

Mean

Lower

Forward run

B. Lower Clear Creek,

 scenario B

A. Lower Clear Creek,

 scenario A

1.5

1.0

0.5

0

Upper

Mean

Lower

Forward

run

1.5

1.0

0.5

0

Upper

Mean

Lower

Forward run

0 20 40 60 80 100 120
Summary and Conclusions

A ground-water change model was constructed for the 

purpose of evaluating possible effects of withdrawals from the 

C aquifer in the area of Leupp, Arizona, on connected surfacewater features including Clear Creek, Chevelon Creek, and the 

Little Colorado River below Blue Spring. The model was constructed using the USGS model code MODFLOW-2000. The 

approximate aquifer thickness was mapped and four zones were 

used for the distribution of hydraulic conductivity. The model 

includes confined and unconfined areas where storage properties are represented with values of the product of specific storage and thickness, and specific yield, respectively. A limited 

calibration was attempted using historical withdrawals and observed head changes. The final model, however, used values of 

hydraulic conductivity and other properties that were deemed 

representative of large zones in the aquifer.

Evaluations of possible stream depletion included evaluations of scenarios A and B, with maximum withdrawal rates of 

about 6,500 and 11,500 acre-ft/yr, respectively. The scenarios 

included a 51-year period of withdrawals followed by a 50-

year period of no withdrawals. The computed volume of stream 

depletion for the 101-year period was about 6 percent of the 

ultimate volume that would occur as a result of the scenarios. 

For both scenarios, computed maximum depletion rates in upper Clear Creek, upper Chevelon Creek, and the Little Colorado 

River below Blue Spring were on the order of hundredths of a 

cubic foot per second; detailed time series of these results were 

not presented. Most of the computed depletion was in lower 

Clear Creek; computed rates reached about 0.4 ft3/s for scenario 

A and about 0.5 ft3/s for scenario B near the end of the simulation. Computed depletion in lower Chevelon Creek is less 

than that in lower Clear Creek because Chevelon Creek is more 

distant from the withdrawal locations, and drawdown reaches 

Clear Creek first.

A Monte Carlo analysis was carried out to compute statistical properties of streamflow depletion using a range of parameter sets. The analysis was constrained by using a change model 

that simulated a period of historical withdrawals from 1961 to 

2001. An objective function that used observations from 3 relevant well locations was computed for each of 1,829 sets of 

7 model parameters, and the 900 sets with the best objective 

function were used to calculate a set of streamflow depletions 

by using the model that simulated the proposed withdrawals 

for scenarios A and B. Statistical results show that the mean 

and mean plus and minus standard deviations are skewed to the 

high side because the streamflow is bounded on the low side 

by zero.

The change-model approach seems to be reasonable for 

calculating the possible effects of withdrawals near Leupp. 

The model is not designed or calibrated for other purposes. 

Evaluation of the cumulative effects of the proposed withdrawals and existing withdrawals would best be done with a 

calibrated ground-water flow model.

27
References

Bills, D.J., Truini, Margot, Flynn, M.E., Pierce, H.A., Cathings, R.D., and Rymer, M.J., 2000, 

Hydrogeology of the regional aquifer near Flagstaff, Arizona, 1994–97: U.S. Geological Survey 

Water-Resources Investigations Report 00-4122, 142 p.

Blakey, R.C., 1990, Stratigraphy and geologic history of Pennsylvanian and Permian rocks, Mogollon Rim 

region, central Arizona and vicinity: Geological Society of America Bulletin, v. 102, p. 1189–1217.

Brown, D.E., Carmony, N.B., and Turner, R.M., compilers, 1981, Drainage map of Arizona showing 

perennial streams and some important wetlands: Phoenix, Arizona Game and Fish Department, scale 

1:100,000.

Cooley, M.E., Harshbarger, J.W., Akers, J.P., and Hardt, W.F., 1969, Regional Hydrogeology of the Navajo 

and Hopi Indian Reservations, Arizona, New Mexico, and Utah: U.S. Geological Survey Professional 

Paper 521-A, 61 p.

Daniel, Debra, 1981, Total dissolved solids content of the aquifers of the Basin and Range and the C-aquifer, Arizona Department of Water Resources, 1 pl.

Fenneman, N.M., 1946, Physical divisions of the United States: U.S. Geological Survey map, scale 

1:700,000.

Harbaugh, A.W., Banta, E.R., Hill, M.C., and McDonald, M.G., 2000, MODFLOW-2000, the U.S. 

Geological Survey modular ground-water model—User guide to modularization concepts and the 

Ground-Water Flow Process: U.S. Geological Survey Open-File Report 00–92, 121 p.

Hart, R.J., Ward, J.J., Bills, D.J., and Flynn, M.E., 2002, Generalized hydrogeology and ground-water 

budget for the C Aquifer, Little Colorado River basin and parts of the Verde and Salt River basins, 

Arizona and New Mexico: U.S. Geological Survey, Water-Resources Investigations Report 02–4026, 

54 p.

Hill, M.C., Banta, E.R., Harbaugh, A.W., and Anderman, E.R., 2000, MODFLOW-2000, the U.S. 

Geological Survey modular ground-water model—User guide to the Observation, Sensitivity, and 

Parameter Estimation Processes and three post-processing programs: U.S. Geological Survey OpenFile Report 00–184, 209 p.

Hoffmann, J.P., Phillips, J.V., Bills, D.J., and Halford, K.J., in press, Geologic, hydrologic, and chemical data from the C aquifer, near Leupp, Arizona: U.S. Geological Survey Scientific Investigations 

Report 2005–5280.

Leake, S.A. and Claar, D.V., 1999, Procedures and computer programs for telescopic mesh refinement using MODFLOW: U.S. Geological Survey Open-File Report 99–238, 53 p.

Mann, L.J., 1976, Ground-water resources and water use in southern Navajo County, Arizona: Phoenix, 

Arizona Water Commission Bulletin 10, 106 p.

Mann, L.J., and Nemecek, E.A. 1983, Geohydrology and water use in southern Apache County, Arizona: 

Phoenix, Arizona Department of Water Resources Bulletin 1, 86 p.

McDonald, M.G. and Harbaugh, A.W., 1988, A modular three-dimensional finite-difference ground-water 

flow model: U.S. Geological Survey Techniques of Water-Resources Investigations, Book 6 Modeling 

Techniques, Chap. A1, 586 p.

McDonald, M.G. and Harbaugh, A.W., 1996, User’s documentation for MODFLOW-96, an update to the 

U.S. Geological Survey modular finite-difference ground-water flow model: U.S. Geological Survey 

Open-File Report 96–485, 56 p.

McGavock, E.H., Anderson, T.W., Moosburner, Otto, and Mann, L.J., 1986, Water resources of southern 

Coconino County, Arizona: Phoenix, Arizona Department of Water Resources Bulletin 4, 53 p.

Peter Mock Groundwater Consulting, Inc., 2003, in HDR Engineering, Inc. Report of Findings for the 

Assessment of Western Navajo-Hopi Water Supply Needs, Alternatives, and Impacts, Volume 3, 

Appendix B-3, variously paged.

28
Reilly, T.E., Franke, O. Lehn, and Bennett, G.D., 1987, The principle of superposition and its application 

in ground-water hydraulics: U.S. Geological Survey Techniques of Water-Resources Investigations 

Book 3, Chap. B6, 28 p.

Southwest Ground-water Consultants, Inc., 2003, in Report of findings for the assessment of Western 

Navajo-Hopi water supply needs, alternatives, and impacts, v. 3, Appendix B-1: HDR Engineering, Inc 

for U.S. Bureau of Reclamation, 2003, 10 p.

Theis, C.V., 1935, The relation between lowering of the piezometric surface and the rate and duration of 

discharge of a well using ground-water storage: Transactions of the American Geophysical Union, 

v. 16, p. 519–524.

U.S. Geological Survey, 2005, Elevation Derivatives for National Applications, accessed June, 15, 2004, at 

URL http://edna.usgs.gov/

29
This page intentionally left blank; 

inside back cover.
Leake and others—Numerical Ground-Water Change Model, Northeastern Arizona—SIR 2005–5277

39