Ground-Water Flow Model of the Sierra Vista 

Subwatershed and Sonoran Portions of the 

Upper San Pedro Basin, Southeastern Arizona, 

United States, and Northern Sonora, Mexico

Prepared in cooperation with the

UPPER SAN PEDRO PARTNERSHIP and BUREAU OF LAND MANAGEMENT

Scientific Investigations Report 2006–5228

U.S. Department of the Interior

U.S. Geological Survey

ARIZONA, UNITED STATES

SONORA, MEXICO

Sierra Vista

P

San de ro R vi er

Cananea

Tombstone

Palominas

Bisbee

Fort Huachuca
Inside Front Cover 

BLANK
Ground-Water Flow Model of the Sierra 

Vista Subwatershed and Sonoran 

Portions of the Upper San Pedro Basin, 

Southeastern Arizona, United States, 

and Northern Sonora, Mexico

By D.R. Pool and Jesse E. Dickinson

Prepared in cooperation with the

UPPER SAN PEDRO PARTNERSHIP and the 

U.S. BUREAU OF LAND MANAGEMENT

Scientific Investigations Report 2006–5228

U.S. Department of the Interior

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

DIRK KEMPTHORNE, Secretary

U.S. Geological Survey

Mark D. Myers, Director

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

For product and ordering information: 

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

Telephone: 1-888-ASK-USGS

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

natural hazards, and the environment: 

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

Telephone: 1-888-ASK-USGS

Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement 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.

Suggested citation:

Pool, D.R., and Dickinson, J.E., 2007, Ground-water flow model of the Sierra Vista Subwatershed and Sonoran 

portions of the Upper San Pedro Basin, southeastern Arizona, United States, and northern Sonora, Mexico: 

U.S. Geological Survey Scientific Investigations Report 2006–5228, 48 p.
iii

Contents

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

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

Purpose and Scope ..............................................................................................................................1

Approach................................................................................................................................................1

Description of Study Area ...................................................................................................................2

Climate...........................................................................................................................................2

Land and Water Use....................................................................................................................2

Acknowledgments................................................................................................................................4

Conceptual Model of the Ground-Water Flow System............................................................................4

Hydrogeologic Framework..................................................................................................................4

Aquifers ..................................................................................................................................................4

Pre-Tertiary Rocks .......................................................................................................................7

Tertiary Pre-Basin Fill Sediments..............................................................................................7

Basin Fill ........................................................................................................................................8

Stream Alluvium...........................................................................................................................8

Ground-Water Flow System................................................................................................................9

Ground-Water Budget..........................................................................................................................9

Predevelopment.........................................................................................................................13

Postdevelopment .......................................................................................................................14

Artificial Recharge............................................................................................................14

Evapotranspiration ...........................................................................................................14

Ground-Water Withdrawals............................................................................................15

Simulation of Ground-Water Flow.............................................................................................................18

Previous Models .................................................................................................................................18

Model Framework...............................................................................................................................18

Model Boundaries .....................................................................................................................18

Spatial and Temporal Discretization.......................................................................................19

Aquifer Hydraulic Properties ...................................................................................................19

Layer Distribution..............................................................................................................19

Transmissive Properties ..................................................................................................20

Storage Properties ...........................................................................................................20

Simulated Inflows and Outflows ......................................................................................................26

Recharge.....................................................................................................................................26

Natural Recharge .............................................................................................................26

Artificial and Incidental Recharge.................................................................................26

Evapotranspiration of Ground Water......................................................................................29

Streamflow Routing...................................................................................................................30

Ground-Water Withdrawals.....................................................................................................30

Initial Conditions..................................................................................................................................31

Model Calibration................................................................................................................................32

Steady State ...............................................................................................................................32

Average-Annual Steady-State Conditions....................................................................32

Oscillatory Steady-State Conditions..............................................................................35
iv

Contents—Continued

Transient-State Conditions.......................................................................................................35

Model Limitations and Possible Improvements.............................................................................42

Summary and Conclusions.........................................................................................................................44

References Cited..........................................................................................................................................45

Figures

1. Model area within the Upper San Pedro Basin, Arizona, United States 

and Sonora, Mexico .....................................................................................................................3

2. Geology and saturated thickness of the regional aquifer, Upper 

San Pedro Basin, United States and Mexico...........................................................................6

3. Generalized hydrogeologic section and extent of numerical model layers, 

Upper San Pedro Basin, United States and Mexico...............................................................7

4. Extent and altitude of thick silt and clay interval in the basin fill, 

Upper San Pedro Basin, United States and Mexico. 4A, Top of 

thick silt and clay interval; 4B, Base of thick silt and clay 

interval..........................................................................................................................................10

5. Ground-water flow system, boundary conditions for the model, streamflow 

routing network, evapotranspiration areas, and predevelopment 

water-level comparison sites, Upper San Pedro Basin, United States 

and Mexico ..................................................................................................................................12

6. Annual ground-water withdrawal rates for selected uses, 1902 to 2002, 

Upper San Pedro Basin, United States and Mexico. 6A, Mine uses; 

6B, Non-mine uses; 6C, Total uses...........................................................................................16

7. Extent and lithology of hydrogeologic units and altitude of the base of 

numerical model layers, Upper San Pedro Basin, United States and 

Mexico. 7A, Model layer 5; 7B, Model layer 4; 7C, Model layer 3; 

7D, Model layer 2; 7E, Model layer 1 .......................................................................................21

8. Simulated total transmissivity of the numerical model layers, Upper 

San Pedro Basin, United States and Mexico.........................................................................27

9. Pre-1961 water-level hydrographs for the Sierra Vista subwatershed of 

the Upper San Pedro Basin, United States ............................................................................33

10. Observed and simulated predevelopment water-level altitude error for 

selected monitoring wells, Upper San Pedro Basin, United States and 

Mexico. 10A, Observed and simulated predevelopment water-level 

altitudes; 10B, Residual of observed and simulated predevelopment 

water-level altitudes...................................................................................................................34

11. Distributions of simulated withdrawal wells and selected hydrograph 

comparison sites, Upper San Pedro Basin, United States and Mexico............................38

12. Hydrographs showing observed and simulated water-level altitudes 

at selected wells, 1940–2003, Upper San Pedro Basin, United States 

and Mexico ..................................................................................................................................39

13. Estimated and simulated seasonal baseflow at streamflow-gaging stations 

on the San Pedro River, Arizona, 1929–2003. 13A, Charleston; 13B, Palominas...............41

14. Observed and simulated water-level altitudes, winter 2002, Upper 

San Pedro Basin, United States and Mexico.........................................................................43
Tables

1. Correlation of hydrogeologic units and model layers.............................................................5

2. Estimated predevelopment ground-water budget................................................................13

3. Simulated hydraulic properties of hydrogeologic units for the ground-water 

flow model of the Upper San Pedro basin, United States and Mexico .............................28

4. Simulated predevelopment ground-water budget................................................................36

5. Simulated ground-water budget for March 2002 to March 2003 .......................................44

Conversion Factors and Datum

Multiply By To obtain

centimeter (cm) 0.3937 inch

meter (m) 3.281 foot

kilometer (km) 0.6214 mile 

square meter (m2

) 0.0002471 acre 

square kilometer (km2

) 0.3861 square mile 

cubic hectometer (hm3

) 810.7 acre-foot 

cubic hectometer per year (hm3

/yr) 811.03 acre-foot per year

cubic meter per day (m3

/yr) 3.378 acre-foot per year

cubic meter per second (m3

/s) 35.31 cubic foot per second

square hectometer (hm2

) 2.471 acre

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

1988 (NAVD 88). 
vi

List of Abbreviations

DEM Digital elevation model

ENSO El Niño Southern Oscillation

ft bls feet below land surface

GIS Geographic Information System

Landsat Landsat satellites

LIDAR Light Detection And Ranging

MF2K MODFLOW 2000

ohm-m ohm-meters

PDO Pacific Decadal Oscillation

SPRNCA San Pedro Riparian Conservation Area

STR1 Streamflow Routing package for MODFLOW

USDA U.S. Department of Agriculture

USGS U.S. Geological Survey

BLM Bureau of Land Management

Arizona Well-Numbering System

6 5 4 3

7 8 11 12

18

19 20 22 24

30 29

31 32 33 34 35 36

b

9

21 c

d

a

c d

d

17

a

b

Quadrant D, Township 21 North, Range 21 East, section 22,

quarter section d, quarter section d, quarter section c

R. 21 E.

T.

21

N.

WELL (D–21–21)22ddc

ARIZONA WELL-NUMBERING SYSTEM

A

C D

B

The well numbers used by the U.S. Geological Survey in Arizona are in accordance with the Bureau of Land 

Management's system of land subdivision. The land survey in Arizona is based on the Gila and Salt River 

meridian and base line, which divide the State into four quadrants and are designated by capital letters A, 

B, C, and D in a counterclockwise direction beginning in the northeast quarter. The first digit of a well 

number indicates the township, the second the range, and the third the section in which the well is situated. 

The lowercase letters a, b, c, and d after the section number indicate the well location within the section. 

The first letter denotes a particular 160-acre tract, the second the 40-acre tract, and the third the 10-acre 

tract. These letters are also assigned in a counterclockwise direction beginning in the northeast quarter. If 

the location is known within the 10-acre tract, three lowercase letters are shown in the well number. When 

more than one well is within a 10-acre tract, consecutive numbers beginning with 1 are added as suffixes. 

In the example shown, well number (D–21–21)22ddc designates the well as being in the SE1/4, SE1/4,

SW1/4, section 22, Township 21 North, and Range 21 East.

2

14

23

25

10

1

16

28 27

22

c

D b a

15 13

26 d
Ground-Water Flow Model of the Sierra Vista and Sonoran 

Portions of the Upper San Pedro Basin, Southeastern 

Arizona, United States, and Northern Sonora, Mexico 

By D.R. Pool and Jesse E. Dickinson

Abstract

A numerical ground-water model was developed to 

simulate seasonal and long-term variations in ground-water 

flow in the Sierra Vista subwatershed, Arizona, United States, 

and Sonora, Mexico, portions of the Upper San Pedro Basin. 

This model includes the simulation of details of the groundwater flow system that were not simulated by previous models, 

such as ground-water flow in the sedimentary rocks that 

surround and underlie the alluvial basin deposits, withdrawals 

for dewatering purposes at the Tombstone mine, discharge to 

springs in the Huachuca Mountains, thick low-permeability 

intervals of silt and clay that separate the ground-water flow 

system into deep-confined and shallow-unconfined systems, 

ephemeral-channel recharge, and seasonal variations in 

ground-water discharge by wells and evapotranspiration.

Steady-state and transient conditions during 1902–2003 

were simulated by using a five-layer numerical groundwater flow model representing multiple hydrogeologic 

units. Hydraulic properties of model layers, streamflow, 

and evapotranspiration rates were estimated as part of 

the calibration process by using observed water levels, 

vertical hydraulic gradients, streamflow, and estimated 

evapotranspiration rates as constraints. Simulations 

approximate observed water-level trends throughout most 

of the model area and streamflow trends at the Charleston 

streamflow-gaging station on the San Pedro River. Differences 

in observed and simulated water levels, streamflow, and 

evapotranspiration could be reduced through simulation of 

climate-related variations in recharge rates and recharge from 

flood-flow infiltration. 

Introduction

Increased water use has increased demand on the water 

supply in the Sierra Vista subwatershed of the Upper San 

Pedro Basin. Ground water is the primary source of water 

in the basin for many uses, including riparian, agricultural, 

private, public, and industrial; and supply for the military 

installation at Fort Huachuca. Baseflows in perennial-stream 

reaches of the San Pedro River are sustained by ground-water 

flow from upgradient recharge areas. Increasing demand 

on the ground-water supply likely will result in decreasing 

amounts of ground-water available to baseflow and riparian 

vegetation. Basic water-budget analysis shows that groundwater withdrawals for upgradient uses will eventually result 

in a reduction in discharge from the basin owing to reduced 

baseflow and evapotranspiration by plants. The amount of the 

reduced discharge will be equivalent to the amount withdrawn 

from the system assuming that inflow to the ground-water 

system does not change. The rate and location of reduced 

discharge, however, is not well known because of a lack of 

information about the hydrogeologic system in the basin. 

Improved knowledge of interactions between the stream and 

ground-water systems is needed to make informed decisions 

about the management of the water resources in the basin.

Purpose and Scope

A new numerical ground-water flow model of the study 

area was constructed to simulate new conceptualizations of 

the ground-water flow system in the Upper San Pedro Basin. 

Hydrologic investigations and monitoring during the mid1990s–2004 have improved understanding of the groundwater flow system beyond that simulated the by previous 

models. The new model can be used to aid management 

of the resources and to assess the impacts of ground-water 

withdrawals and artificial recharge on streamflow and 

riparian vegetation in the basin. This report documents the 

development and calibration of the new numerical groundwater flow model. 

Approach

Information from extensive hydrologic investigations 

and monitoring by multiple agencies and investigators during 

the mid-1990s–2004 was used to refine the conceptual 

model of the ground-water flow system in the study area. 

Investigations provided new information to better define 

the spatial and temporal distribution of recharge, vertical 
hydraulic-head gradients, interactions between the stream 

and the aquifer, and the subsurface distribution of silt and 

clay layers in the primary aquifer. Definition of the spatial 

and temporal distribution of recharge was improved by new 

data from several streamflow-gaging stations, unsaturatedzone monitoring beneath several ephemeral streams and 

in interdrainage areas, monitoring ground-water levels and 

ground-water storage, and geochemical analysis of runoff 

and ground water in recharge areas. Definition of vertical 

hydraulic-head gradients near the San Pedro River and 

stream–aquifer interactions was improved by using water 

levels measured bimonthly or more frequently at multiple 

piezometers. Geophysical investigations were used to augment 

drill-log information defining the distribution of silt and clay 

layers in the primary aquifer. Rates of evapotranspiration 

for several important vegetation types were better defined 

for the period 2002-03. Much of the new information was 

acquired in the Arizona portion of the watershed. New 

information in Mexico also was provided by recent geologic 

and hydrogeologic investigations and ground-water flow 

models. New and existing information was integrated to 

develop a multilayered, numerical model of steady-state and 

transient ground-water flow conditions. Hydraulic and flow 

properties of the numerical model were constrained by using 

observations of hydraulic head and streamflow.

Description of Study Area

The study area includes the upper drainage area of 

the San Pedro River in northern Sonora, Mexico, and the 

Sierra Vista subwatershed in southeastern Ariz., United 

States (fig. 1), and includes about 4,500 km2 extending from 

the southern boundary near Cananea, Mexico, to about 

Fairbank, Ariz. The watershed is bounded by mountains and 

intervening broad alluvial passes. Mountains that bound the 

study area on the west include the Mustang Mountains and 

Huachuca Mountains, but most of both mountain ranges are 

completely within the San Pedro River drainage and study 

area. Mountains to the east include the Mule Mountains, 

Tombstone Hills, and Sierra San Jose. Mountains to the south 

include the Sierra Mariquita and Sierra Los Ajos. Altitude of 

the watershed ranges from about 1,052 m at the streamflowgaging station near Tombstone, Ariz., to more than 2,700 m in 

the Huachuca Mountains. 

The major streams are intermittent within the study 

area and include the San Pedro River, Babocomari River, 

and Rio Los Fresnos. The San Pedro River flows perennially 

near Hereford, Ariz., and from south of Arizona Highway 90 

into the Boquillas area (fig. 1), but perennial reaches prior to 

development extended south of Palominas. Major perennial 

reaches along the Babocomari River occur between the 

Mustang Mountains and Huachuca Mountains and near the 

confluence with the San Pedro River. Perennial reaches occur 

along Rio Los Fresnos. Other tributary streams are ephemeral. 

Major ephemeral streams include Walnut Gulch near 

Tombstone and Greenbush Draw between the Mule Mountains 

and Sierra San Jose.

Streamflow and shallow ground water support a 

narrow corridor of riparian vegetation along the floodplain 

of the intermittent streams. Riparian vegetation includes 

phreatophytes, such as cottonwood, willow, and mesquite. 

The riparian area is a valued resource that supports several 

endangered species and is an important migratory bird habitat. 

Much of the riparian area in the United States has been 

protected by designation as the San Pedro Riparian National 

Conservation Area (SPRNCA) and is managed by the U.S. 

Bureau of Land Management. 

Climate

Climate of the study area is arid to semiarid. Precipitation 

data are available during 1905–2002 from a gage at 

Tombstone, at an altitude of 1,384 m, and during 1956–2002 

from a gage at the southern end of the Huachuca Mountains 

(Coronado National Memorial), at an altitude of 1,692 m 

(fig. 1). Historically, the wettest season is July through 

September, averaging 21 cm of precipitation at the Tombstone 

gage and 27 cm of precipitation at the Huachuca Mountains 

gage. Precipitation during October through March averages 

12 and 23 cm at the Tombstone and Huachuca Mountain 

gages, respectively. Precipitation during April through June 

averages less than 3 cm at both gages. Long-term monitoring 

indicates that annual precipitation rates were above average 

before about 1940 and during the early to mid-1980s, below 

average from about 1940 to 1980, and about average after 

the mid-1980s (Pool and Coes, 1999). Trends in seasonal 

precipitation, however, are different from trends in annual 

precipitation. “Winter” (November–February) precipitation 

and streamflow varies with identifiable multi-year and 

decadal-scale fluctuations in temperature and precipitation 

(Webb and Betancourt, 1992; Dickinson and others, 2004; 

Pool, 2005). High rates of winter precipitation tend to occur 

during El Niño periods (Pool, 2005). Decadal-scale variations 

include above-average winter precipitation rates before 

1940 and during 1956–97. “Summer” (July–September) 

precipitation rates and streamflow runoff generally decreased 

during the later half of the 1900s (Pool and Coes, 1999; 

Thomas and Pool, 2006).

Land and Water Use

Water use in the study area changed from use 

primarily for livestock grazing and mining before 1940 to 

primarily agricultural use during the mid-1900s and mining, 

industrial, and domestic use during the late 1900s. Some 

of the water that was derived from dewatering of the mine 

workings beginning in 1902 at Tombstone (Hollyday, 1963) 

and in 1906 at Bisbee, Ariz., was used for domestic, industrial, 

and irrigation purposes. The remainder of the water derived 

from the mines was allowed to evaporate or recharge the 

aquifer. Dewatering at the Tombstone mine ceased in 1911, but 

continued at the Bisbee mine until 1986. Extensive water use 

began at the Cananea mine during the mid-1980s (Consultores 

en Agua Subterranea S.A. por Mexicana de Cananea, S.A. 

de C.V., 2000). The earliest agriculture likely began with 

the irrigation south of Bisbee by using water derived from 

mine dewatering. Agricultural water use near the San Pedro 

River began before 1935 near Hereford by using supplies 

that were derived from streamflow and wells near the river. 

 Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Figure 1. Model area within the Upper San Pedro Basin, Arizona, United States and Sonora, Mexico.

Introduction 

30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

PIMA COUNTY 

SANTA CRUZ 

COUNTY 

Figure 1. Map showing model area within the Upper San Pedro Basin, Arizona, United States, and Sonora, Mexico.

ARIZONA

SONORA

UNITED STATES

MEXICO

Cananea

Bisbee

Tombstone

B a ri

abocom R vi

R o

nos o L s

e s r F i

er

Palominas

Naco

Huachuca City 

MULE

MOUNTAINS

SIERRA

MARIQUITA

Fort Huachuca

East

Range

COCHISE

COUNTY

110o30’

15’

31o

110o

15’

Hereford 

Lewis Springs

San Pedro

Riparian

National

Conservation

Area

Charleston

BEN ON

TE

S

S UB RSHED

WA

TE

S UB RSHED

WA

TERSH D

NORTHERN BOUNDARY

OF GROUND-WATER

FLOW MODEL

BOUNDARY

E

WA

P

de

ro

R

vi

er

Sa

n

HUACHUCA

MOUNTAINS

90

80

92

82

Elgin 

Fairbank 

TOMBSTONE HILLS

MUSTANG

MOUNTAINS

DRAGOON

MOUNTAINS

SIERRA SAN JOSE

SIERR

LO

A

S AJOS

2

UNITED STATES

MEXICO

River

San edro P

ARIZONA

Model area

Upper

San Pedro

Basin

Tucson

Sierra Vista

Subwatershed

WHETSTONE

MOUNTAINS

Tombstone

precipitation gage

Sierra Vista

Streamflowgaging station

at Palominas

Coronado National

Memorial

precipitation gage

Streamflow-gaging

station at Charleston

Streamflow-gaging

station near Tombstone

SIERRA TA V IS

Copper Queen

mine

Tombstone mines

Gr

W Gulch

alnut

r

ee us n

D

w ra

b h

Boquillas

Cottonwood

EXPLANATION

BOUNDARY OF GROUND-WATER

 FLOW MODEL
Agricultural water use near the San Pedro River in the Sierra 

Vista subwatershed increased throughout the 1970s in areas 

south of Highway 90. Agriculture was eliminated in the area 

north of Hereford during the mid-1980s and further decreased 

after 1998 to only a few irrigated fields in the Palominas 

area. Agriculture in the Mexico portion of the study area 

is not well documented; however, satellite remote-sensing 

information (Kepner and Edmonds, 2002) indicates that 

several agricultural areas existed near the Rio Los Fresnos and 

San Pedro River in 1973. Industrial and domestic water use 

began to be significant as the population increased in the Fort 

Huachuca and the Sierra Vista areas after the mid-1960s and 

near Cananea beginning in the 1980s.

Acknowledgments

Member agencies and individuals of the Upper San Pedro 

Partnership were helpful in collecting and disseminating 

information used to construct the ground-water flow model. 

Individual knowledge of historical water use was often 

useful in developing model input when documentation was 

insufficient. The Morris K. Udall Center was instrumental in 

arranging for the dissemination of hydrogeologic information 

that was collected in support of mining operations in 

Mexico. Recent studies and data made available by the U.S. 

Department of Interior, U.S. Bureau of Land Management; 

University of Arizona Department of Hydrology; and U.S. 

Department of Agriculture, Agricultural Research Service 

were used in the model construction. The Phelps Dodge 

Corporation made available hydrogeologic and ground-water 

withdrawal information for areas near the Copper Queen mine 

near Bisbee. 

Conceptual Model of the Ground-Water 

Flow System

A conceptual model of the ground-water flow system is 

a general understanding of where and at what rate recharge 

occurs, how water flows through the aquifer, and where water 

discharges. Though conceptual models are constrained by 

observations of water levels in wells, discharge rates and 

distributions, and estimates of recharge rates and distributions, 

observed ground-water conditions may be explained by many 

different conceptual models. Conceptual models evolve with 

improved knowledge about ground-water flow systems and 

provide a basis for construction of numerical models.

Previous concepts of ground-water flow in the 

Upper San Pedro Basin were simple because of limited 

data availability. A more complex conceptual model was 

developed as more details of the system became known. 

Important new advances include improved quantification 

of evapotranspiration rates and distributions and a better 

understanding of the spatial distribution of recharge, vertical 

ground-water flow, stream-aquifer interactions, spatial 

distribution of silt and clay in the alluvial aquifer, the influence 

of natural variations in recharge and discharge, and response 

of the ground-water flow system to variations in recharge and 

discharge. A general description of the conceptual groundwater flow system includes descriptions of the aquifers and 

quantification of recharge rates and distributions, groundwater flow to the discharge areas, and discharge rates and 

distributions.

Hydrogeologic Framework

The hydrogeologic framework includes descriptions 

of the thickness, distribution, and hydraulic properties of 

water-bearing rocks within a ground-water flow system. 

The hydrogeologic framework of the Upper San Pedro Basin 

has been described in several reports (Brown and others, 

1966; Pool and Coes, 1999; Fleming and Pool, 2002; Condor 

Consulting Inc., 2003; and Coes and Pool, 2005). The primary 

aquifer is as much as 520 m of basin fill, alluvial deposits that 

overlie relatively impermeable crystalline and sedimentary 

rocks (table 1, fig. 2 and fig. 3). Secondary aquifers of local 

importance include sedimentary rocks that crop out in the 

mountains and underlie the basin fill and stream alluvium 

along the major streams. Ground water generally flows 

from recharge areas near the mountains, through sand and 

gravel layers in the basin fill, and toward the San Pedro 

or Babocomari Rivers. Ground water discharges near the 

San Pedro and Babocomari Rivers as baseflow, near the 

San Pedro River as small springs, as ground-water outflow, 

and as evapotranspiration in riparian areas. A portion of the 

ground-water flow is intercepted upgradient of the streams by 

pumping wells and phreatophytes where depths to water are 

shallow along ephemeral streams. 

Aquifers

The primary regional aquifer includes upper and 

lower basin fill, described by Brown and others (1966), that 

accumulated in the structural depression between mountain 

ranges during the late Tertiary and early Pleistocene (table 

1, fig. 2 and fig. 3). The saturated thickness of the regional 

aquifer is 200–400 m within two structural depressions. 

The largest depression occurs throughout much of the extent 

of the basin fill in the United States. A smaller structural 

depression is in Mexico between Sierra San Jose and Sierra 

Mariquita (fig. 1). The saturated thickness of the regional 

aquifer is less than 100 m in small structural basins that 

coincide with the upper parts of Walnut Gulch, Babocomari 

River, and Greenbush Draw drainage basins (fig. 1). 

Secondary aquifers include Quaternary terrace and 

alluvial deposits that are less than about 10 m thick and 

generally coincide with the floodplains of the San Pedro and 

Babocomari Rivers and tributary streams (Pool and Coes, 

1999). Tertiary prebasin fill sediment and Mesozoic and 

Paleozoic limestone that crop out in the mountains also are 

secondary aquifers in places. Other rocks that crop out in the 

mountains and hills surrounding the basin are insignificant 

aquifers and include Tertiary and older granitic and volcanic 

rocks and Mesozoic sedimentary rocks consisting of 

mudstone, quartzite, and conglomerate. 

 Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Table 1. Correlation of hydrogeologic units and model layers.

Hydrogeologic unit Lithologic description

Thickness, in meters1

Model 

layer

Range in 

thickness of 

individual unit

Range in 

thickness of 

combined 

model layers

Post-entrenchment stream alluvium sand and gravel 0–10

0–520

Pre-entrenchment stream alluvium clay, silt, and fine sand 0–10 1

Upper basin fill sand and gravel facies 0–100

sand and gravel facies, and silt and clay facies 0–300 2

Lower basin fill sand and gravel facies, and siltstone and mudstone facies 0–170 3

sand and gravel facies 0–400 4

sand and gravel facies on the perimeter of the alluvial basin 1,500+

1,500 5 Pantano formation siltstone and conglomerate 1,500+

Consolidated rocks limestone, sandstone, mudstone, quartzite, conglomerate, 

granite, and volcanic rocks

1,500+

1

Thickness includes the unsaturated material between the land surface and top of the upper-most layer.

Conceptual Model of the Ground-Water Flow System 
Figure 2. Geology and saturated thickness of the regional aquifer, Upper San Pedro Basin, United States and Mexico.

0 10 KILOMETERS

0 6 MILES

San

Geology modified from Reynolds, 1988

30’

31o45’

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

ARIZONA

SONORA

UNITED STATES

MEXICO

Sierra Vista

Cananea

Tombstone

110o30’

15’

15’

31o

110o

TERSH D BOUNDARY

E

WA

GEOLOGY

SATURATED THICKNESS OF THE

REGIONAL ALLUVIAL AQUIFER, IN METERS

Volcanic rock

Sandstone, siltstone, mudstone,

 and conglomerate

Limestone

Granite and metamorphic rock

Boundary of alluvium

0–100

100–200

200–300

300–400

400–500

EXPLANATION

Figure 2. Geology and saturated thickness of the regional aquifer, Upper San Pedro basin, United States and Mexico.

P

de r

River

o

Bisbee

NORTHERN BOUNDARY

OF GROUND-WATER

FLOW MODEL

 Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Figure 3. Generalized hydrogeologic section and extent of numerical model layers, Upper San Pedro Basin, United States and Mexico.

Ground water is transmitted primarily through layers of 

permeable sand and gravel within the basin fill and Quaternary 

terrace and alluvial deposits. Poorly permeable silt and 

clay layers are interbedded with the sand and gravel layers 

and significantly influence ground-water flow. The silt and 

clay layers result in confined ground-water flow within the 

underlying permeable sand and gravel layers; and limit the 

ability of recharge water to infiltrate downward to the aquifer 

where the silt and clay layers lie above the water table.

Pre-Tertiary Rocks

Outcrops of pre-Tertiary rocks in the study area are 

predominantly limestone and clastic rocks (fig. 2) with local 

areas of crystalline rocks. Ground-water flow through these 

rocks may be an important part of the overall ground-water 

flow system, but permeability of the rocks is not well known. 

Much of the limestone is likely permeable because springs 

are common near the contact with low-permeability rock. 

Clastic sediments include sandstone, siltstone, and mudstone. 

Sandstone likely has variable permeability, and the siltstone 

and mudstone have low permeability. Outcrops of lowpermeability crystalline rocks are barriers to ground-water 

flow in pre-Tertiary sedimentary rocks.

Limestone occurs in the Huachuca Mountains, Whetstone 

Mountains, Mustang Mountains, Dragoon Mountains, Mule 

Mountains, and Sierra Los Ajos. Limestone likely underlies 

the alluvial-basin deposits in much of the region south of the 

Tombstone Hills and between the Huachuca Mountains and 

Mule Mountains. 

Clastic sedimentary rocks are common in the Mustang 

Mountains, along the western and southern flanks of the 

Huachuca Mountains, in the southern Mule Mountains, and 

throughout much of the portion of the basin in Mexico. Clastic 

sedimentary rocks underlie much of the alluvial-basin deposits 

in the area between the Tombstone Hills and the Mustang 

Mountains and in Mexico.

Granite and metamorphic crystalline rocks occur on the 

eastern flank of the Huachuca Mountains, Mule Mountains, 

Sierra Los Ajos, and the Tombstone Hills. Two types of 

crystalline rocks occur in the Tombstone Hills—granidiorite 

exists in a small region on the northern edge of the Tombstone 

Hills, and volcanic rocks exists, in the western part of the 

Tombstone Hills. The crystalline rock has low permeability. 

Little is known about the permeability of the volcanic 

rocks; however, some volcanic rocks may have moderate 

permeability. In some areas, such as the southwestern flank of 

the Huachuca Mountains, volcanic rocks are interbedded with 

clastic sediments, and ground water occurs in caves and issues 

from springs.

Tertiary Pre-Basin Fill Sediments

The Pantano Formation is a conglomerate that was 

deposited in an alluvial basin during low-angle tectonic 

extension of the region and before the Basin and Range 

structural disturbance. The unit was tilted and faulted by later 

tectonism, resulting in a southwestern dip of 15–45 degrees 

(Brown and others, 1966). The unit is generally of low 

permeability because of cementation (Brown and others, 

1966), but can yield water to wells through fractures and 

may be an important local water-bearing unit. Gravity studies 

(Gettings and Houser, 1995; Halverson, 1984) indicate the 

Pantano Formation is likely several hundred meters thick 

in two structural depressions in the west-central part of the 

Sierra Vista subwatershed. The two depressions are separated 

by an east-west trending ridge in the subsurface near Sierra 

Vista (Gettings and Houser, 1995). 

Conceptual Model of the Ground-Water Flow System 

Figure 3. Generalized hydrogeologic section and extent of numerical model layers, Upper San Pedro basin, United States and Mexico. 

PRE- AND POST-ENTRENCHMENT

 STREAM ALLUVIUM

BASIN FILL SAND AND GRAVEL

BASIN FILL SILT AND CLAY

PRE-BASIN FILL ALLUVIAL DEPOSITS

SEDIMENTARY ROCK

CRYSTALLINE ROCK

MODEL LAYER 1

MODEL LAYER 2

MODEL LAYER 3

MODEL LAYER 4

MODEL LAYER 5

Water table

EXPLANATION WEST

EAST

Lower basin fill

basin fill

Upper
Basin Fill

The basin fill in the San Pedro River Basin is similar 

in lithology and structure to basin fill throughout southern 

Arizona (Anderson and others, 1992). The unit was deposited 

during and after Basin and Range structural deformation that 

created the current distribution of basins. Coarse-grained 

and fine-grained facies of the alluvial sediments generally 

occur at the basin margins and basin centers, respectively. 

Vertically, the sediments generally grade from coarsegrained sediments of sand and gravel in the lowest part to 

fine-grained sediments of silt and clay in the middle part and 

interbedded coarse-grained and fine-grained sediments in the 

upper part. The sediments that accumulated in the basins are 

generally divided into upper and lower hydrostratigraphic 

units. The lower basin fill tends to be more massively bedded 

and indurated and includes coarse-grained parts commonly 

described as conglomerate and fine-grained parts commonly 

described as siltstone and mudstone. The upper basin fill is 

more heterogeneous and poorly indurated. The basin fill in 

the study area differs only in detail from basin fill in nearby 

basins. The upper and lower basin fill units are well defined 

only in the northern part of the Sierra Vista subwatershed 

for which good subsurface information from lithologic and 

electric logs is available (Pool and Coes, 1999). The top and 

bottom of the basin fill units and the contact between the units 

in the southern part of the subwatershed were defined for the 

ground-water flow model on the basis of a few lithologic and 

electric logs and electrical geophysical surveys. The basin fill 

in Mexico is undifferentiated because lithologic information is 

insufficient to distinguish between the two basin fill units.

The lower basin fill is an important water-bearing unit 

throughout most of the basin. The unit unconformably overlies 

the Pantano Formation and older rocks, consists of interbedded 

sand and gravel of variable cementation on the basin margins, 

and includes thick facies of siltstone and mudstone (Brown 

and others, 1966; table 1). Information from several test wells 

indicates that the thickness of the lower basin fill ranges from 

about 45 m to 100 m near Fort Huachuca. The unit is much 

more thick, however, in the two major structural depressions.

The upper basin fill lies above a depth of 120 m in all 

wells where it has been recognized. Saturated thickness is 

generally greatest near the basin center and thinnest near 

areas of bedrock outcrop and ranges from 0 to 57 m at test 

wells on Fort Huachuca where it has been identified. Other 

lithologic and geophysical data indicate that the thickness 

of the unit generally is less than 100 m, but is as much as 

300 m thick near Hereford and Palominas (fig. 1). The unit is 

conformable with the lower basin fill and consists of weakly 

cemented and compacted soft reddish brown clay, gravel, 

sand, and silt (Brown and others, 1966; table 1). The upper 

basin fill includes a relatively permeable fan gravel facies 

near the mountains that grades laterally to a poorly permeable 

silt and clay facies with interspersed sand and caliche beds 

near the basin center. The unit is primarily a confining bed 

of silt and clay where it is saturated between Sierra Vista and 

the San Pedro River and between Hereford and Highway 

90 along the San Pedro River (fig. 1). The upper basin 

fill is an important aquifer where the fan gravel facies is 

saturated. Sand beds within the silt and clay facies also may 

transmit substantial amounts of water provided individual 

beds are sufficiently interconnected; however, the extent of 

individual beds is not well known because of a lack of detailed 

subsurface information. The upper basin fill is equivalent 

to the St. David Formation described by Gray (1965), 

which crops out extensively north of the study area near 

St. David, Ariz.. 

Distribution of silt and clay in the basin fill is defined 

by subsurface data from lithologic logs, vertical electrical 

soundings (Pool and Coes, 1999; Consultores en Agua 

Subterranea S.A. Por Mexicana de Cananea, S.A. de C.V., 

2000; Fleming and Pool, 2002), and other surface and aerial 

geophysical surveys (University of Arizona Geophysics Field 

Camp 2001, 2002, 2004; Condor Consulting, 2003; fig. 4). 

The vertical and lateral extent of significant silt and clay 

intervals was estimated as the extent of the basin fill where the 

aerial electromagnetic survey mapped an electrical resistivity 

of 12 ohm-m or less (fig. 4). The data were discretized at 10 m 

intervals; therefore, intervals of silt and clay of less than about 

10 m thickness were excluded from the results. The base of 

the silt and clay is uncertain in areas where the basin fill is 

underlain by electrically conductive sedimentary rocks, such 

as siltstone and mudstone, in much of the southern half of the 

study area and near the Babocomari River. The extent of the 

silt and clay defines the extent of confined parts of the aquifer. 

The silt and clay intervals are generally 10–300 m thick. Local 

unconfined aquifers occur where as much as 100 m of sand 

and gravel facies overlie silt and clay in the south-central part 

of the Sierra Vista subwatershed. 

Stream Alluvium

Stream alluvium along the San Pedro and Babocomari 

Rivers is a locally important water-bearing unit that forms a 

local unconfined aquifer. The unit unconformably overlies 

lower basin fill and volcanic rocks in the Charleston, Ariz., 

area and overlies upper basin fill above and below Charleston 

and along the Babocomari River. The oldest deposits of stream 

alluvium are clay, silt, and fine sand, with interbedded coarse 

sand and pebble to cobble gravel that were deposited before 

the river was entrenched in about 1890 (Hereford, 1993). 

The pre-entrenchment deposits are as much as 6 m thick 

and 1.5 km wide. Sand and gravel in the post-entrenchment 

alluvium were deposited after stream entrenchment of 

1–10 m in a narrow channel within the floodplain and preentrenchment alluvium. The younger deposits are a few 

meters thick. The post-entrenchment alluvium is permeable. 

Pre-entrenchment alluvium transmits water, but generally is 

less permeable than overlying and adjacent post-entrenchment 

alluvium and underlying deposits of basin fill.

 Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Ground-Water Flow System

Water recharges the ground-water flow system of 

the Upper San Pedro Basin in the higher elevations of the 

basin and discharges as spring flow and stream flow, and 

through riparian vegetation, in regions at lower elevations 

and as underflow to the downgradient basin. No ground 

water is assumed to enter the basin from adjacent basins. 

The distribution of ground-water flow and locations of 

discharge areas (fig. 5) are influenced by locations and 

rates of recharge, basin topography, and the distribution of 

hydraulic properties in the subsurface. The flow system has 

been modified by changes in discharge, including changes 

in ground-water withdrawals by wells and vegetation, and 

changes in recharge rates due to changes in climate and 

land use. 

Recharge occurs through deep percolation of 

precipitation that is in excess of evaporation, transpiration by 

vegetation, and runoff. Recharge is concentrated in two areas: 

(1) at high elevations and (2) in ephemeral stream channels, 

where permeable near-surface materials, soil or rock, facilitate 

infiltration. At high elevations, precipitation rates are greatest, 

and evaporation and transpiration demands are low, allowing 

high potential rates of infiltration through permeable rocks, 

such as limestone or sandstone. Surface runoff of precipitation 

in excess of infiltration is concentrated in ephemeral channels 

that commonly have highly permeable deposits that readily 

transmit excess infiltrated water below depths accessible by 

vegetation. The resulting recharge distribution is concentrated 

in areas of permeable rocks in the mountains and along stream 

channels underlain by sand and gravel and the alluvial aquifer. 

The greatest rates of infiltration are along ephemeral stream 

channels near the mountains and are commonly referred to as 

mountain-front recharge.

Ground water flows through the permeable rocks and 

sediments toward the low-elevation discharge areas along the 

major streams (fig. 5). Ground water discharges naturally to 

streams, springs, and riparian vegetation. Springs in incised 

canyons where erosion has exposed an underlying rock of 

low permeability are common discharge areas for water 

that flows in the limestone and sandstone aquifers in the 

mountains. These types of springs are common at the lower 

elevations of the Huachuca Mountains, Mule Mountains, 

and Sierra San Jose. Not all of the water flowing through the 

sedimentary-rock aquifers discharges at high-elevation springs. 

Some of the water flows downgradient into the alluvial aquifer 

at the southern end of the Huachuca Mountains, below much 

of the Mule Mountains, and throughout much of the alluvial 

basin in Mexico. Common areas of discharge from the alluvial 

aquifer occur where stream channels incise the top of the 

saturated zone (water table). Several tributary channels have 

incised the water table where it closely overlies thick intervals 

of silt and clay, creating springs and perennial stream reaches. 

Discharge also occurs through evapotranspiration near the 

springs and along perennial-stream reaches.

Ground water in the Sierra Vista subwatershed also 

discharges to the adjacent ground-water flow system in the 

Benson subwatershed as underflow through aquifers that are 

common to both basins (fig. 5). Most of the ground water 

flowing to the Benson subwatershed flows through stream 

alluvium near the streamflow-gaging station near Tombstone 

and through basin fill west of the Dragoon Mountains and east 

of outcrops of crystalline rock adjacent to the streamflowgaging station (fig. 5). Much of the recharge within the Walnut 

Gulch watershed likely discharges to the Benson subwatershed 

ground-water system through basin fill in this area. Ground 

water flowing to the Benson subwatershed also may flow 

through basin fill between the Tombstone streamflow-gaging 

station and outcrops of crystalline rock about 5 km northwest 

of the streamflow-gaging station. Little ground water probably 

flows northward across the boundary between the two 

basins beneath the broad alluvial plain east of the Whetstone 

Mountains. Water-level elevations indicate that flow is parallel 

to the boundary in this area (fig. 5).

The rate of ground-water flow between the Sierra Vista 

and Benson subwatersheds likely varies with changes in 

recharge rates and ground-water withdrawals in both basins. 

Hydrographs of wells in the basin fill near the boundary, 

however, indicate that water levels changed little in the 

area (Barnes and Putman, 2004), and that variations in 

ground-water discharge through the basin fill to the Benson 

subwatershed are small.

Ground-Water Budget

Ground-water budgets were constructed for both 

steady-state and transient conditions. Estimates of the 

steady-state budget determine the average rate of natural 

flow through the aquifer system and include the naturally 

occurring components of recharge and discharge through 

ground-water flow to adjacent aquifers, streams, drains, and 

evapotranspiration. The transient budget includes variations in 

ground-water withdrawals, evapotranspiration, and artificial 

and incidental recharge. The steady-state ground-water budget 

can be estimated by using water-balance methods because 

most recharge to the aquifer system discharges as stream 

baseflow. Baseflow from about 70 percent of the watershed 

has been measured continuously since the early 1930s at the 

streamflow-gaging station at Charleston, Ariz. Steady-state 

discharge through evapotranspiration is uncertain, but was 

estimated for the riparian area along the San Pedro River 

within the United States for 2002 and 2003 (Leenhouts and 

others, 2005). The transient budget is estimated from available 

data on ground-water withdrawals, irrigated area, riparian 

vegetation, and artificial recharge.

Conceptual Model of the Ground-Water Flow System 
Figure 4. Extent and altitude of thick silt and clay interval in the basin fill, Upper San Pedro Basin, United States and Mexico. 

4A, Top of thick silt and clay interval; 4B, Base of thick silt and clay interval.

30’

31o45’

0 10 KILOMETERS

0 6 MILES

PIMA COUNTY

SANTA CRUZ

COUNTY

ARIZONA

SONORA

UNITED STATES

MEXICO

Sierra Vista

Cananea

Tombstone

B ari

abocom

R vi

R o

nos oL s e s r F

i

er

Palominas

Huachuca City

Ft. Huachuca

COCHISE

COUNTY

110o30’

15’

31o

110o

15’

Hereford

San Pedro

Riparian

National

Conservation

Area

Charleston

R vi e r

P

de ro

San

Elgin

Fairbank

Figure 4. Extent and altitude of: A, The top of thick silt and clay interval in the basin fill; B, The base of the thick silt and clay interval in the

basin fill. 

A. The top of the thick silt and clay interval

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

1,400–1,500

1,300–1,400

1,200–1,300

1,100–1,200

1,000–1,100

ALTITUDE OF TOP OF

SILT AND CLAY,

IN METERS ABOVE

NAVD 88:

EXPLANATION

Bisbee

10 Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

PIMA COUNTY

SANTA CRUZ

COUNTY

ARIZONA

SONORA

UNITED STATES

MEXICO

Tombstone

B ari

abocom

R vi

R o

nos oL s e s r F

i

er Huachuca City

Ft. Huachuca

COCHISE

COUNTY

110o30’

15’

31o

110o

15’

San Pedro

Riparian

National

Conservation

Area

R vi e r

P

de ro

San

Figure 4. Continued.

B. The base of the thick silt and clay interval

Sierra Vista

1,300–1,400

1,200–1,300

1,100–1,200

1,000–1,100

ALTITUDE OF BASE OF

SILT AND CLAY,

IN METERS ABOVE

NAVD 88:

EXPLANATION

Bisbee

Figure 4. Continued.

Conceptual Model of the Ground-Water Flow System 11
30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

UNITED STATES

MEXICO

Sierra Vista

1,

1,160

280

1,3 0

006 1,

003 1,

1, 1,540 560

1,480

1,36 0

1,320

1,2 20

1,2 6 0

1,480

1,420

1,400

1,420

1,540

1,520

1,440

1,440

1,380

1,420 1,460

1,440

1,460

1,400

1,380

1,340

1,300

1,560

1,

1,580

260

1,340

1,400 1,3 80

1,400

1, 4 20

1,280

1,380

1,460

1,500

1,460

1,360 1,320

1,280

1,240

1,480

1,500

1,520

1,520

1,460

1,440

1,420

1,420

1,520

110o30’

15’

15’

31o

110o

P

San

de ro

River

No-flow boundary

based on estimated

ground-water flow path

No-flow boundary based on

estimated ground-water flow path

MODEL GRID ROW 1

MODEL GRID COLUMN 1

MODEL GRID COLUMN 320

MODEL GRID ROW 440 Ground-water flow model boundary

Cananea

NOTE: Water levels in Mexico are based on

 post-1960 data (Consultores en Agua

 Subterranea S.A. for Mexicana de

 Cananea, S.A. de C.V., 2000)

BOUNDARY CONDITIONS

Evapotranspiration area

No flow

Specified head

PREDEVELOPMENT WATER-LEVEL

 ALTITUDE—Contour interval

 20 meters. NAVD 88.

WELL —Pre-1961 water level

EXPLANATION

1,300

Spring simulated as a drain

Streamflow-routing network

Spring simulated as a stream

WELL WITH HYDROGRAPH

STREAMFLOW-GAGING

 STATION WITH HYDROGRAPH

(D–19–21)12dbb2

(D–20–20)02ddd

1, 580

00 1, 5

0

Tombstone

1, 4 60

1,600

Figure 5. Ground-water flow system, boundary conditions for the numerical model, streamflow-routing network, evapotranspiration areas,

and predevelopment water-level comparison sites. Upper San Pedro basin, United States and Mexico. 1,200

1,180

Palominas

Cottonwood

Lewis Springs

Bisbee

Figure 5. Ground-water flow system, boundary conditions for the model, streamflow-routing network, evapotranspiration areas, and 

predevelopment water-level comparison sites, Upper San Pedro Basin, United States and Mexico.

12  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Predevelopment

A predevelopment steady-state period is difficult to 

define in the Upper San Pedro Basin because of documented 

stream-channel incision, observed variations in stream 

baseflow, estimated variations in evapotranspiration rates, and 

the uncertain transient effects of early withdrawals at mines 

for dewatering purposes. Stream-channel incision prior to 

1900 and subsequent widening of the stream channel through 

the mid-1950s (Hereford, 1993) likely induced ground-water 

level decline and increased rates of baseflow discharge from 

the ground-water system for an undetermined period. Previous 

studies of the aquifer system have assumed that steady-state 

conditions existed during 1940, prior to extensive development 

of the ground-water supply (Freethey, 1982; Corell and others, 

1996; Goode and Maddock, 2000). Variations in baseflow at 

Charleston since the mid-1930s, however, indicate that rates 

of discharge from the aquifer system changed throughout the 

period in response to changing rates of recharge and discharge 

(Pool and Coes, 1999). One likely cause of variations in 

discharge rates is climate-induced variations in recharge that 

occur on decadal (PDO) and more frequent (ENSO/El Niño) 

scales (Hanson and others, 2004; Dickinson and others, 2004; 

Pool, 2005). The steady-state period assumed by previous 

investigations to occur in 1940 was a wet year at the end of a 

decadal wet period. Mine dewatering at Tombstone and near 

Bisbee prior to 1940 also caused changes to the ground-water 

flow system; however, the effects were largely unmonitored 

and uncertain. As a result of mine dewatering before 1940, 

hydrologic information for 1940 may not represent steadystate conditions. The predevelopment ground-water budget 

developed for the ground-water flow model described in this 

report includes the earliest data available as well as data from 

the streamflow-gaging station at Charleston during 1935–39 

and estimates of evapotranspiration rates during 1935. 

Estimates of predevelopment discharge rates are assumed to 

be balanced by an equivalent recharge rate.

Discharge as ground-water underflow was estimated for 

three locations along the northern extent of the model area; the 

combined rate was about 4,000 m3/d (table 2). Ground-water 

discharge through basin fill east of the San Pedro River was 

estimated to be equivalent to recharge in the Walnut Gulch 

Watershed, about 2,700 m3/d. Ground-water discharge from 

the Sierra Vista subwatershed through stream alluvium at 

the streamflow-gaging station near Tombstone is 1,000 m3/d 

on the basis of previous estimates (Corell and others, 1996). 

Spring discharge measurements of 300 m3/d (John Sottilare, 

Hydrologic Technician, U.S. Bureau of Land Management, 

written commun., 2005) at the northern extent of the modeled 

area were used to estimate minimum rates of ground-water 

discharge through basin fill west of the San Pedro River. 

Additional unaccounted ground water may flow through basin 

fill south of the springs. 

Natural discharge of ground water occurs as baseflow to 

the San Pedro and Babocomari Rivers, at several springs near 

the San Pedro River, and at several springs in the Huachuca 

Mountains (fig. 5). San Pedro River baseflow at the Charleston 

and Palominas streamflow-gaging stations (fig. 5) varied 

seasonally with evapotranspiration demands and near-stream 

agricultural withdrawals. Stream baseflow at the Charleston 

streamflow-gaging station varied seasonally during 1936–40 

from about 35,750 m3/d during the winter to about 8,661 m3/d 

during June and averaged about 19,900 m3/d annually. 

Stream baseflow to the Babocomari River is not well 

documented but is estimated to be about 8,300 m3/d on the 

basis of an analysis by Schwartzman (1990). Discharge to 

springs that issue from the basin fill near Lewis Springs and 

springs in the Huachuca Mountains (fig. 5) was much less 

than that in the other discharge areas and averaged about 

1,500 m3/d (table 2). 

The predevelopment rate of ground-water discharge is 

poorly defined. Previous investigations of the Upper San Pedro 

basin in the United States and Mexico have estimated about 

26,400 m3/d of ground-water discharge through riparian 

evapotranspiration (Freethey, 1982; Esparza, 2002; table 2). 

Previous investigations of changes in the riparian area in the 

San Pedro Valley (Reichardt and others, 1978) and recent 

comparisons of repeat aerial photographs in the Sierra Vista 

subwatershed (Russ Scott, Research Hydrologist, Agriculture 

Research Service, written commun., 2005) indicate that 

predevelopment evapotranspiration rates may be as little as 

40 percent of post-1970 rates.

Total annual discharge from the predevelopment system 

was about 60,600 m3/d, which included estimated discharge 

components of 4,000 m3/d for ground-water discharge, 

28,700 m3/d for baseflow discharge to streams and springs, 

and 26,400 m3/d for evapotranspiration through riparian 

vegetation (table 2). The annual rate of recharge to the groundwater flow system was equivalent to the annual rate of groundwater discharge, 60,600 m3

/d, provided that steady-state 

conditions existed during predevelopment.

Table 2. Estimated predevelopment ground-water budget. 

 [cubic-meters per day]

Water-budget component

Average annual

Inflow Outflow

Recharge1 60,600 0

Stream baseflow2 0 28,700

Evapotranspiration 0 26,400

Ground-water underflow 0 4,000

Drains3 0 1,500

Storage4 0 0

TOTAL 60,600 60,600

1

Recharge is estimated as equivalent to the estimated outflow components.

2

Stream baseflow represents net baseflow because predevelopment gaining and losing reaches are poorly defined.

3

Drains include discharge at several springs in the Huachuca Mountains.

4

Storage inflow represents removal from storage and input to the groundwater flow system. Storage outflow represents removal from the groundwater flow system and input to storage. No storage change is assumed to 

occur during predevelopment.

Conceptual Model of the Ground-Water Flow System 13
Postdevelopment

The postdevelopment ground-water budget includes 

variations in natural and artificial recharge rates, 

evapotranspiration, and ground-water withdrawals. Variations 

in natural recharge (Hanson and others, 2004; Dickinson 

and others, 2004; Pool, 2005) were not considered in this 

analysis, but may contribute significantly to water-level 

and stream-baseflow trends. Variations in artificial recharge 

were considered and result from changes in the amounts of 

infiltration of excess irrigation water, sewage effluent near 

municipalities, waste water at individual domestic septic 

systems, and water discharged from mines for dewatering 

purposes. Variations in evapotranspiration rates are not 

well known for most of the postdevelopment period, but 

variations were estimated on the basis of vegetative mapping 

that was done at different times. Variations in ground-water 

withdrawals were estimated on the basis of both monitored 

and estimated withdrawals. Except for documented 

withdrawals at mines and at Fort Huachuca, ground-water 

withdrawals prior to about 1970 were estimated by using data 

from previous publications. After about 1970, withdrawals 

were documented for the mines and for Fort Huachuca, and 

deliveries of ground-water withdrawals were documented 

for municipal water systems and private water companies. 

Withdrawals for other uses, including irrigation and private 

use, were estimated on the basis of previous investigations. 

Variations in all water-budget components, except natural 

recharge, were estimated for 1902–2002.

Artificial Recharge

Artificial recharge in the study area occurs as the result 

of infiltration of excess irrigation water; sewage effluent 

near Sierra Vista, Fort Huachuca, Tombstone, Bisbee, and 

Naco; waste water at individual domestic septic systems; and 

water discharged for dewatering purposes from mines near 

Bisbee and Tombstone. Incidental recharge also occurs in the 

distribution of potable water from supply wells to individual 

users. Water providers report deliveries to individual users 

rather than well withdrawals, therefore, incidental recharge 

from the distribution systems is implicitly included by 

reported deliveries.

Recharge of excess irrigation water occurred primarily 

in the agricultural areas near the San Pedro River near 

Palominas and in Sonora, Mexico (fig. 1). Recharge occurred 

also at Warren Ranch near Bisbee beneath agricultural fields 

that were irrigated with water pumped from mine workings 

during 1905–87. Recharge of irrigation water occurred also 

beneath golf courses. Rates of recharge from excess irrigation 

depend on irrigation practices and efficiencies. Previous 

ground-water flow models of the area assumed that irrigation 

efficiencies were about 70 percent, which resulted in 30 

percent of the irrigation withdrawals returning to the aquifer 

through deep percolation (Freethey, 1982; Corell and others, 

1996; Goode and Maddock, 2000). Deep percolation beneath 

golf courses near Naco, Sierra Vista, and Fort Huachuca was 

estimated as 1.4 x 10-3 m/d to the turf area (Phelps Dodge 

Corporation, 1998).

Recharge from deep percolation of sewage effluent 

occurred at individual septic systems and below sewage 

treatment facilities at Sierra Vista, Naco, Bisbee, Fort 

Huachuca, Huachuca City, and Tombstone (fig. 1). Effluent 

treated prior to 1970 was typically discharged into ephemeral 

channels. Improved facilities often used treated effluent 

to irrigate turf and crops, or discharged it to evaporation 

or infiltration ponds. Since early treatment records were 

unavailable, the average of recent estimates of multiple 

municipal effluent-recharge rates in the Sierra Vista 

subwatershed (Arizona Department of Water Resources, 

2005), which was 14 percent of total delivered withdrawals, 

was used to estimate recharge rates. The City of Sierra Vista 

began treating sewage effluent in 1967, increased the capacity 

with an additional facility in 1978, and converted the treatment 

facility to a recharge facility in 2001. During 1967–79, 

effluent was discharged into a nearby ephemeral channel. 

Beginning in 1980, effluent was disposed of through irrigation 

of 3.4 x 106 m2 near the facility. Beginning in July 2002, 

effluent was discharged to recharge basins and recharged at a 

rate of 2,960 m3/d (Arizona Department of Water Resources, 

2005). Effluent was discharged to ponds near Naco and 

Bisbee beginning in about 1930 (Phelps Dodge Corporation, 

1998); however, the Bisbee effluent ponds were outside of the 

watershed. Effluent was recharged at Fort Huachuca during 

1978–95 at two turf sites, a golf course and parade grounds, 

and below the treatment facility at rates of about 1,360 m3/d 

(United States Department of Defense, 2002). During 2001, 

about 1,825 m3/d, or approximately 33 percent of withdrawals, 

was estimated to have been recharged. Treatment facilities at 

Huachuca City discharge to evaporation ponds where little 

recharge occurs (Arizona Department of Water Resources, 

2005). Treatment facilities at Tombstone began discharging 

effluent into Walnut Gulch in 1972. Recharge at individual 

septic systems was estimated to be 14 percent of withdrawals 

on the basis of estimates made for a decision support system 

developed for the Upper San Pedro Partnership (Kevin 

Lansey, professor, University of Arizona Department of 

Civil Engineering and Engineering Mechanics, written 

commun., 2005). 

Evapotranspiration

Evapotranspiration from riparian vegetation near 

the San Pedro and Babocomari Rivers was defined by 

investigations that mapped historical extents of riparian 

vegetation and recent extents and evapotranspiration rates 

for several riparian-vegetation types. Unpublished maps of 

riparian-vegetation extent in 1935 were available from the 

Arizona Department of Water Resources. Several studies of 

riparian vegetation were used to estimate change, including 

those by Reichardt and others (1978), U.S. Army Corps 

of Engineers (2001), Kepner and Edmonds (2002), and 

Leenhouts and others (2005). 

14  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
The extent of riparian vegetation changed during the 

20th century. Reichardt and others (1978) evaluated totalacreage changes for dense riparian vegetation and other land 

uses by using aerial photographs taken in 1935, 1966, and 

1977. Total dense riparian land along the San Pedro River 

in the United States increased by 79 percent from 1935 to 

1966 and increased by more than 100 percent from 1935 to 

1978. Comparable results were derived from comparing 2001 

maps (U.S. Army Corps of Engineers, 2001) with 1935 aerial 

photographs, which indicated the extent of riparian vegetation 

had increased by 150 percent within the Sierra Vista 

subwatershed (Russ Scott, research hydrologist, Agriculture 

Research Service, written commun., 2005). Results of 

vegetation-type mapping across the Upper San Pedro Basin 

from Landsat satellite photos taken in 1973, 1986, 1992, 

and 1997 (Kepner and Edmonds, 2002) indicate mesquite 

landcover increased significantly and riparian landcover 

decreased about 20 percent during 1973–86 but recovered to 

106 percent of the 1973 coverage during 1992–97. 

Changes in extent of riparian vegetation imply variations 

in rates of evapotranspiration. Estimates of changes in rates 

of riparian evapotranspiration, however, require information 

on changes in the extent and density of several vegetation 

types that use ground water from different depths and at 

different rates. Major vegetation types within the riparian 

area that use ground water at different rates include mesquite, 

cottonwood and willow, and sacaton grass. Detailed maps of 

each vegetation type are not available for the 1935 period; 

therefore, estimates of change in riparian evapotranspiration 

can only be derived from available information that documents 

change in the extent of riparian vegetation. The most recent 

and accurate estimates of evapotranspiration rates are 

available from Leenhouts and others (2005). Estimated rates 

of evapotranspiration within the Sierra Vista subwatershed 

range from 32,400 to 40,700 m3/d during 2003, but rates vary 

annually by as much as 30 percent on the basis of 3 years of 

data. An estimate of the rate of evapotranspiration in the part 

of the basin in Mexico is 4,500 to 5,700 m3/d, assuming rates 

of riparian evapotranspiration per unit area in the Mexico 

part of the basin are the same as rates in the Sierra Vista 

subwatershed. Total estimated evapotranspiration in the model 

area during 2003 ranges from 37,000 to 46,000 m3/d, which 

is 50–90 percent greater than the predevelopment estimate 

of 26,400 m3/d. The percent difference between estimated 

predevelopment and 2003 rates of evapotranspiration is 

similar to the estimated percent change in extent of riparian 

vegetation; both values suggest that recent evapotranspiration 

rates within the model area may be about twice the 

evapotranspiration rates during 1935.

Ground-Water Withdrawals

Demand for ground water in the study area is divided 

into five main categories: (1) mine withdrawals from mines 

near Tombstone, Bisbee, and Cananea, which also includes 

withdrawals for municipal and domestic use; (2) agricultural 

withdrawals for irrigated areas primarily near the San Pedro 

River; (3) municipal and industrial withdrawals by water 

purveyors for deliveries within incorporated cities and some 

adjacent unincorporated areas; (4) Fort Huachuca withdrawals 

for military use; and (5) withdrawals in the Sierra Vista 

subwatershed from domestic wells in unincorporated areas 

and stock wells (figs. 1 and 6). Withdrawals prior to about 

1940 were dominated by mine dewatering near Tombstone 

and Bisbee. Large volumes of ground water withdrawn for 

other purposes became possible with the widespread use of 

high-powered pumps around 1940. Subsequent increases in 

agriculture, municipal, military, and industrial water demands 

have been met by drilling new wells and installing well pumps 

with larger capacities. Total ground-water withdrawals in 

the study area varied from about 38 to 55 hm3/yr (30,800 to 

46,000 acre-ft/yr) during 1979–2002 and peaked in 1988 

(fig. 6).

Ground-water withdrawals were estimated by using 

historical data from previous studies. Withdrawals in the 

Sierra Vista subwatershed generally were estimated by using 

indirect estimates for agricultural use and reported population 

to indicate pumping volumes for domestic withdrawals in 

unincorporated areas. Reported- and metered-delivery records 

from water suppliers were used when available. Records of 

withdrawals at the Tombstone and Copper Queen mines were 

available from Hollyday (1963) and Southwest Ground-water 

Consultants (2004), respectively. Annual withdrawals for all 

uses were estimated in previous studies for the period 

1940–97. Withdrawals for different periods were apportioned 

for this study by using information from several previous 

studies, including those for 1940–77 by Freethey (1982), 

1978–88 by the Arizona Department of Water Resources 

(1990), 1989–90 by Corell and others (1996), 1991–97 by 

Goode and Maddock (2000), and 1998–2002 by Arizona 

Department of Water Resources (2005). Withdrawals for 

use at the Cananea mine and domestic and agricultural uses 

in Mexico were reported by Esparza (2002). Withdrawals 

for the Cananea mine were used, in part, for the City of 

Cananea. Reported and metered withdrawals in the Sierra 

Vista subwatershed for 1998–2002 were estimated by 

deliveries of public-water purveyors reported to the Arizona 

Corporation Commission. Fort Huachuca withdrawals have 

been documented since 1940 (Mike Shaughnessey, Directorate 

of Public Works Real Property Office, U.S Army Garrison, 

Fort Huachuca, Arizona, written commun., 2003). Data for 

domestic use in unincorporated areas of the Sierra Vista 

subwatershed was provided by the Arizona Department of 

Water Resources for 1990 and 2002 (Arizona Department of 

Water Resources, 2005).

Conceptual Model of the Ground-Water Flow System 15
Figure 6. Annual ground-water withdrawal rates for selected uses, 1902 to 2002, Upper San Pedro Basin, United States and Mexico. 

6A, Mine uses; 6B, Non-mine uses; 6C, Total uses.

16  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico

Figure 6. Annual ground-water withdrawal rates for selected uses, 1902 to 2002, Upper San Pedro basin, United States and Mexico: A, 

Mine uses; B, Non-mine uses; C, Total uses.

B. Non-mine uses GROUND-WATER WITHDRAWALS, IN CUBIC HECTOMETERS

GROUND-WATER WITHDRAWALS,

IN ACRE-FEET

GROUND-WATER WITHDRAWALS,

IN ACRE-FEET

40

35

30

25

20

15

10

5

30,000

20,000

25,000

15,000

10,000

5,000

0

10,000

20,000

0

GROUND-WATER WITHDRAWALS,

IN ACRE-FEET

20,000

15,000

10,000

5,000

30,000

40,000

45,000

0

A. Mine uses GROUND-WATER WITHDRAWALS, IN CUBIC HECTOMETERS

Tombstone mine

Copper Queen mine in Bisbee

Cananea mine

Total mine withdrawals

30

25

20

15

10

5

0

C. Total uses GROUND-WATER WITHDRAWALS, IN CUBIC HECTOMETERS

60

50

40

30

20

10

0

1900190519251930

19101915

19201935194019451950195519601965197019751980198519901995

20002005

0

1900190519251930

19101915

19201935194019451950195519601965197019751980198519901995

20002005

1900190519251930

19101915

19201935194019451950195519601965197019751980198519901995

20002005

U.S. agriculture including vineyards

U.S. incorporated

U.S. unincorporated

Fort Huachuca

Mexico domestic and agriculture

Total non-mine withdrawals

NOTE: The portions of Copper Queen mine withdrawals from the

 Upper San Pedro and adjacent basin are uncertain
Withdrawals for dewatering of mines and other mining 

uses exceeded all other withdrawals in the study area during 

1902–60 (fig. 6). The earliest mining withdrawals began 

with pumping 44 hm3 (35,700 acre-ft) of water from the 

underground mine workings at Tombstone during 1902–11. 

Additional small amounts of withdrawals occurred from 

the mines at Tombstone during later periods. Water was 

withdrawn from the Copper Queen mine at Bisbee to 

dewater the mine workings and for leaching operations 

during 1906–87 at rates that ranged from 1 hm3/yr to more 

than 62 hm3/yr (800–16,400 acre-ft/yr). The withdrawals, 

however, were likely derived from both the Upper San Pedro 

and adjacent ground-water flow systems because the 

withdrawals straddle the surface-water divide between the 

basins. The part of the withdrawals that were derived from 

each ground-water basin is uncertain. Furthermore, lowpermeability faults may isolate the water in rocks of the 

mine from the regional aquifer (Southwest Ground-water 

Consultants, 2004), therefore, withdrawal of water from 

storage at the mine should have little influence on the regional 

aquifer in the Upper San Pedro Basin.

Ground-water withdrawals from the Cananea mine and 

City of Cananea began in 1940 and increased to 7.8 hm3/yr 

(6,300 acre-ft/yr) in 1985 (fig. 6). Withdrawal rates increased 

significantly during 1986–2002 and ranged from 12 to 

27 hm3/yr (9,700-21,800 acre-ft/yr, Esparza, 2002). Ground 

water was withdrawn near Cananea from many wells that tap 

the basin fill throughout a large area between Sierra San Jose 

and Sierra Mariquita (fig. 1).

Agricultural withdrawals were determined by using 

previous estimates, crop consumptive-use rate, irrigated area 

determined from aerial photographs, remote sensing, and 

agricultural surveys in the United States portion of the study 

area. The most accurate data are available from field surveys 

conducted in the United States in the late 1980s and in 2002 

by the Arizona Department of Water Resources (Arizona 

Department of Water Resources, 1990 and 2005). Other useful 

data included aerial photographs in the United States from 

1935 and Landsat scenes of the entire study area for 1973, 

1986, 1992, and 1997 (Kepner and Edmonds, 2002). Total 

rates of irrigation withdrawals during 1930–86 were derived 

from previous estimates in the United States (Freethey, 1982; 

Corell and others, 1996; Goode and Maddock, 2000). Total 

rates of agricultural withdrawal were estimated on the basis 

of 70-percent irrigation efficiency for nondeficit irrigation 

and 100-percent efficiency for deficit irrigation. The residual 

30 percent of withdrawals was assumed to infiltrate and return 

to the aquifer. Consumptive-use rates in the United States 

were assumed to be 1.00 m/yr for all agriculture, except for 

some fields in the Palominas and Hereford area that were 

determined to be deficit irrigated in 2002 (Arizona Department 

of Water Resources, 2005). Deficit irrigation implies that all 

applied irrigation water is transpired by the crops, and no 

water recharges the aquifer. Consumptive-use rates for fields 

under deficit irrigation in 2002 were assumed to be 0.79 m/yr 

during 1986–2002 (Arizona Department of Water Resources, 

2005) and 1.00 m/yr prior to 1986.

The earliest agriculture in the Sierra Vista subwatershed 

was likely in the Warren Ranch area south of Bisbee, which 

was irrigated with water pumped from the Copper Queen 

mine at Bisbee beginning in 1905 (Southwest Groundwater Consultants, 2004). Aerial photos indicate that several 

hundred hectares (a few hundred acres) of agricultural land 

near Hereford may have been actively irrigated in 1935 

(Arizona Department of Water Resources, unpublished aerial 

photographs, 2005) and, therefore, that withdrawals from the 

regional aquifer in this area began before 1935. The initiation 

of irrigation in the Hereford area is uncertain, and for this 

analysis agricultural withdrawals were assumed to have begun 

in 1930. Irrigation of other areas in the United States and 

Mexico was assumed to have begun in 1940. Rates of total 

ground-water withdrawals for agriculture in the United States 

portion of the study area was about 2 hm3/yr (1,600 acre-ft/yr) 

or less before 1950, about 3–4 hm3/yr (2,400 to 

3,200 acre-ft/yr) during 1950–63, about 3–7 hm3/yr 

(2,400– 5,700 acre-ft/yr) during 1964–85, and about 3 hm3/yr 

(2,400 acre-ft/yr) during 1986–2002 (fig. 6).

Agricultural withdrawal rates in the Mexico portion of 

the study area were derived from Esparza (2002) for the post1980 period and were distributed across agricultural areas near 

Rio Los Fresnos and the San Pedro River defined by Kepner 

and Edmonds (2002). Mexico agriculture during 1940–80 was 

assumed to be in the same locations documented by Kepner 

and Edmonds (2002) and to increase at the same rate as 

agriculture in the United States. Rates of total ground-water 

withdrawals for domestic and agricultural uses in Mexico 

increased about 3.5 hm3/yr (2,800 acre-ft/yr) in 1980 to more 

than 10 hm3/yr (8,100 acre-ft/yr) during 1999–2002 (fig. 6).

Deliveries of withdrawals for domestic and industrial uses 

in incorporated areas of the United States increased from less 

than 2 hm3/yr (1,600 acre-ft/yr) before 1960, to about 5 hm3

(4,100 acre-ft) in 1960, to 12 hm3 (9,700 acre-ft) in 1989, and 

to 14.9 hm3 (12,100 acre-ft) in 2002 (fig. 6). Incorporated 

areas in the United States include Sierra Vista, Huachuca City, 

Tombstone, and Bisbee (fig. 1).

Fort Huachuca withdrawals have been documented 

since 1940. Most of the wells are near Highway 90 and 

a few are across the eastern range of the facility (fig. 1). 

Early withdrawals included wells east of Highway 90. 

Later withdrawals were at wells near Highway 90. Annual 

rates of withdrawals increased from less than 1 hm3/yr 

(800 acre-ft/yr) before 1951, to 3.6 hm3 (2,900 acre-ft) 

in 1951, to a maximum of 4 hm3 (3,200 acre-ft) in 1989, 

and decreased thereafter to less than 2 hm3 (1,600 acre-ft) in 

2002 (fig. 6).

Domestic withdrawals in unincorporated areas in most 

of the United States portion were estimated for later years 

by using Arizona Department of Water Resources reports 

(1990–2005) and estimated for earlier years by calculating the 

difference between estimated withdrawals for other uses and 

Conceptual Model of the Ground-Water Flow System 17
total withdrawals. Unincorporated withdrawals in the United 

States increased from less than 1 hm3/yr (800 acre-ft/yr) before 

1950 to about 2 hm3/yr (1,600 acre-ft/yr) during 1968–77 and 

4.9 hm3/yr (4,000 acre-ft/yr) in 2002 (fig. 6).

Simulation of Ground-Water Flow

Ground-water flow was simulated by using MODFLOW 

2000 (MF2K; Harbaugh and others, 2000 a,b). MF2K uses 

finite-difference numerical methods to solve the ground-water 

flow equation throughout a discretized three-dimensional 

grid that represents the aquifer system. MF2K is a modular 

modeling system that has many processes and packages 

available to simulate many possible conditions that might 

affect ground-water flow. Simulation of the ground-water 

flow system in the study area required several packages for 

the processes of recharge, drains, general-head boundaries, 

streams (Prudic, 1989), and evapotranspiration.

Previous Models

Several published numerical models simulate groundwater flow in the Upper San Pedro Basin. The earliest was 

a two-dimensional model of ground-water flow within 

the alluvial aquifer of the Sierra Vista subwatershed during 

1940–78 (Freethey, 1982). This model was a generalized 

model that was one of many constructed to compare groundwater flow systems in alluvial basins in southern Arizona as 

part of the U.S. Geological Survey Southwest Alluvial Basins 

Regional Aquifer-System Assessment (Anderson and others, 

1992). All later numerical models use the same or similar 

assumptions and conceptual model of the ground-water flow 

system as the initial model by Freethey (1982). Later models 

extended the simulation period beyond 1978, however, and 

included a more detailed horizontal and vertical discretization 

of the system. Later numerical ground-water models of 

the Sierra Vista subwatershed included those developed by 

Vionnet and Maddock (1992) and Corell and others (1996). 

Goode and Maddock (2000) simulated ground-water flow 

in the alluvial aquifer of the entire Upper San Pedro Basin, 

including Mexico. Esparza (2002) simulated ground-water 

flow in the Mexico portion of the alluvial aquifer.

Lack of information required that previous models 

incorporate a simplified conceptualization of the ground-water 

flow system. All recharge was assumed to occur near the 

mountains and to be constant through time. The influence of 

silt and clay layers on ground-water flow and the distribution 

of gaining and losing reaches of the rivers was not recognized. 

Little information was available to describe the vertical 

distribution of hydraulic head and ground-water flow. As a 

result, the simulation of the vertical distribution of groundwater flow was not important in previous models. The vertical 

component of ground-water flow was simulated, however, by 

using as many as four layers. Ground-water flow through the 

permeable rocks in the mountains and flow through basin fill 

that underlies the silt and clay was not simulated in previous 

models. Springs that discharge from the regional aquifer 

and carbonate rocks were not simulated. Stress conditions, 

including ground-water withdrawals and stream baseflow, 

were simulated as average-annual or average conditions during 

multiple years.

Model Framework

Model Boundaries

The area of the simulated ground-water flow system 

includes the portion of the Upper San Pedro Basin in Mexico, 

the Sierra Vista subwatershed in the United States, and a 

small portion of the Benson Subwatershed that lies north of 

the Babocomari River and between the Mustang Mountains 

and the San Pedro River. Simulated boundary conditions 

include no-flow at the lower model boundary and along 

most of the lateral model boundaries and specified head at 

the subwatershed outflow (fig. 5). The upper boundary is 

simulated as a water-table aquifer with evapotranspiration, 

streams, and drains. The upper boundary of the model is the 

land surface interpolated from 60-m DEMs for Mexico, and 

1-m, 10-m, and 30-m DEMs for the United States. One-meter 

DEMs were only available for areas near the San Pedro River 

in the United States. Ten-meter DEMs were used for most of 

the area near the Babocomari River and along some major 

tributaries to the San Pedro and Babocomari Rivers.

No-flow conditions were assumed to occur at more than 

1,500 m below land surface, at the surface-water divide of 

the upper San Pedro Basin, and at outcrops of impermeable 

rock near the San Pedro River. No-flow conditions also were 

assumed to occur along two ground-water flow paths at the 

northern boundary of the model. The region between the 

outcrops of impermeable rock near the San Pedro River and 

the drainage divide in the Whetstone Mountains was simulated 

as a no-flow boundary because water-level data indicated 

that the boundary approximates a ground-water flow path 

(fig. 5), and gravity data indicated a bedrock high exists in the 

region (Gettings and Houser, 1995). A portion of the northern 

subwatershed boundary that extends from near the Dragoon 

Mountains to the streamflow-gaging station near Tombstone 

was simulated as a no-flow boundary by using water-level 

data that indicated it approximates a ground-water flow path 

(fig. 5).

A portion of the boundary at the northern extent of the 

model was simulated as a specified-head flow boundary. 

Interpolated water-level data were used to estimate the 

specified head. Head values were assumed to be constant 

throughout the simulation period because water levels varied 

less than 2 m at wells (D-19-21)12dbb2 and (D-20-20)02ddd 

(fig. 5) during 1960–2001 with no apparent trend (Barnes and 

Putman, 2004). 

18  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Spatial and Temporal Discretization

The model domain was discretized into grids in the 

horizontal and vertical dimensions. Horizontal discretization 

was 250 m throughout the model extent resulting in grids 

that included 440 rows and 320 columns (fig. 5). The model 

domain was discretized vertically among five layers where 

ground-water flow is primarily horizontal through laterally 

extensive hydrogeologic units (fig. 3). The spatial extent of 

model layers was determined by the saturated extent of the 

hydrogeologic units. Extents generally were limited to regions 

of each hydrogeologic unit where the saturated thickness was 

a minimum of 10 m. Saturated thickness was based on winter 

2002 water-level data for the United States and on water 

levels in Mexico reported in a hydrogeology report of the 

area (Consultores en Agua Subterranea S.A. for Mexicana de 

Cananea, S.A. de C.V., 2000). Layer 5 encompassed the entire 

simulated area. Each successively shallower layer included a 

shallower hydrogeologic unit and a less extensive area. The 

extent of each layer was nested within the extent of the next 

deeper layer to allow ground-water flow between vertically 

adjacent layers. 

The simulation of transient ground-water conditions 

included the earliest known withdrawals from the aquifer that 

began in December 1902 (Hollyday, 1963) and continued 

through February 2003. Transient stresses included artificial 

recharge and seasonal ground-water withdrawals by 

evapotranspiration and wells. Each of the transient stresses 

varied in magnitude, which was represented by temporal 

discretization throughout the transient simulation. Temporal 

discretization in MF2K included stress periods and time 

steps. Stress periods for this model included two seasons 

per year. Spring/summer conditions were from March 12 to 

October 15 of each year. The remaining period represented 

fall/winter conditions. The greatest annual rate of groundwater withdrawal for evapotranspiration and wells was during 

the spring/summer period. The fall/winter period had low 

rates of evapotranspiration, no agricultural ground-water 

withdrawals, and reduced rates of withdrawal for domestic and 

municipal supply. The ground-water flow equation was solved 

by MF2K for each of 10 time steps in each seasonal stress 

period. Multiple time steps allowed more accurate simulation 

of the transient system response to variations in stresses and 

reduced computational errors in solving the ground-water 

flow equation.

Aquifer Hydraulic Properties

Key properties of the aquifer that control groundwater flow and response of the ground-water flow system 

to stresses were saturated thickness, horizontal hydraulic 

conductivity (Kh), vertical hydraulic conductivity (Kv), 

specific storage, and specific yield. Each property, except 

Kv, was spatially distributed throughout the model based on 

lithology represented by polygons in each model layer. Kv

was calculated in the model by assuming a vertical anisotropy 

(Kh/Kv) for each lithology. Properties of each rock type vary 

within a range of values developed on the basis of published 

values. Additional detail beyond the lithologic polygons in the 

distribution of hydraulic properties was added by subdividing 

some polygons when hydrologic information, such as spatial 

variation in hydraulic gradients, indicated a need for additional 

variation in hydraulic properties.

Layer Distribution

The five model layers represent the primary 

hydrogeologic units and several lithologies within each layer 

(figs. 3 and 7). Layer 5 represents pre-Tertiary granite and 

metamorphic rock, limestone, sandstone, siltstone, mudstone, 

conglomerate, volcanic rock, and Tertiary Pantano formation 

and basin fill in areas along the margin of the alluvial basin 

where subsurface data are insufficient to define the base of the 

alluvium. Layer 4 represents sand and gravel of lower basin 

fill that lie adjacent to and underlie the thick interval of silt 

and clay within basin fill. Layer 4 includes stream alluvium 

and the silt and clay facies, sand, and gravel of basin fill in 

Mexico. Layers 2 and 3 represent the silt and clay facies of 

upper and lower basin fill, respectively, and include adjacent 

sand and gravel facies and interbedded facies within the Sierra 

Vista subwatershed. Layers 1, 2, and 3 are not defined in 

Mexico because of a lack of subsurface information. Layer 1 

includes the stream alluvium and shallow unconfined ground 

water in the sand and gravel overlying the silt and clay of the 

upper basin fill within the Sierra Vista subwatershed. Each 

layer, except layer 1, was convertible between confined and 

unconfined ground-water flow conditions where no overlying 

layer exists, or where overlying layers were desaturated during 

the transient simulation. Layer 1 was simulated as a watertable aquifer throughout the extent of the layer.

The thickness of each layer was a function of the base 

and top surfaces of each layer and the water table. The land 

surface forms the top surface of each layer where no overlying 

layer exists; otherwise, the base of the overlying layer forms 

the top of each layer. The simulated base of layer 5 was 

uniformly 1,500 m below land surface. The base of layer 4 

includes a large trough in the Sierra Vista subwatershed and 

a small trough in Mexico (fig. 7B). The large trough was 

asymmetric and had its lowest areas near the San Pedro River 

and west of the river where the river traverses the Tombstone 

Hills. The deepest areas of the large trough in layer 4 are as 

much as 500 m below the margins of the trough. The smaller 

trough in Mexico lies between the Sierra Mariquita and Sierra 

San Jose and has deep areas that are 350 m below the margins. 

The base of layer 3 includes three troughs (fig. 7C). The 

largest trough in layer 3 falls near and west of the San Pedro 

River in two areas: the Palominas-Hereford area and the area 

near Highway 90 (fig. 1), where the base of the trough was as 

much as 200 m lower than along the margins. Other troughs in 

layer 3 are north of the Babocomari River and near Fairbank 

where the lowest parts of the troughs were as much as 130 m 

lower than along the margins. The base of layer 2 includes 

Simulation of Ground-Water Flow 19
three troughs that are similar to those in layer 3 (fig. 7D). The 

largest trough in layer 2, near and west of the San Pedro River, 

includes the lowest areas that are as much as 200 m lower in 

the base than along the trough margins. The two other troughs, 

north of the Babocomari River and near Fairbank, include 

low areas where the base is as much as 130 m lower than the 

base along the trough margins (fig. 7D). The base of layer 1 

(fig. 7E) slopes gradually upward from north to south along 

the San Pedro River at about the gradient of the stream. The 

base of layer 1 outside of the San Pedro River generally was 

flat where the layer represents sand and gravel of the basin fill 

near Palominas. Near Highway 90 and west of the San Pedro 

River, the base of layer 1 was as much as 60 m lower than it 

was near Palominas.

Saturated thicknesses of the layers are less than or equal 

to the maximum simulated thicknesses of 1,500 m in layer 5, 

400 m in layer 4, 170 m in layer 3, less than 300 m in layer 

2, and generally less than 100 m in layer 1. Where layer 5 

was overlain by layer 4, the simulated thickness was between 

1,100 and 1,490 m. Where layer 4 was overlain by layer 3, the 

simulated thickness was between 100 and 300 m. Elsewhere, 

layer 4 ranges from 10 to 400 m in thickness. The simulated 

thickness of layer 3 was limited by the base of layer 2 

throughout the layer extent and varies between 10 and 130 m. 

Layer 2 was generally less than 180 m in thickness except in 

a region near Palominas and Hereford where the layer was 

as much as 310 m thick. The extent of layer 1 was generally 

less than 50 m in thickness, except where the layer represents 

sand and gravel overlying the silt and clay facies in the area 

near Highway 90 and west of the San Pedro River where the 

layer was more than 100 m thick. Layer 1 was 3–10 m thick 

along the San Pedro River where the layer represents the 

stream alluvium.

Transmissive Properties

Calibrated hydraulic conductivity values range from 

0.0001 to 12.5 m/d (table 3). Granite and metamorphic 

rocks were assigned values ranging from 0.0001–0.050 m/d. 

Hydraulic conductivity estimates for the pre-Tertiary 

sedimentary rocks resulted from investigations in the Bisbee 

(Phelps Dodge Corporation, 1998) and Tombstone (Hollyday, 

1963) areas. In the Bisbee area, hydraulic conductivity of 

the sandstone was simulated as 6.0 m/d, and the simulated 

hydraulic conductivity of the mudstone, siltstone, and 

conglomerate ranged from 0.18–0.9 m/d throughout 

an effective thickness of about 1,500 m (Phelps Dodge 

Corporation, 1998). In the Tombstone area, the simulated 

transmissivity of the sedimentary rocks ranged from about 

900–11,000 m2/d. An equivalent hydraulic conductivity, 

assuming an effective saturated thickness of about 1,500 m, 

was about 0.6–0.7 m/d. Calibrated hydraulic conductivity 

of pre-Tertiary sedimentary and volcanic rocks for this 

model included a range of values from 0.0001 to 0.625 m/d. 

Hydraulic conductivity of sedimentary rocks was 0.0006–

0.625 m/d for limestone and 0.0001–0.15 m/d for other 

sedimentary rocks that include low-permeability mudstone and 

permeable conglomerate. Pre-Tertiary volcanic rock hydraulicconductivity values ranged from 0.0001 to 0.063 m/d. Basin 

fill included a range of values from 0.0002 to 10.0 m/d for 

the sand and gravel facies. Silt and clay facies of the lower 

basin fill were assigned values of 0.001–0.01 m/d. Intervals in 

the upper basin fill dominated by silt and clay were assigned 

values of 0.05–1.0 m/d. Regions of interbedded basin fill were 

assigned hydraulic-conductivity values of 0.01–4.0 m/d in 

the lower basin fill and 0.02–4.0 m/d in the upper basin fill. 

Sand and gravel within the lower basin fill (conglomerate) 

were assigned values of 0.0001–6.25 m/d. Sand and gravel in 

the upper basin fill were assigned values ranging from 0.05 to 

7.0 m/d. Stream alluvium was assigned values of 7.5 m/d for 

both the pre-entrenchment and post-entrenchment alluvium.

No horizontal anisotropy was simulated throughout the 

model, and vertical anisotropy values varied with lithology. 

Vertical anisotropy values (Kv/ Kh) ranged from 3.5 in some 

parts of the aquifer system dominated by sand and gravel to 

122.5 in silt and clay facies of basin fill (table 3). Sedimentary, 

granitic, and metamorphic rocks were assigned vertical 

anisotropy values ranging from 5–17.5. Basin fill vertical 

anisotropy values ranged from 3.5 to 75 for sand and gravel 

facies, 8.8 to 122.5 for interbedded facies, and 3.5 and 122.5 

for silt and clay facies, respectively. Stream alluvium was 

assigned anisotropy values ranging from 3.5 to 22.5.

The overall ability of the aquifer system to transmit 

ground water was determined by the spatial distribution of 

transmissivity (fig. 8), which was the product of saturated 

layer thickness and hydraulic conductivity. The overall 

transmissive properties of the regional aquifer estimated for 

this model are similar to those of previous models. Values of 

total transmissivity range from about 1 to 2,000 m2/d. Low 

values, less than 50 m2/d, generally coincide with regions of 

sedimentary and crystalline rock and parts of the basin fill 

that have large thicknesses of silt and clay. Regions having 

transmissivity values that range from 50 to 1,000 m2/d 

are associated with regions of basin fill surrounding the 

regions of silt and clay and thick sequences of sandstone 

and conglomerate. The highest values of transmissivity, 

1,000–3,000 m2/d, are associated with the thick regions of 

sand and gravel that surround the regions of silt and clay.

Storage Properties

Two types of storage properties were assigned to model 

cells to simulate the movement of water into and out of storage 

in pore spaces with variation in hydraulic head. Specific-yield 

values were assigned to represent the ratio of drainable pore 

spaces to total cell volume when the water table intersects the 

cell. Specific-storage values also were assigned to each cell 

to represent the change in storage of water that results from 

expansion and contraction of the aquifer skeleton and water 

within the aquifer. 

20 Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Figure 7. Extent and lithology of hydrogeologic units and altitude of the base of numerical model layers, Upper San Pedro Basin, 

United States and Mexico. 7A, Model layer 5; 7B, Model layer 4; 7C, Model layer 3; 7D, Model layer 2; 7E, Model layer 1.

30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

UNITED STATES

MEXICO

Sierra Vista

Cananea

Tombstone

110o30’

15’

15’

31o

110o

TERSH D BOUNDARY

E

WA

P

de ro

R

vi

er

San

Geology modified from Reynolds, 1988

STREAM ALLUVIUM

BASIN FILL SAND, GRAVEL, SILT,

 AND CLAY

VOLCANIC ROCK

SANDSTONE, SILTSTONE, MUDSTONE,

 AND CONGLOMERATE

LIMESTONE

GRANITE AND METAMORPHIC

 ROCK

EXPLANATION

Figure 7. Extent of lithologies and altitude of the base of: A, Model layer 5; B, Model layer 4; C, Model layer 3; D, Model

layer 2; E, Model layer 1, Upper San Pedro Basin, United States and Mexico.

MODEL GRID ROW 1

MODEL GRID ROW 440

MODEL GRID COLUMN 1

MODEL GRID COLUMN 320

A. Model layer 5

Extent of model grid

NOTE: Base of layer 5 is 1,500 meters below land surface

Bisbee

Simulation of Ground-Water Flow  21
30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

UNITED STATES

MEXICO

Sierra Vista

Bisbee

Tombstone

110o30’

15’

15’

31o

110o

TERSH D BOUNDARY

E

WA

STREAM ALLUVIUM AND BASIN

 FILL

BASIN FILL SAND AND GRAVEL

BASIN FILL SILT AND CLAY

STRUCTURE CONTOUR—Shows altitude

 of base of layer 4. Contour interval 100

 meters. Datum is NAVD 88

EXPLANATION

Figure 7. Continued.

MODEL GRID ROW 1

MODEL GRID COLUMN 1

MODEL GRID COLUMN 320

B. Model layer 4

900

800

1

1,200

1,200 1,200

1,3 0 0

1,3 00

9 0 0

,000

1,100

1,200

900

1,100

1,000

1,200

1,200

1,10 0

01, 00

MODEL GRID ROW 440

Cananea

900

P

de ro

River

San

Extent of model grid

Figure 7. Continued.

22  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Figure 7. Continued.

30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

UNITED STATES

MEXICO

Sierra Vista

Bisbee

Tombstone

110o30’

15’

15’

31o

110o

TERSH D BOUNDARY

E

WA

BASIN FILL SAND AND GRAVEL

BASIN FILL SAND, GRAVEL, SILT

 AND CLAY

BASIN FILL SILT AND CLAY

STRUCTURE CONTOUR—Shows altitude

 of base of layer 3. Contour interval 100

 meters. Datum is NAVD 88

EXPLANATION

Figure 7. Continued.

MODEL GRID ROW 1

MODEL GRID COLUMN 1

MODEL GRID COLUMN 320

C. Model layer 3

1,2 0

1,0 0 0

1,100

1,200

1,100

1,100

1 1, 00

1, 1 0 0

MODEL GRID ROW 440

Cananea

1,100

0

P

de ro

River

San

Extent of model grid

Simulation of Ground-Water Flow  23
Figure 7. Continued.

30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

UNITED STATES

MEXICO

Sierra Vista

Tombstone

110o30’

15’

15’

31o

110o

TERSH D BOUNDARY

E

WA

P

de ro

River

San

BASIN FILL SAND AND GRAVEL

BASIN FILL SAND, GRAVEL,

 SILT, AND CLAY 

BASIN FILL SILT AND CLAY

STRUCTURE CONTOUR—Shows altitude

 of base of layer 2. Contour interval 100

 meters. Datum is NAVD 88

EXPLANATION

Figure 7. Continued.

MODEL GRID ROW 1

MODEL GRID COLUMN 1

MODEL GRID COLUMN 320

D. Model layer 2

MODEL GRID ROW 440

Cananea

1,100

,100 1

1,200

1,200

1,200

1,200

1,10 0

1,2 00

11, 00

1,000 Bisbee

Extent of model grid

24  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Figure 7. Continued.

30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

UNITED STATES

MEXICO

Sierra Vista

Tombstone

110o30’

15’

15’

31o

110o

TERSH D BOUNDARY

E

WA

P

San

dro

R

vi

er

STREAM ALLUVIUM AND BASIN

 FILL

Unsaturated zone

STRUCTURE CONTOUR—Shows altitude

 of base of layer 1. Contour interval 50

 meters. Datum is NAVD 88

EXPLANATION

Figure 7. Continued.

MODEL GRID ROW 1

MODEL GRID COLUMN 1

MODEL GRID COLUMN 320

E. Model layer 1

MODEL GRID ROW 440

Cananea

1,250

1,300

1,200

1,250

1,300

1,300

e

1,250

200

1,

Bisbee

Extent of model grid

Simulation of Ground-Water Flow  25
Calibrated storage properties (table 3) ranged from a 

specific storage coefficient of 1.0 x 10-6 m-1 for confined 

portions of the model layers to specific-yield values of 0.30 

in unconfined portions of the model layers. A specificstorage coefficient value of 1 x 10-6 m-1 was assigned to the 

sedimentary and crystalline rocks in layer 5. Lower basin 

fill, simulated by layers 3 and 4, was assigned values of 

1 x 10-6 m-1–1 x 10-5 m-1. Upper basin fill, simulated by 

layer 2, was assigned values of 1 x 10-6 m-1–2 x 10-5 m-1. 

Stream alluvium, simulated by layer 1, was assigned values 

of 1 x 10-6 m-1–5 x 10-6 m-1. Specific yield for crystalline 

rocks and low-permeability sedimentary rocks in layer 5 

was assigned values of 0.001–0.01. Limestone was assigned 

a specific-yield value of 0.01 throughout the extent of the 

layer except in the area near the Tombstone mine, which was 

assigned a value of 0.02. Sandstone was assigned specificyield values of 0.01–0.20. Specific–yield values for lower 

basin fill in layer 4 ranged from 0.05 in areas of consolidated 

basin fill to 0.20 in areas of unconsolidated coarse-grained 

basin fill. Specific-yield values for basin fill in layers 2 and 

3 ranged from 0.05 in areas of silt and clay to 0.30 in areas 

of coarse-grained basin fill. The specific yield of layer 1 

was 0.25 in areas where the layer represents basin fill that 

lies to the west of the San Pedro River, 0.20–0.30 in the 

pre-entrenchment stream alluvium, and 0.30 in the postentrenchment stream alluvium. 

Simulated Inflows and Outflows

Recharge

Recharge to the aquifer system includes natural sources 

and artificial sources that result from human water use. 

Natural recharge occurs in areas where rates of infiltration of 

precipitation and runon are greater than rates of evaporation, 

transpiration, and runoff. Rates of infiltration are influenced 

by slope, near-surface materials, and rates of precipitation. 

Rates of evaporation are influenced by vegetation type 

and elevation. Artificial recharge occurs as excess applied 

irrigation water and as incidental infiltration of water in the 

distribution and sewage systems and as intended infiltration 

of sewage effluent and discharge from mine dewatering. 

Rates of recharge varied throughout the investigation because 

of variations in climate, vegetation, land use, and mining 

practices. Variations in natural recharge were not included in 

the simulation of ground-water flow, but variations in artificial 

recharge were included.

Natural Recharge

Natural recharge has been assumed by previous 

investigations to occur primarily near the mountain fronts 

through infiltration of runoff in ephemeral channels. 

Mountain-front recharge often has been assumed to occur 

within 1.6 km of the bedrock of the mountains. Recent 

investigations, however, indicate that significant infiltration 

also occurs in the limestone in the mountains, along ephemeral 

channels that cross the alluvial slope below the mountain-front 

recharge areas (Coes and Pool, 2005; Goodrich and others, 

2004) and along ephemeral reaches of the San Pedro River 

(Pool and Leenhouts, 2002). Predevelopment natural recharge 

rates were distributed across the watershed for the groundwater flow model by using water-budget methods. The overall 

amount of water flowing through the ground-water system was 

determined by the sum of estimated predevelopment discharge 

to stream baseflow, evapotranspiration, and ground-water flow 

to the Benson subwatershed. The total recharge rate along 

ephemeral stream channels was estimated as 15 percent of 

the total rate of flow through the ground-water system (Coes 

and Pool, 2005). Recharge along ephemeral channels was 

approximated by uniformly distributing 15 percent of the 

total recharge throughout the alluvial deposits that lie below 

mountain-front recharge areas. A spatial average rate of annual 

recharge to sedimentary rocks was estimated as the average 

annual baseflow discharge per unit area from four watersheds 

in the Huachuca Mountains: Garden Canyon, Huachuca 

Canyon, Ramsey Canyon, and the Upper Babocomari River. 

Mountain-front recharge was distributed within 1.6 km of the 

mountains as the residual of the total recharge and the sum 

of total estimated recharge along ephemeral channels below 

the mountain-front recharge area and recharge to sedimentary 

rock. Simulated rates of predevelopment natural recharge 

were set to 0 m/d in areas of crystalline rock and ranged from 

1.0 x 10-6 to 2 x 10-5 m/d in areas of sedimentary rocks, 

2.0 x 10-5 to 1.35 x 10-4 m/d within 1.6 km of the mountain 

fronts. The remainder of the alluvial surface was assigned a 

rate of 1.0 x 10-6 m/d.

Artificial and Incidental Recharge

Artificial recharge in the study area results from of 

infiltration of sewage effluent near Sierra Vista, Naco, Bisbee, 

Fort Huachuca, and Tombstone; waste water from individual 

domestic septic systems; and water discharged from the 

mines near Bisbee and Tombstone for dewatering purposes. 

Incidental recharge occurs as excess irrigation water and 

during the distribution of potable water withdrawn by wells. 

Each source of artificial and incidental recharge was simulated 

either explicitly as a transient recharge rate or implicitly as 

reduced withdrawal rates. Incidental recharge from distribution 

systems was implicitly simulated through the use of reported 

deliveries to customers by water providers rather than actual 

total well withdrawals.

Recharge though deep percolation of sewage effluent 

was simulated for sewage treatment facilities and individual 

septic systems. Sewage treatment facilities exist for Sierra 

Vista, Naco, Bisbee, Fort Huachuca, Huachuca City, and 

Tombstone. Recharge was distributed near each treatment 

facility on the basis of the method of discharge. Early treated 

effluent typically was discharged into ephemeral channels. 

Improved facilities often used treated effluent to irrigate turf 

and crops, or discharged it to evaporation or infiltration ponds. 

26  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

UNITED STATES

MEXICO

Tombstone

110o30’

15’

15’

31o

110o

TERSH D BOUNDARY

E

WA

Figure 8. Simulated total transmissivity of the numerical model layers, Upper San Pedro Basin, United States and Mexico.

MODEL GRID ROW 1

MODEL GRID COLUMN 1

MODEL GRID COLUMN 320

MODEL GRID ROW 440

Cananea

0–10

10–50

50–200

200–1,000

1,000–3,000

EXPLANATION

TOTAL TRANSMISSIVITY, IN METERS SQUARED PER DAY

Sierra Vista

Bisbee

Extent of model grid

Figure 8. Simulated total transmissivity of the numerical model layers, Upper San Pedro Basin, United States and Mexico.

Simulation of Ground-Water Flow  27
In the absence of early treatment records, recharge rates 

were estimated as 14 percent of total delivered withdrawals on 

the basis of the average of recent effluent recharge estimates 

in the Sierra Vista subwatershed (Arizona Department of 

Water Resources, 2005). The City of Sierra Vista began 

treating sewage effluent in 1967, increased the capacity with 

an additional facility in 1978, and converted to a recharge 

facility in 2001. During 1967–79, effluent was discharged 

into a nearby ephemeral channel. Beginning in 1980, effluent 

was disposed of through irrigation of 3.4 x 106 m2 near the 

facility. Beginning in July 2002, effluent was discharged 

to recharge basins and recharged at a rate of 2,960 m3/d 

(Arizona Department of Water Resources, 2005). Effluent was 

discharged to ponds near Naco and Bisbee beginning about 

1930 (Phelps Dodge Corporation, 1998); however, the Bisbee 

effluent ponds were outside of the watershed. Effluent at Fort 

Huachuca during 1978–95 was recharged at two turf sites—a 

golf course and parade grounds—and below the treatment 

facility at about 1,360 m3/d (U.S. Department of Defense, 

2002). During 2001 about 1,825 m3/d, about 33 percent of 

withdrawals, was estimated to have recharged. Treatment 

facilities at Huachuca City discharge to evaporation ponds 

where little recharge occurs. Treatment facilities at Tombstone 

began discharging effluent into Walnut Gulch in 1972 at a 

rate of about 300 m3/d. Recharge at individual septic systems 

was simulated assuming recharge wells return 14 percent of 

withdrawals to the uppermost model layer near each well 

site. The 14 percent value was based on estimates made for a 

decision support system developed for the Upper San Pedro 

Partnership (Kevin Lansey, Professor, University of Arizona 

Department of Civil Engineering and Engineering Mechanics, 

written commun., 2005). Recharge in municipal areas before 

Table 3. Simulated hydraulic properties of hydrogeologic units for the ground-water flow model of the Upper San Pedro basin, 

United States and Mexico.

[m, meter; m/d, meter per day; Kh, horizontal hydraulic conductivity; Kv, vertical hydraulic conductivity; max, maximum; min, minimum; average values 

represent the average of values assigned to several polygons; average values are not area weighted]

Hydrogeologic unit

 Hydraulic conductivity Vertical anisotropy Specific storage Specific yield m/d Kh

/Kv m-1

Average max min Average max min Average max min Average max min

Limestone 0.072 0.625 0.0006 7.9 17.5 5.0 1.0x10-6 1.0x10-6 1.0x10-6 0.011 0.020 0.010

Sedimentary rocks .039 .300 .0001 7.3 17.5 5.0 1.0x10-6 1.0x10-6 1.0x10-6 .088 .200 .010

Granitic and 

metamorphic rocks

.006 .050 .0001 8.8 17.5 5.0 1.0x10-6 1.0x10-6 1.0x10-6 .006 .010 .001

Volcanic rocks .018 .063 .0001 7.1 17.5 5.0 1.0x10-6 1.0x10-6 1.0x10-6 .002 .010 .001

Undifferentiated basin fill

Undifferentiated sand 

and gravel1

.795 10.000 .0003 9.4 17.5 3.5 2.4x10-6 5.0x10-6 1.0x10-6 .120 .200 .010

Undifferentiated silt 

and clay1

.285 1.250 .0013 3.5 3.5 3.5 5.0x10-6 5.0x10-6 5.0x10-6 .250 .250 .250

Upper basin fill

Sand and gravel 3.459 7.000 .0500 26.8 75.0 8.8 1.5x10-5 2.0x10-5 1.0x10-6 .177 .300 .100

Interbedded .887 4.000 .0200 27.3 87.5 8.8 2.0x10-5 2.0x10-5 2.0x10-5 .150 .250 .050

Silt and clay .229 1.000 .0500 65.0 87.5 8.8 2.0x10-5 2.0x10-5 2.0x10-5 .057 .100 .050

Lower basin fill

Sand and gravel .979 6.250 .0002 10.8 36.1 3.5 3.5x10-6 5.0x10-6 1.0x10-6 .119 .200 .010

Interbedded .785 4.000 .0100 38.2 122.5 12.3 6.7x10-6 1.0x10-5 5.0x10-6 .092 .100 .050

Silt and clay .005 .010 .0010 122.5 122.5 122.5 6.3x10-6 1.0x10-5 5.0x10-6 .050 .050 .050

Stream alluvium

Undifferentiated1 4.929 12.500 2.5000 8.9 22.5 3.5 3.9x10-6 5.0x10-6 1.0x10-6 .264 .300 .250

Pre-entrenchment 7.500 7.500 7.5000 22.5 22.5 22.5 1.0x10-6 1.0x10-6 1.0x10-6 .291 .300 .200

Post-entrenchment2 7.500 7.500 7.5000 7.5 7.5 7.5 1.0x10-6 1.0x10-6 1.0x10-6 .300 .300 .300

1

Hydrogeologic units designated as undifferentiated represent units in layers 4 and 5. 

2

Values represent a single polygon.

28  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
construction of treatment facilities also was assumed to occur 

as deep percolation at individual septic systems and was 

simulated by assuming that a recharge well was near each 

withdrawal well.

Recharge of excess irrigation water has primarily 

occurred in the agricultural areas near the San Pedro River 

near Palominas and in Sonora, Mexico. Recharge of irrigation 

water also occurs beneath golf courses. Rates of recharge 

depend on irrigation practices and efficiencies. Previous 

ground-water flow models of the area have assumed that 

agricultural irrigation efficiencies are about 70 percent, 

which results in 30 percent of the irrigation withdrawals 

returning to the aquifer through deep percolation (Freethey, 

1982; Corell and others, 1996; Goode and Maddock, 2000). 

Deep percolation of excess applied irrigation water was 

simulated in the numerical model as recharge to the upper 

most layer at each well by using an additional recharge well. 

Deep percolation beneath golf courses near Naco, Sierra 

Vista, and Fort Huachuca was simulated by applying a 

recharge rate of 1.4 x 10-3 m/d to the turf area (Phelps Dodge 

Corporation, 1998).

Evapotranspiration of Ground Water

Evapotranspiration from riparian vegetation was 

simulated near the San Pedro and Babocomari Rivers on the 

basis of recent studies that improve definition of extents of 

riparian vegetation and evapotranspiration rates of several 

riparian vegetation types (Leenhouts and others, 2005). 

Vegetation type across the model area has been mapped to 

50 m grids from Landsat satellite photographs taken in 1973, 

1986, 1992, and 1997 (Kepner and Edmonds, 2002) and GIS 

coverages of vegetation type near the San Pedro River in the 

United States developed by using aerial photos in 2001 (U.S. 

Army Corps of Engineers, 2001).

Simulation of evapotranspiration by using MODFLOW 

requires three types of input data: evapotranspiration surface 

elevation, maximum evapotranspiration rate, and extinction 

depth. Evapotranspiration surface elevation was commonly 

taken to be the land surface derived from DEMs. Maximum 

evapotranspiration rates and extinction depths for various 

phreatophyte types are available from published reports. 

Extinction depth was the maximum depth to water below 

the evapotranspiration surface where the phreatophytes 

can withdraw water from the aquifer. Evapotranspiration 

rate was calculated in MODFLOW by using the modelcalculated water-level altitude and depth to water below 

the evapotranspiration surface. Where the water-level 

altitude was at or above the evapotranspiration surface, the 

simulated evapotranspiration rate equaled the maximum 

rate. No evapotranspiration occurred where the depth of the 

calculated water level below the evapotranspiration surface 

was greater than the extinction depth. The evapotranspiration 

rate was calculated from a linear relation between the two 

extremes where the depth of the calculated water level below 

the evapotranspiration surface was between zero and the 

extinction depth.

The spatial distribution of evapotranspiration (fig. 5) 

was derived from mapped distributions of the dominant 

phreatophyte vegetation types of mesquite and riparian 

woodland, which was predominantly cottonwood, willow, and 

sycamore. An additional phreatophyte type, sacaton grass, was 

not explicitly simulated because maps of the distribution were 

not available for areas in Mexico or along the Babocomari 

River. In addition, the evapotranspiration extinction depth 

for sacaton grass was 3 m below land surface, which was 

greater than the depth to water across most of the area. The 

phreatophyte types were mapped as polygonal areas. Each 

polygon was assigned maximum evapotranspiration rates 

and extinction depths that were mapped to the model grid. 

The evapotranspiration surface was assigned to the model 

grid as 1.5 m below average land-surface altitude within 

each 250 x 250-m model grid cell. The evapotranspiration 

surface was below the average land surface within a grid cell 

for two reasons: (1) the evapotranspiration rate was greatest 

where the water table was more than a meter below the land 

surface for most phreatophyte vegetation types, and (2) most 

phreatophytes were at the lower elevations within grid 

cells. High resolution LIDAR data (1 m2) within the area of 

evapotranspiration indicate that minimum land surface within 

each model cell was 1–4 m below the average land surface in 

each model cell. Each evapotranspiration cell also required 

an assigned model layer from which evapotranspiration rates 

were calculated. Evapotranspiration rates were assigned to the 

upper-most active model layer at each evapotranspiration cell, 

which was predominantly layer 1 along the San Pedro River 

in the Sierra Vista subwatershed; layers 2, 4, and 5 along the 

Babocomari River; and layers 4 and 5 in Mexico.

Maximum evapotranspiration rates and extinction 

depths were assigned on the basis of data from recent 

investigations within the SPRNCA (Leenhouts and others, 

2005). Available evapotranspiration rates represent rates for 

100 percent coverage for each vegetation type. Maximum 

rates, therefore, need to be scaled for the percentage of area 

covered by phreatophytes within each model grid cell. Cover 

percentage was not known, however, throughout most of the 

simulation period. Estimation of cover percentage would 

require detailed analysis of aerial photographs taken at several 

times during the simulation period and was not done for 

this analysis. A multiplying factor was, therefore, uniformly 

applied to the maximum rate of evapotranspiration to match 

the estimated total evapotranspiration rates. Maximum rates 

of evapotranspiration for each phreatophyte type for transient 

conditions were 1.3 x 10-3 and 5.2 x 10-5 m/d for mesquite 

and riparian woodland, respectively. Multiplying factors 

and maximum evapotranspiration rates were scaled upward 

by a factor of 3 through the simulated periods to match 

increasing rates of estimated evapotranspiration. Extinction 

depths were 6 m and 14 m for riparian woodland and 

mesquite, respectively.

Simulation of Ground-Water Flow  29
Streamflow Routing

Ground-water discharge to streams was simulated by 

using the MODFLOW streamflow-routing package STR1 

(Prudic, 1989). Input to the package requires stage, elevation 

of the bottom of the streambed, elevation of the top of 

the streambed, conductance, width, slope, sinuosity, and 

Manning’s roughness coefficient. STR1 routes the streamflow 

through a network of channels and uses Manning’s equation 

to calculate stream stage at a hypothetical rectangular channel 

cross section. The application of STR1 does not allow for 

simulation of surface runoff from individual or seasonal runoff 

events; therefore, no stream inflow from upland channels was 

simulated. Improved simulation of streamflow routing requires 

the additional capability of simulating more complex streamchannel geometry, including overbank flow and temporal 

routing of flow.

Streams in the model area were simulated by using a 

network of channels along the San Pedro and Babocomari 

Rivers and several tributaries that include springs that could 

contribute to surface flow in the San Pedro River. The stream 

network extends from the model outflow near the streamflowgaging station near Tombstone on the San Pedro River to north 

of Cananea on the San Pedro River. The simulated portion of 

the Babocomari River extends from the confluence with the 

San Pedro River to near Elgin. Major simulated tributaries 

include the lower reaches of Lyle Creek, a Babocomari 

River tributary, and Rio Los Fresnos in Mexico. Springs 

that issue from the regional aquifer along several ephemeral 

channels, including Coyote Wash (Moson Spring), Murray 

Wash (Murray Spring), Bakarich-McCool Wash (HorseThief Spring), and Government Draw (Lewis Springs), were 

simulated along channels that are tributary to the San Pedro 

River. Streams also were simulated along the lower reaches 

of Garden Canyon Wash and Miller Creek where riparian 

vegetation indicates shallow ground water and possible past 

springflow.

Simulated streams were divided into ninety-nine 

segments and reaches within segments. Each segment 

was assigned top and bottom elevations of the streambed 

at the upper and lower extent of the segment, streambed 

conductance, stream width, sinuosity, and Manning’s 

roughness coefficient. Slope was calculated by STR1 for 

each reach within segments. Elevations of the top of the 

streambed were estimated using DEMs. The most accurate 

DEMs are available from LIDAR surveys completed during 

2003 along the San Pedro River (David Goodrich, Research 

Hydrologist, Agriculture Research Service, written commun., 

2003). The streambed elevation was estimated on the basis 

of the minimum elevation within any 250-m reach of the 

stream. The bottom of the streambed was set at 0.5 m below 

the streambed elevation. Stage initially was set to 1 m above 

the streambed and allowed to vary with head in the adjacent 

model cell. Streambed conductance varied from 0.10–10 m/d. 

Low values were assigned in areas where the stream overlies 

layer 5. Higher values were assigned where the stream overlies 

layer 1. Stream width varied from narrow widths of 1 m in the 

upper stream reaches, to 3 m along the middle reaches of the 

San Pedro River, to 5 m below Highway 90. The number of 

stream vertices specified in STR1 did not adequately define 

the stream sinuosity. Sinuosity was, therefore, estimated for 

each segment by using the ratio of the length of the simulated 

stream segment to the actual stream length and applied to the 

stream length in each cell. Manning’s roughness coefficient 

was 0.22 for all segments.

Ground-Water Withdrawals

Different methods were used to spatially distribute 

simulated withdrawals of ground water for various uses. 

Withdrawals could be assigned to particular well locations 

where withdrawals were reported by well, such as withdrawals 

at Fort Huachuca wells, the Tombstone mines and Copper 

Queen mine, and after 1970 at wells of many large water 

purveyors. Other withdrawals were distributed on the basis 

of available information, such as locations of irrigated 

areas for agricultural uses and existing domestic wells for 

most unincorporated uses. Withdrawals in Mexico for mine 

uses were assigned equally to the approximate locations 

of withdrawal wells. More information is available for the 

distribution of ground water withdrawn before 1940, when 

withdrawals were predominantly for mine use, and for recent 

decades when most withdrawals were reported by large water 

purveyors. 

Rates and distributions of agricultural withdrawals were 

generalized on the basis of irrigated areas prior to the mid1980s and are uncertain for any particular year. Agricultural 

withdrawal rates in Mexico were distributed across 

agricultural areas near Rio Los Fresnos and the San Pedro 

River on the basis of agricultural areas defined by Kepner and 

Edmonds (2002). Irrigation withdrawals in the United States 

before 1987 were assigned to known irrigation wells near the 

irrigated areas. Irrigation withdrawals during 1987–98 were 

assigned to individual wells near agricultural areas that were 

documented by the Arizona Department of Water Resources 

(Arizona Department of Water Resources, 1990). Agricultural 

withdrawals during 1999 through 2002 were assigned to 

wells on the basis of irrigation practices documented in 2002 

(Arizona Department of Water Resources, 2005).

The rates and distributions of withdrawals in 

unincorporated areas are poorly defined throughout the 

simulation period. Withdrawal estimates must be made from 

indirect sources, including population and well-construction 

records. Unincorporated withdrawals in the area between 

Bisbee and Naco were derived from a ground-water flow 

model of the area (Phelps Dodge Corporation, 1998).

Rates of ground-water withdrawals were distributed 

among the two simulated seasons, spring/summer and fall/

winter. The withdrawal seasons were determined on the basis 

of variations in evapotranspiration rates, recent agricultural 

irrigation practices, and monthly variations in deliveries during 

1998–2002 for some of the largest water suppliers: Bella Vista 

30 Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Water Company, Arizona Water Company, Pueblo del Sol 

Water Company, and Fort Huachuca. Public withdrawal rates 

varied seasonally from 84 percent of the average-annual rate 

during the fall/winter season to 115 percent of the annual rate 

during the spring/summer season. Agricultural withdrawals 

were simulated for the spring/summer season. Withdrawals for 

mine use did not vary seasonally.

Knowledge of the vertical distribution of withdrawals is 

important for proper simulation of the transient response of 

the aquifer system to pumping stress, especially in regions 

of confined aquifers. Withdrawals were assigned to model 

layer 5 in areas where layers 1, 2, 3, and 4 do not exist or the 

saturated thickness was thin. Withdrawals from model layer 

5 include withdrawals at the Tombstone mines and Copper 

Queen mine, some withdrawals for domestic use near the 

base of the Huachuca Mountains and in the Bisbee area, and 

some mine withdrawals in Mexico. In parts of the aquifer 

system that are dominated by unconfined ground-water flow 

and where a good hydraulic connection exists among the 

upper and lower basin fill units, withdrawals were assumed to 

derive from the coarse-grained facies of lower basin fill, model 

layer 4. These unconfined areas lie outside of the region of silt 

and clay (fig. 4) and include most withdrawals in the Sierra 

Vista and Fort Huachuca area, withdrawals in the Bisbee 

and Tombstone areas, and most of the domestic withdrawals 

in the region between the base of the Huachuca Mountains 

and the western extent of the silt and clay. The vertical 

distribution of withdrawals was less certain in regions of 

confined ground water where silt and clay exist. Withdrawals 

in areas of confined aquifers were distributed among model 

layers 1, 2, 3, and 4 according to the hydraulic conductivity 

and saturated thickness of the layers with a few exceptions. 

Withdrawals for domestic purposes were assigned to layer 2 

outside of the regions dominated by silt and clay (fig. 4). 

Agricultural withdrawals between Highway 90 and Hereford 

can reasonably be assumed to originate from confined 

portions of the ground-water system on the basis of significant 

variations in hydraulic head in abandoned irrigation wells 

that are 100–300 m deep. Agricultural withdrawals between 

Highway 90 and Hereford were, therefore, assigned to model 

layer 4. Pre-1940 agricultural withdrawals near Hereford 

were assigned to model layer 2. Construction information 

for some agricultural wells near Palominas indicated that the 

withdrawals were derived exclusively from beneath the silt and 

clay layers and, therefore, were assigned to model layer 4.

Initial Conditions

Initial conditions for simulating transient ground-water 

flow approximate steady-state conditions, which were defined 

by a period of negligible changes in inflow, outflow, or 

storage. Hydrologic records indicate that the ground-water 

system in the Upper San Pedro Basin has undergone storage 

change in response to several types of changes in recharge or 

discharge since about 1900. Observed variations in baseflow 

at the Charleston and Palominas streamflow-gaging stations 

during the period of record indicate that the ground-water 

system has changed (Pool and Coes, 1999) and that true 

steady-state conditions have not existed. A period of minimal 

change must, therefore, be chosen as the initial conditions 

to calibrate the steady-state aquifer hydraulic properties and 

flow terms and to simulate the transient variations in the 

ground-water system. Ideally, initial conditions also would 

represent a period before the system changed significantly. 

The temporal and spatial distribution of hydrologic data, 

however, prevented definition of the pre-ground-water 

withdrawal period and required approximation of steady-state 

conditions by using later data. Initial hydrologic conditions 

that existed prior to mine withdrawals beginning in 1902 

were used for this ground-water flow model. The hydrologic 

system in 1902 was assumed to be approximated by later 

data, including streamflow conditions in 1940 and water-level 

distributions during 1940–60. Initial conditions in 1902 allow 

the simulation of all known ground-water withdrawals.

Previous models also assumed that hydrologic conditions 

existing in about 1940 were adequate for the initial conditions 

of the steady-state ground-water flow model. Several changes 

in the ground-water flow system, however, are known to have 

occurred, including lowering of stream elevation through 

stream-channel incision along the major streams (Hereford, 

1993) and associated changes in riparian vegetation and 

evapotranspiration; mine dewatering beginning in 1902 

(Phelps Dodge Corporation, 1998; Hollyday, 1963); groundwater withdrawals for agricultural, stock, domestic, and 

municipal use; and climate-induced recharge variations 

(Dickinson and others, 2004; Pool, 2005). The uncertain 

transient effects of these changes decreases the confidence of 

the adequacy of 1940 as the initial conditions for the steadystate model.

The steady-state hydrologic conditions were defined by 

baseflow and hydraulic gradients, indicated by water levels in 

wells, throughout the ground-water flow system. Streamflow 

monitoring at the Charleston and Palominas streamflowgaging stations documented conditions since 1904 and 1930, 

respectively. The spatial distribution of water-level data, 

however, was not sufficient to define hydraulic gradients 

before about 1961 throughout most of the United States 

portion of the study area and later in Mexico. Hydrographs of 

water levels for eight wells during 1945–60 (fig. 9) indicate, 

however, that water-level declines were minimal during that 

period and, therefore, water levels measured before 1961 can 

be used to approximate 1940 conditions. Water-level data 

for the wells generally varied within a range of about 3 m 

and showed a slight decline, less than 2 m per decade. The 

water level in well (D-23-21)06ccc2 (fig. 9) increased a few 

meters. Only the earliest water levels at any well that had 

multiple measurements before 1961 were used to approximate 

1940 conditions. Water levels in wells in Mexico prior to the 

mid-1980s were not available. Mexico data used to estimate 

steady-state hydraulic gradients were derived from published 

data (Consultores en Agua Subterranea S.A. por Mexicana 

Simulation of Ground-Water Flow  31
de Cananea, S.A. de C.V., 2000). Steady-state water levels 

in wells were insufficient to define vertical ground-water 

gradients throughout the model area. Recent data collected 

in nested piezometers near the river, however, limit the 

magnitude of vertical hydraulic gradients during steady-state 

conditions in those areas.

Model Calibration

Model calibration included adjustment of hydraulic 

properties and some flux terms to match hydrologic 

observations. Hydraulic properties that were adjusted include 

horizontal and vertical hydraulic conductivity, stream-bed 

hydraulic conductivity, and storage properties—specific 

storage and specific yield. Riparian evapotranspiration 

rates were adjusted to match overall estimated rates of 

riparian evapotranspiration throughout the model area. 

The steady-state simulation was used to calibrate hydraulic 

properties—hydraulic conductivity and vertical anisotropy—

by matching simulated with observed steady-state hydraulic 

gradients, and the average-annual steady-state water-budget 

components of recharge rates, stream baseflow, estimated 

ground-water underflow, and ground-water discharge with 

riparian vegetation. Transient calibration included adjusting 

storage properties to match hydrographs of selected wells 

and the overall magnitude of evapotranspiration rates to 

match estimated changes in evapotranspiration. Calibration of 

aquifer-storage properties was achieved by matching slopes of 

long-term water-level trends and magnitude of seasonal waterlevel variations. In addition, vertical anisotropy was adjusted 

to match the seasonal timing and magnitude of the water-level 

and streamflow variations.

Steady State

Methods and data used for steady-state calibration of 

this model are similar to those used for previous models. 

The model was calibrated to average-annual steady-state 

conditions. Ground-water response to seasonal variations in 

evapotranspiration also was simulated. Seasonal variations 

can be in a dynamic oscillatory steady state when storage and 

other flux components vary seasonally, but storage does not 

change during a period of 1 or more years. 

Average-Annual Steady-State Conditions

The estimated (table 2) and simulated (table 4) averageannual steady-state ground-water budgets are in agreement. 

Effects of transient ground-water withdrawals and artificial 

recharge near the Bisbee and Tombstone mines were 

assumed to have had no impact on basin discharge in 1935. 

Simulated recharge was 58,100 m3/d (17,200 acre-ft/yr). 

Simulated net ground-water discharge included 58,700 m3/d 

(17,400 acre-ft/yr) to streams, 1,500 m3/d (400 acre-ft/yr) to 

drains, 2,700 m3/d (800 acre-ft/yr) as ground-water underflow, 

and 22,000 m3/d (6,500 acre-ft/yr) as evapotranspiration. 

Simulated stream baseflow at Charleston was 19,600 m3/d 

(5,800 acre-ft/yr or 8.0 ft3/s) in comparison with the estimated 

annual baseflow of about 19,900 m3/d (5,900 acre-ft/yr 

or 8.1 ft3/s). 

The distribution of predevelopment hydraulic gradients 

and stream baseflow was approximated in the calibration 

process by adjusting the hydraulic properties of streams and 

the hydraulic conductivity of polygonal regions representing 

different hydrogeologic conditions in each layer. The resulting 

distribution of hydraulic properties produced the simulated 

water-level altitudes at predevelopment observation wells 

(fig. 9). Simulated predevelopment hydraulic head was 

generally within 10 m of seventy-seven pre-1961 observed 

values (mean absolute error of 9.0 m and average error of 

-0.7 m) and twenty-nine pre-1950 values (mean absolute error 

of 7.5 m and average error of -3.5 m, fig. 10). The greatest 

error occurred above elevations of 1,440 m near the Huachuca 

Mountains where heads were simulated as 10–40 m lower 

than observed. Hydraulic heads in this part of the model 

are sensitive to variations in hydraulic properties because 

hydraulic gradients are steep and water levels are close to 

the land surface. Water levels in this area varied as much as 

10–30 m with climate variations during 1940–2002 (Dickinson 

and others, 2004; Pool, 2005).

Vertical and horizontal hydraulic properties in three 

areas near Lewis Springs, Cottonwood, and Palominas were 

calibrated by using observed vertical hydraulic gradients near 

the San Pedro River during 1995–2005 (fig. 5). Simulated 

steady-state hydraulic gradients were less than gradients 

that were observed at each site. Observed gradients near 

Lewis Springs were derived from several piezometers that 

were screened at depths of 7–10 m and 60 m within the 

pre-entrenchment alluvium and basin fill, respectively. 

Observed water-level altitudes in the deep wells during 

1995–2005 were 1.5–1.9 m greater than those in shallow 

wells, and simulated steady-state water-level altitudes in 

layers 2 and 4 were 0.7 and 1.1 m greater than those in layer 

1, respectively. Observed gradients near Cottonwood were 

derived from piezometers that were screened at depths of 5 m 

and 60 m within the pre-entrenchment alluvium and basin 

fill, respectively. Observed deep water-level altitudes near 

Cottonwood were 6 m greater than those in shallow wells, 

and simulated steady-state water-level altitudes in layers 2 

and 4 were about 4.2 and 6.1 m greater than those in layer 1, 

respectively. Observed gradients near Palominas were derived 

from piezometers that were screened at depths of 23 m and 

60 m within the basin fill. Observed water-level altitude in 

the deep well near Palominas was more than 2 m greater than 

the water-level altitude in the shallow well during 2001–03, 

and simulated steady-state water-level altitudes in layers 

2 and 4 were 0.6 and 1.7 m greater than those in layer 1, 

respectively.

32  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
W

ATER-LEVEL ALTITUDE, IN METERS

Figure 9. Pre-1961 water-level hydrographs for the Sierra Vista subwatershed of the Upper San Pedro basin, United States.

1,270

1,255

1,265

1,260

1,275

1,290

1,285

1,280

1,295

1,450

1,455

1,460

1,465

1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961

(D–23–22)33dcd1, 2

(D–24–23)06aaa1

(D–21–21)22ddc

(D–21–21)27cbd

(D–23–22)21cda

(D–21–21)29cca

(D–23–21)06ccc2

(D–24–23)06aaa1

(D–23–22)18bbb

31o30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

ARIZONA

SONORA

UNITED STATES

MEXICO

Bisbee

Tombstone

110o30’

110o

15’

P de ro R vi er

San

Sierra Vista

R. 17 E.

R. 18 E.

R. 19 E. R. 20 E. R. 21 E. R. 22 E. R. 23 E. R. 24 E.

T.

21

N.

T.

20

N.

T.

19

N.

T.

18

N.

T.

22

N.

T.

23

N.

T.

24

N.

(D–21–21)27cbd

(D–21–21)29cca

(D–21–21)22ddc

(D–23–22)21cda

(D–23–22)33dcd1, 2

(D–23–22)18bbb

(D–23–21)06ccc2

Figure 9. Pre-1961 water-level hydrographs for the Sierra Vista subwatershed of the Upper San Pedro Basin, United States.

Simulation of Ground-Water Flow  33
Figure 10. Observed and simulated predevelopment water-level altitude error for selected monitoring wells, Upper San Pedro Basin, 

United States and Mexico. 10A, Observed and simulated predevelopment water-level altitudes; 10B, Residual of observed and simulated 

predevelopment water-level altitudes.

34  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico

Figure 10. Observed and simulated water-level altitude error for selected monitoring wells, Upper San Pedro Basin, United States and 

Mexico. A, Observed and simulated predevelopment water-level altitudes; B, Residual of observed and simulated predevelopment 

water-level altitudes.

B. Residual of observed and simulated predevelopment water-level altitudes RESIDUAL OF OBSERVED AND SIMULATED WATER-LEVEL ALTITUDE, IN METERS

40

30

20

10

0

–10

–20

–30

–40

1,150 1,200 1,250 1,300 1,350 1,400

OBSERVED WATER-LEVEL ALTITUDE, IN METERS

1,450 1,500 1,550 1,600

A. Observed and simulated predevelopment water-level altitudes SIMULATED 

WATER-LEVEL ALTITUDE, IN METERS

1,550

1,500

1,450

1,400

1,350

1,300

1,250

1,200

1,150

1,150 1,200 1,250 1,300 1,350 1,400

OBSERVED WATER-LEVEL ALTITUDE, IN METERS

NOTE: Residual of observed hydraulic-head difference during

 2000 and simulated predevelopment difference

1,450 1,500 1,550 1,600

NO ERROR

PRE-1961 DATA

PRE-1950 DATA

VERTICAL HYDRAULIC-HEAD

 DIFFERENCE

EXPLANATION
Oscillatory Steady-State Conditions

An oscillatory steady-state simulation was completed 

to simulate predevelopment seasonal evapotranspiration 

conditions and to establish initial conditions for transient 

simulations of seasonal pumping and evapotranspiration. 

Simulation of average-annual steady-state conditions 

approximated the oscillatory steady-state conditions, but 

did not include the effects of seasonal evapotranspiration on 

streamflow and storage change that would be expected during 

the initiation of the seasonal transient simulation. Lack of 

inclusion of these seasonal initial conditions could result in 

errors in the transient simulation.

The oscillatory steady-state simulation included seasonal 

evapotranspiration stresses—fall/winter and spring/summer—

for a 30-year period by using initial conditions from the 

steady-state simulation. Water levels near the streams attained 

oscillatory steady-state conditions after 30 years of simulation, 

but water levels continued to change at low rates across the 

regional aquifer. 

The oscillatory steady-state simulation resulted in 

hydraulic head distributions that were similar to those of 

the average-annual steady-state simulation, but the averageannual water budget was slightly different (table 4). 

Recharge rates were identical to the average-annual 

steady-state rate of 58,100 m3/d. Seasonal variation in 

evapotranspiration resulted in seasonal variations in 

streamflow and storage. Inflow to the ground-water system 

from streams ranged from 24,000 m3/d in fall/winter to 

29,000 m3/d in spring/summer and averaged 27,000 m3/d 

annually. Inflow from streams represents streamflow that 

reinfiltrates into the ground-water flow system after 

previously discharging to the stream in upstream reaches 

and was not inflow to the ground-water system that was 

in addition to recharge. Average seasonal outflow to 

streams ranged from 54,800 m3/d in spring/summer to 

62,400 m3/d in fall/winter and averaged 57,900 m3/d annually. 

Seasonal evapotranspiration rates ranged from an average 

of 36,500 m3/d during spring/summer to 0 m3/d during 

fall/winter and averaged 21,700 m3/d annually. Underflow to 

the adjacent downgradient basin was 2,700 m3/d during both 

spring/summer and fall/winter seasons. Discharge to drains 

was 1,500 m3/d during spring/summer and fall/winter. Water 

was removed from storage (inflow to the ground-water flow 

system) at an average seasonal rate of 10,000 m3/d during 

spring/summer and 300 m3/d during fall/winter. Water flowed 

into storage (outflow from the ground-water flow system) at 

average seasonal rates of 16,100 and 1,500 m3/d during the 

fall/winter and spring/summer, respectively. 

The oscillatory simulation of seasonal predevelopment 

conditions resulted in seasonal baseflow variation at 

Charleston and Palominas that can be compared with observed 

baseflow during the late 1930s. Minimum seasonal baseflows 

at Charleston were simulated as 26,800 m3/d (10.9 ft3/s) and 

12,600 m3/d (5.2 ft3/s) for fall/winter and spring/summer, 

respectively. Observed seasonal baseflows at Charleston were 

40,300 m3/d (16.5 ft3/s) and 9,000 m3/d (3.7 ft3/s) in the late 

winter and early summer (June), respectively. A portion of the 

greater observed seasonal variation may be related to groundwater withdrawals near Palominas that may have removed 

9,000 m3/d (3.7 ft3/s) from the ground-water system or directly 

from the stream during the spring/summer. A portion of 

observed baseflow also was supplied by near-stream recharge 

during flood events and was not simulated by the model. 

Minimum seasonal baseflows at Palominas were simulated as 

12,300 m3/d (5.0 ft3/s) and 8,300 m3/d (3.4 ft3/s) for fall/

winter and spring/summer, respectively. Observed seasonal 

baseflows at Palominas were 11,700 m3/d (4.8 ft3/s) and 

1,700 m3/d (0.7 ft3/s) for late winter and early summer (June), 

respectively.

Transient-State Conditions

Transient conditions resulting from ground-water 

withdrawals and variations in evapotranspiration rates 

were simulated for 1902–2002. Water-level hydrographs 

representing several hydrogeologic conditions across the basin 

were used to calibrate transient conditions by adjusting storage 

and vertical hydraulic properties to approximate observed 

trends (fig. 11 and 12). Observed seasonal and long-term low 

flows at Palominas and Charleston and simulated seasonal 

baseflow (fig. 13) were compared. Long-term water-level 

hydrograph records are few. Records of two to four decades 

in length are, however, available for several wells (fig. 11). 

Primary water-level hydrographs for the Palominas-Hereford 

agricultural area include a long-term record for well 

(D-23-22)33dcd1,2 and records of monthly values for 

2000–2002 for wells (D-23-22)17acc, PAL-SUM, and PALUWD. Other long-term records near the agricultural area 

include those for wells (D-23-22)18bbb and (D-24-23)06aaa1. 

Primary hydrographs for the Sierra Vista area include pre1980 records for well (D-21-21)29cca, post-1980 records for 

well (D-22-22)06dac, bi-monthly records beginning in 1995 

for several observation wells at Fort Huachuca, and continuous 

records for a well cluster near the San Pedro River at Highway 

90 [(D-22-22)06aaa1,2,3,4,5 and (D-22-22)06abd, also 

known as Lewis Springs MW1–5 and MW6, respectively]. 

Fort Huachuca water-level data compared with simulated 

water-level trends include MW1, (D-21-21)10ada; MW3, 

(D-21-21)19ddc; MW4, (D-21-20)15acd; MW5, (D-21-

21)16aaa; MW7, (D-20-21)33ddc; TW4, (D-20-20)35ccb; and 

TW6, (D-20-20)35ccb. Hydrographs for the Sonoran portion 

of the model area include data at several wells beginning about 

1980 (Consultores en Agua Subterranea S.A. por Mexicana 

de Cananea, S.A. de C.V., 2000), but the well locations are 

not known. 

Simulation of Ground-Water Flow  35
Table 4. Simulated predevelopment ground-water budget.

 [cubic-meters per day]

Water-budget 

component

Steady-state simulation Oscillatory steady-state simulation

Average annual Average spring/summer Average fall/winter Average annual

Inflow Outflow Inflow Outflow Inflow Outflow Inflow Outflow

Recharge 58,100 0 58,100 0 58,100 0 58,100 0

Stream baseflow1 27,000 58,700 29,000 54,800 24,000 62,400 27,000 57,900

Evapotranspiration 0 22,000 0 36,500 0 0 0 21,700

Ground-water under flow 0 2,700 0 2,700 0 2,700 0 2,700

Drains2 0 1,500 0 1,500 0 1,500 0 1,500

Subtotal3 85,100 84,900 87,100 95,500 82,100 66,600 85,100 83,800

Storage4 0 0 10,000 1,500 300 16,100 6,100 7,400

TOTAL5 85,100 84,900 97,100 97,000 82,400 82,700 91,200 91,200

1

Inflow from stream baseflow represents stream baseflow that reinfiltrated the ground-water system after discharging to the stream and is not additional 

inflow to the ground-water system. 

2

Drains include discharge to several springs in the Huachuca Mountains.

3

Subtotal of ground-water flow components.

4

Storage inflow represents removal from storage and input to the ground-water flow system. Storage outflow represents removal from the ground-water flow 

system and input to storage.

5

Total flow is the sum of all flow and storage components.

Observed water-level trends are well simulated for most 

hydrograph comparison sites (fig. 12). Simulated waterlevel records include data for the model layer that represents 

water-table conditions—layer 1, 2, or 4—and the model 

layer that represents the deeper confined aquifer—usually 

layer 4. Trends that were similar but showed no significant 

differences in hydraulic head between the water-table and 

confined aquifers were simulated for several wells in regions 

of thick silt and clay, including Fort Huachuca wells MW4, 

TW6, MW4, MW5, (D-22-22)06dac, (D-23-22)33dcd1,2, 

and (D-23-22)17acc; piezometers near Lewis Springs (Lewis 

Springs MW1 through 6); and piezometers near Palominas, 

PAL-SUM, PAL-UWD, and PAL-UWS. Water-level declines 

of 2–4 m observed before 1970 in wells (D-21-21)29cca, 

(D-23-22)18bbb, and (D-24-23)06aaa1 were slightly greater 

than simulated water levels. Differences in the water-level 

decline observed before 1970 and simulated water-level 

decline may be caused by no simulation of low recharge 

rates during the mid-century drought from 1945 to 1960, 

not simulating possible changes in evapotranspiration rates 

and changes in stream elevation, or simulating the spatial 

distribution of withdrawals inaccurately. Simulated declines 

in the confined portions of the aquifer near Palominas were 

greater than the observed water-level decline at well 

(D-23-22)33dcd1,2. This well was open to both the basin fill 

and the stream alluvium, and water levels in the well likely 

are more representative of the higher hydraulic head in the 

water-table aquifer within the stream alluvium than the lower 

hydraulic head in the basin fill.

The accuracy of the simulation of post-1990 water-level 

trends was an important measure of model calibration because 

distributions of seasonal ground-water withdrawals are well 

known for the period. Withdrawals are best defined during 

2001–03 on the basis of recent surveys that documented 

withdrawal rates and locations (Arizona Department of Water 

Resources, 2005). Water levels during the 2001–03 period 

were monitored at several observation wells in the Sierra Vista 

and Palominas areas, including deep and shallow piezometers 

near Lewis Springs (Lewis Springs MW1 through 6) and 

Palominas (wells PAL-SUM and PAL-UWD; fig. 11). 

Simulated water-level trends match observed trends 

during 1995–2003 at several wells in the Sierra Vista area. 

At wells that are believed to represent water-table conditions, 

Fort Huachuca wells TW4, MW1, MW3, MW4, and MW5, 

simulated water-level altitudes in the shallowest active layer, 

normally layer 2, are within about 6 m of observed altitudes, 

and simulated trends are similar to observed trends (fig. 12). 

Wells that have screened intervals beneath significant 

thicknesses of silt and clay, Fort Huachuca wells TW6, 

MW7, and(D-22-22)06dac, are believed to represent confined 

conditions at wells. Simulated water-level altitudes in layer 

4 at these wells are within about 6 m of those observed, and 

simulated trends are similar to those observed. 

Downward ground-water flow was simulated among 

layers in the region of thick silt and clay between the 

San Pedro River and Sierra Vista. Significantly lower water-

36  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
level altitudes of 5–18 m were simulated in layer 4 with 

respect to layer 2 for wells TW6, MW4, and MW5. Upward 

ground-water flow was simulated for well (D-22-22)06DAC 

and Lewis Springs piezometers MW1–5 with water-level 

altitudes in layer 4 that were about 2–3 m greater than the 

water-level altitude in layer 1. Agricultural withdrawals from 

layer 4 at several nearby wells during 1973–86 resulted in 

strong seasonal water-level variations of about 1–2.5 m at 

(D-22-22)06dac and in the deep piezometers at Lewis Springs, 

Lewis Springs MW5 and MW6. Vertical directions of groundwater flow in the Sierra Vista area did not change during the 

simulation, but downward gradients increased for TW6, MW4, 

and MW5 after about 1950, and upward gradients increased 

for well (D-22-22)06dac and Lewis Springs piezometers 

MW1–6 following the cessation of nearby agricultural 

withdrawals after 1986.

Water-level trends and seasonal variations in water levels 

in the Palominas area during 2001–03 were well simulated for 

wells (D-23-22)33dcd1,2, (D-23-22)17acc, PAL-SUM, PALUWD, and PAL-UWS. The greater seasonal withdrawal during 

spring/summer simulated in layers 2 and 4 was caused 

by nearby agricultural withdrawals. Only well (D-23-

22)33dcd1,2 had included a long-term water-level record, 

which began about 1950. A slightly declining simulated 

water-level trend in the unconfined aquifer before 1980, 

layer 1 at (D-23-22)33dcd1,2, was similar to the observed 

trend; but simulated water-level altitudes were about 3–4 m 

below observed values. Simulation of the declining waterlevel trend before 1970 in layer 2 resulted in a significant 

reduction in upward hydraulic gradients toward layer 1 

from about 3 m in 1940 to less than 1 m after about 1970. 

Stronger upward hydraulic gradients were reestablished 

after 1986 when withdrawals in nearby agricultural areas 

were reduced. Simulated seasonal water-level variations 

of about 2.5 m in layer 4 at well (D-23-22)17acc during 

2001–03 were slightly greater than observed variations of 

1.5–2 m. Well (D-23-22)17acc is more than 400 m deep. 

No lithologic log exists, but the well likely penetrates a thick 

sequence of silt and clay, and water levels likely represent 

confined conditions. Simulated water-level altitude, upward 

vertical gradients, and seasonal water-level variations in 

layers 1, 2, and 4 at shallow and deep piezometers PALUWS, PAL-UWD, and PAL-SUM, respectively, are similar to 

observed conditions during 2001–03. Piezometer PAL-SUM 

is screened across about 3 m about 60 m beneath a thick 

sequence of silt and clay, and water levels likely represent 

confined conditions. Piezometer PAL-UWD is screened 

across about 3 m about 20 m beneath an interbedded sequence 

of clay, silt, sand, and gravel, and water levels likely represent 

confined conditions. Piezometer PAL-UWS is screened 

across about 1.5 m at a depth of about 9 m within sand and 

gravel, and water levels likely represent water-table conditions.

Observed summer and winter stream baseflow at 

the Charleston and Palominas streamflow-gaging stations 

are important calibration constraints. Discharge from the 

regional ground-water flow system to the stream—baseflow—

was approximated by using the seasonal 7-day minimum 

flow at each station. Seven-day minimum flows, however, 

normally include a surface-water component, especially 

following periods of high runoff. Seven-day minimum 

flows were, therefore, considered as high estimates of groundwater discharge except for periods of extended drought. 

Simulated stream baseflow at each station was approximated 

by using rates of streamflow during the time steps of each 

seasonal stress period. The magnitude of summer baseflow 

was well simulated for both stations. Simulated winter 

baseflow at Palominas was similar to that observed; 

simulated winter baseflow at Charleston was less than that 

observed (fig. 13).

Simulated rates of summer and winter baseflow at 

Charleston are similar to estimated rates during 1936–40 

(fig. 13A). Simulated minimum summer baseflow for the 

period was about 13,900 m3/d (5.7 ft3/s). The observed 

average 7-day minimum low flow was about 9,000 m3/d 

(3.7 ft3/s). Simulated winter baseflow during 1936–40, 

28,500 m3/d (10.9 ft3/s), was less than the minimum 7-day 

low flow, 30,400 m3/d (12.4 ft3/s). Simulated and estimated 

baseflow decline followed the initiation of significant nearstream ground-water withdrawals for agricultural use in the 

Palominas-Hereford area and growth for riparian woodland 

vegetation after 1940. Simulated and estimated minimum 

summer baseflows were similar in magnitude and generally 

declined throughout the simulation period; both baseflows 

declined from about 8,800 m3/d (3.6 ft3/s) during the early 

1940s to about 2,800 m3/d (1.1 ft3/s) near the end of the 

simulated period. Overall declines in rates of the maximum 

recovery of baseflow during the winter at Charleston were 

estimated and simulated as about 17,000 m3/d (7.0 ft3/s) and 

18,500 m3/d (7.6 ft3/s), respectively. Much of the simulated 

decline in winter baseflow, however, occurred during 1940–46. 

Similar estimated declines occurred during 1940–50, but 

estimated values were highly variable after 1950. Much of 

the variability was likely caused by including significant 

amounts of runoff in the estimated 7-day minimum seasonal 

flow during and following extended wet periods. Declines 

in the lowest estimated winter baseflow values that occurred 

during periods of drought in 1938, 1950, 1952, and 1997 

approximately matched the simulated declines.

Simulated rates of winter baseflow at Palominas during 

1936–40, 12,100 m3/d (4.9 ft3/s), are similar to estimated 

rates; however, simulated summer rates, 8,100 m3/d (3.3 ft3/s), 

are greater than observed rates of about 1,700 m3/d (0.7 ft3/s; 

fig.13B). Simulated and estimated baseflow decline followed 

the initiation of significant near-stream ground-water 

withdrawals for agricultural use in the Palominas-Hereford 

area and growth of riparian woodland vegetation after 1940.

Simulation of Ground-Water Flow  37
30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

UNITED STATES

MEXICO

Sierra Vista

Cananea

Tombstone

110o30’

15’

15’

31o

110o

TERSH D BOUNDARY

E

WA

P de

ro

R

vi

er

San

Geology modified from Reynolds, 1988

EXPLANATION

Figure 11. Distribution of simulated withdrawal wells and selected hydrograph comparison sites, Upper San Pedro Basin,

United States and Mexico.

MODEL GRID ROW 1

MODEL GRID ROW 440

MODEL GRID COLUMN 1

MODEL GRID COLUMN 320

TW6

MW7

MW1 Charleston

Palominas

MW5

MW6

MW3 MW4

MW3, 5

TW4

(D–21–21)29cca

(D–22–22)06dac

(D–23–22)33dcd1, 2

(D–23–22)18bbb

(D–23–22)17acc

(D–24–23)06aaa1

PAL-SUM

PAL-UWD and PAL-UWS

OBSERVED AND SIMULATED WATER-

 LEVEL HYDROGRAPH

ESTIMATED AND SIMULATED STREAM

 BASEFLOW HYDROGRAPH

AGRICULTURAL WITHDRAWAL WELL

PUBLIC-USE WITHDRAWAL WELL

FORT HUACHUCA WITHDRAWAL WELL

OTHER WITHDRAWAL WELL

Lewis Springs

Bisbee

Extent of model grid

Figure 11. Distributions of simulated withdrawal wells and selected hydrograph comparison sites, Upper San Pedro Basin, 

United States and Mexico.

38  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
W

ATER-LEVEL ALTITUDE, IN METERS

1,205

1,210

1,215

1940 1950 1960 1970 1980 1990 2000 2010

1,220

1,225

1,215

1,230

1940 1950 1960 1970 1980 1990 2000 2010

1,255

1,260

1,265

1,250

1,270

1,275

W

ATER-LEVEL ALTITUDE, IN METERS

1,260

1,255

1,265

1,270

1,275

1,260

1,265

1,270

1,275

1,265

1,260

1,270

1,275

1,280

1,260

Well (D-21-21)29cca

Well Fort Huachuca TW4

Well Fort Huachuca MW4

Well Fort Huachuca TW6

Well Fort Huachuca MW3

Well Fort Huachuca MW5

Well Fort Huachuca MW1

Well Fort Huachuca MW7

1,265

1,270

1,275

1,240

1,245

1,250

1,235

1,230

1,255

1,260

Figure 12. Observed and simulated water-level altitudes at selected wells, 1940–2005, Upper San Pedro basin, United States and Mexico. Figure 12. Hydrographs showing observed and simulated water-level altitudes at selected wells, 1940–2003, Upper San Pedro Basin, 

United States and Mexico.

Simulation of Ground-Water Flow  39
W

ATER-LEVEL ALTITUDE, IN METERS

1940 1960 1970 1980 1990 2000 2010 1,265

1,270

1,260

1,255

1,250

1,275

1940 1950 1960 1970 1980 1990 2000 2010

1,275

1,280

1,290

1,295

1,300

W

ATER-LEVEL ALTITUDE, IN METERS

1,225

1,230

1,235

1,240

1,285

1,290

1,295

1,300

1,280

1,285

1,275 1950

1,230

1,235

1,240

1,245

1,225

Well (D-22-22)06dac

Well Lewis Springs MW6 (55 meters deep)

Well (D-24-23)06aaa1

Wells Lewis Springs MW5 and MW3,

MW5 (55 meters deep),

MW3 (7 meters deep)

Well (D-23-22)18bbb

Well (D-23-22)33dcd1, 2

Well (D-23-22)17acc

Well near Palominas

PAL-SUM (60 meters deep),

PAL-UWD (20 meters deep), and

PAL-UWS (9 meters deep)

1,230

1,235

1,240

1,275

1,280

1,290

1,285

1,290

Figure 12. Continued.

SIMULATED L1

SIMULATED L2

SIMULATED L4

EXPLANATION

OBSERVED VALUE (SYMBOL OR 

 CURVE)—Color of curve varies

 for graphs with multiple piezo-

 meters. Different color of curves

 are explained for individual graphs.

Figure 12. Continued.

40 Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Figure 13. Estimated and simulated seasonal baseflow at streamflow-gaging stations on the San Pedro River, Arizona, 1930–2005: A, 

Charleston; B, Palominas.

1935 1940 1945 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2005 2010

B. Palominas

0

1930 1935 1940 1945 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 2010

2000

30,000

20,000

10,000

40,000

50,000

A. Charleston

0

1930 1935 1940 1945 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2005 2010

SEASONAL LOW FLOW, IN CUBIC METERS PER D

AY

30,000

20,000

10,000

40,000

50,000

OBSERVED LATE-WINTER 7-DAY LOW FLOW—Flood-

 dominated flow eliminated

OBSERVED EARLY-SUMMER (JUNE) 7-DAY LOW FLOW

SIMULATED SEASONAL BASEFLOW

EXPLANATION

SEASONAL LOW FLOW, IN CUBIC METERS PER D

AY

Figure 13. Estimated and simulated seasonal baseflow at streamflow-gaging stations on the San Pedro River, Arizona, 1929–2003. 

13A, Charleston; 13B, Palominas.

Simulation of Ground-Water Flow  41
The simulated decline in the maximum baseflow recovery 

during winter of about 9,300 m3/d (3.8 ft3/s) occurred more 

gradually than the estimated decline in the maximum winter 

7-day low flow, which mostly occurred during 1940–50 when 

the streamflow-gaging station was not operating. Estimated 

summer baseflow declined to near zero during 1940–50. 

Simulated summer baseflow declined to zero more gradually 

during 1940 to about 1970. Differences in the simulated 

and estimated baseflow trends at Palominas may have been 

caused by inaccurately simulating the poorly known spatial 

distribution of ground-water withdrawals, or by simulated 

vertical hydraulic properties that are lower than in the real 

system.

A large amount of water-level data was available for 

comparison of observed and simulated water-level altitudes 

during the winter of 2002 (fig. 14). Observed and simulated 

water-level altitudes generally are within about 10 m of each 

other throughout the regional alluvial aquifer. Differences 

are greater than 10 m in the higher altitudes of Walnut Gulch 

where simulated values were much less than those observed. 

Simulated water-level altitude in the Huachuca Mountains also 

was less than those observed.

Rates of outflow from the ground-water flow system 

during the final year simulated (March 2002–March 2003) 

were greater than rates of inflows (table 5, subtotal), and the 

difference is equivalent to net change in ground-water storage. 

Inflow from recharge through natural sources was equivalent 

to the steady-state rate of 58,000 m3/d (17,200 acre-ft/yr). 

Additional recharge through artificial sources included 

6,500 m3/d (1,900 acre-ft/yr) at turf and recharge facilities and 

7,200 m3/d (2,100 acre-ft/yr) in agricultural areas and at septic 

systems that were simulated at recharge wells. Discharge to 

wells, evapotranspiration, drains, and ground-water underflow 

to the adjacent basin averaged 171,700 m3/d (50,800 acreft/yr). Annual ground-water withdrawals were 140,000 m3

/d 

(41,400 acre-ft/yr). A portion of annual withdrawals was 

returned to the ground-water flow system as artificial 

recharge at several municipal facilities, irrigated areas, and 

septic systems. Artificial recharge at municipal facilities was 

simulated as 6,500 m3/d (1,900 acre-ft/yr). Artificial recharge 

in irrigated areas and at septic systems was simulated as 

7,200 m3/d (2,100 acre-ft/yr). Net ground-water withdrawals, 

after accounting for artificial recharge, were 126,300 m3/d 

(37,400 acre-ft/yr), of which 73,800 m3/d (21,800 acre-ft/yr) 

were in Mexico and 52,500 m3/d (15,500 acre-ft/yr) were in 

the United States. Discharge from the ground-water system 

to stream baseflow was about 26,600 m3/d (7,900 acre-ft/yr), 

and inflow to the ground-water system from stream baseflow 

was about 17,000 m3/d (5,000 acre-ft/yr). Discharge through 

evapotranspiration was simulated as about 27,600 m3/d 

(8,200 acre-ft/yr), about 2,200 m3/d (700 acre-ft/yr) of which 

occurred in Mexico. The simulated rate of evapotranspiration 

in the Sierra Vista subwatershed was about 25,400 m3/d 

(7,500 acre-ft/yr), which was near the estimated rate of 

25,000–30,000 m3/d (7,300-9,000 acre-ft/yr; Leenhouts and 

others, 2005). Simulated discharge through ground-water 

underflow across the northern boundary was 2,600 m3/d 

(800 acre-ft/yr), nearly the same as the estimated rate, which 

indicates that the simulated ground-water flow system near the 

outflow boundary was minimally affected by transient changes 

in the system. 

The imbalance between ground-water outflows and 

inflows was offset by a decrease in ground-water storage 

(table 5). Net removal from ground-water storage (inflow from 

storage minus outflow from storage) averaged 109,400 m3/d 

(32,400 acre-ft/yr) for the final year of the simulation. 

The storage deficit was distributed among the Sierra Vista 

subwatershed and Mexico portions of the model area as 

about 38,300 and 71,100 m3/d (11,300 and 21,100 acre-ft/yr), 

respectively.

Model Limitations and Possible Improvements

The numerical ground-water flow model described 

here simulates an improved conceptualization of the groundwater system in the Upper San Pedro Basin; however, 

further improvements can be implemented. Additional 

improvements could be made in the simulation of boundary 

conditions, streamflow, variations in recharge rates, and 

evapotranspiration. Improved knowledge of the spatial 

distribution of ground-water withdrawals is needed for the 

unincorporated areas of the United States and in Mexico. 

Improved knowledge of the vertical distribution of groundwater withdrawals also could improve simulations, especially 

for confined parts of the aquifer system, such as in the 

Palominas-Hereford area. 

All numerical ground-water flow models are limited 

by boundary conditions. Models will produce invalid results 

when changes within the system result in significant changes 

at incorrectly simulated boundaries. All boundaries in this 

model, except some flow boundaries along the northern extent, 

were assumed to be no-flow. There were several simulated 

areas where ground water could flow to adjacent ground-water 

flow systems; therefore, significant changes would render 

the no-flow assumption invalid. These conditions may occur 

where saturated permeable rocks, such as limestone or alluvial 

deposits, are continuous across the model boundary as in 

several areas, including the Mule Mountains, the Huachuca 

Mountains, and areas near Elgin and along the northern 

and southern model boundaries. These boundaries could be 

modified by simulating a larger extent, or by simulating the 

boundary as a general head boundary. Observed changes near 

boundaries are likely significant in the current simulations 

only where withdrawals for mine use have occurred near the 

Copper Queen and Cananea mines. Withdrawals near the 

Copper Queen mine were not simulated because of a lack of 

understanding of the ground-water flow system in the area. 

These withdrawals could only be simulated by using a greater 

extent in the area. Withdrawals near the Cananea mine were 

simulated, but induced inflow from the adjacent basin was 

assumed to be minimal because withdrawals also occur south 

of the boundary. Errors in the simulation of the regional 

ground-water flow system that were caused by inaccurate 

boundaries near both mines likely are small because of the 

great distance from the nearest areas of natural discharge 

along streams. 

42  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
30’

31o45’

0 10 KILOMETERS

0 6 MILES

Base from U.S. Geological Survey

digital data, 1:100,000, 1982

Universal Transverse Mercator

projection, Zone 12

UNITED STATES

MEXICO

Sierra Vista

Cananea

110o30’

15’

15’

31o

110o

TERSH D BOUNDARY

E

WA

P

de ro

R

vi

er

San

Geology modified from Reynolds, 1988

EXPLANATION

Figure 14. Observed and simulated water-level altitude, winter 2002, Upper San Pedro Basin, United States and Mexico.

MODEL GRID ROW 1

MODEL GRID ROW 440

MODEL GRID COLUMN 1

MODEL GRID COLUMN 320

WELL WITH WINTER 2002 

 WATER-LEVEL DATA

Tombstone

1,60

1,600

0

1,40 0

1,36 0

1,4 2 0

1,3 6 0

1,4 4 0

1,4 4 0

1,5 0 0

1,400

1,380

1,280

1,340

1,380

1,280

1,320 1,280

1,500

1,400

1,440

1,320

1,300

1,420

1,440

1,440

1,460

1,480

1,400

1,

1,380

440

1,420

1,420

1,380

1,34 0

1,400

1,300

1,260

1,220

1,200

1,180

1,160

1,240

1,

1,300

1,340

1,380

1,280

1,320

1,480 1,36 0

1,4 0 0

260

1,200

1,280

1,2 6 0

1,3 60

1,400

1,460

1,7 00

1,600

1,440

1,420

1,

1,

1

1

8

6 0 0

1,220

1,240

1,240

1,500

1,2 80 1,3 00

1,320

1,400

OBSERVED WATER-LEVEL ALTITUDE,

 WINTER 2002—Contour interval 20

 (lighter color) and 100 meters (darker)

1,360 SIMULATED WATER-LEVEL ALTITUDE,

 WINTER 2002—Contour interval 20

 (lighter color) and 100 meters (darker)

1,220

1,400 1,400

Bisbee

Extent of model grid

Figure 14. Observed and simulated water-level altitudes, winter 2002, Upper San Pedro Basin, United States and Mexico.

Simulation of Ground-Water Flow  43
Boundaries along the northern extent of the model were 

simulated as no-flow and constant-head boundaries throughout 

the period of transient stress. The flow system could, however, 

change near these boundaries. In this case, these boundaries 

may be simulated as general-head boundaries, or the model 

domain may be extended to an impermeable boundary. 

Inaccurate simulation of the flow boundaries may be 

significant locally, but should have minimal influence on the 

simulation of the regional ground-water flow system because 

estimated ground-water flow across the flow boundaries is less 

than 5 percent of the simulated water budget. 

 Recent hydrologic investigations and this numerical 

simulation have shown that near-stream recharge during 

floods is an important process that maintains streamflow 

for several months. The streamflow-routing package used 

in this simulation is inadequate for simulation of flood 

flows. More complex and appropriate stream-aquifer 

interaction can be simulated by using a new streamflowrouting package SFR1 (Prudic and others, 2004) that 

would allow simulation of floods and resulting streamflow 

infiltration. Climate-based recharge variations also could 

improve the simulations. Recharge rates in the area have 

been shown to vary significantly on decadal and seasonal 

scales with observed climate variations (Dickinson and 

others, 2004; Pool 2005). Variation in recharge rates likely 

causes some of the observed multiyear variation in stream 

baseflow. Recharge variations can now be estimated by 

using climate data. Simulation of evapotranspiration can 

be improved. Rates of evapotranspiration were simulated 

in this model as a single stress period with constant 

evapotranspiration parameters and a linear depth-to-water 

versus evapotranspiration relation. Evapotranspiration rates, 

however, vary throughout the growing season, and the depth 

versus evapotranspiration relation is nonlinear. Information on 

seasonal evapotranspiration rates and depth-to-water versus 

evapotranspiration is now available for several plant types. 

Improved simulation of streamflow, variations in recharge 

rates, and evapotranspiration should significantly improve 

simulation of seasonal and multiyear ground-water flow.

Summary and Conclusions

A numerical ground-water model was constructed to 

simulate seasonal and long-term variations in ground-water 

flow in portions of the Upper San Pedro Basin in the Sierra 

Vista subwatershed, Ariz., United States, and Sonora, Mexico. 

The simulation used a numerical MODFLOW model for 

1902–2003. Several improvements from previous simulations 

of ground-water flow were incorporated. New concepts of 

ground-water flow that were developed on the basis of data 

Table 5. Simulated ground-water budget for March 2002 to March 2003. 

 [cubic-meters per day]

Water-budget component

Simulated transient ground-water budget, March 2002 to March 2003

Spring/summer Fall/winter Average annual

Inflow Outflow Inflow Outflow Inflow Outflow

Natural recharge 58,000 0 58,000 0 58,000 0

Artificial recharge 6,500 0 6,400 0 6,500 0

Stream baseflow1 17,000 22,400 17,000 32,700 17,000 26,600

Evapotranspiration 0 46,500 0 0 0 27,600

Ground-water underflow 0 2,600 0 2,600 0 2,600

Drains2 0 1,500 0 1,500 0 1,500

Wells3 7,200 165,400 7,200 102,500 7,200 140,000

Subtotal4 88,700 238,400 88,600 139,300 88,700 198,300

Storage5 156,100 6,400 111,800 61,400 138,100 28,700

TOTAL 244,800 244,800 200,400 200,700 226,800 227,000

1

Inflow from stream baseflow represents stream baseflow that reinfiltrated the ground-water system after discharging to the stream and is not additional 

inflow to the ground-water system. 

2

Drains include discharge to several springs in the Huachuca Mountains.

3

Well inflow is from recharge wells that simulate deep percolation and recharge of excess irrigation water and at septic systems.

4

Subtotal of ground-water flow components.

5

Storage inflow includes removal from ground-water storage and input to the ground-water flow system. Storage outflow includes removal from the groundwater flow system and input to storage.

44  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
collected during 1995–2002 were simulated. Ground-water 

flow in the entire watershed was simulated, including flow 

in the sedimentary rocks within the mountains that surround 

and underlie the alluvial basin deposits. Simulation of flow 

in the sedimentary rocks allowed for inclusion of groundwater withdrawals for dewatering purposes at the Tombstone 

mine and discharge to springs in the Huachuca Mountains. 

The model also was constructed to allow for simulation of 

significantly different water levels and vertical ground-water 

flow between deep-confined and shallow-unconfined systems 

separated by thick silt and clay intervals. Recharge was 

distributed throughout the model area rather than only along 

the mountain fronts, which was the paradigm used in previous 

models, on the basis of recent investigations that estimated 

ephemeral channel recharge and new streamflow records 

in upland areas. Variations in discharge from ground-water 

withdrawals for agricultural, public, domestic, and military 

uses were simulated by using two seasons, fall/winter and 

spring/summer, each year.

A 5-layer model was constructed by using water-level 

data, geologic mapping, and geophysical surveys. Boundary 

conditions included no-flow boundaries that mostly coincided 

with the watershed boundaries, which were no-flow at a depth 

of 1,500 m below land surface throughout the model area, and 

specified-head conditions where ground water discharges to 

the north. The upper boundary of the flow system includes 

water-table conditions, streams, evapotranspiration, and drains 

at springs. The deepest layer, layer 5, extends throughout 

the model area and represents the sedimentary, volcanic, and 

crystalline rocks that outcrop in the mountains and underlie 

the alluvial basin. Layer 4 represents basin fill throughout the 

Mexico portion of the basin and sand and gravel in the lower 

basin fill unit in the Sierra Vista subwatershed. Layers 2 and 

3 include the silt and clay facies of upper and lower basin fill, 

respectively, within the Sierra Vista subwatershed and adjacent 

areas that include sand and gravel and interbedded facies in 

the basin fill. Low permeability regions of silt and clay were 

mapped by using available subsurface electrical surveys and 

lithologic logs. Layer 1 represents stream alluvium along 

the San Pedro and Babocomari Rivers in the Sierra Vista 

subwatershed and saturated sand and gravel facies of basin 

fill that overlie the silt and clay facies of upper basin fill. 

Variations in hydraulic properties were distributed throughout 

each layer with polygons representing different rock types, 

facies within the basin fill, stream alluvium, and subareas 

within parts of each rock type.

Steady-state and transient conditions were simulated. 

No true steady-state initial conditions exist according to 

available records that indicate the ground-water flow system 

has changed in response to changes in discharge through 

well withdrawals and evapotranspiration rates and changes 

in stream-channel morphology, including channel incision 

and widening of the floodplain. Initial conditions for both 

the steady-state and transient simulation were determined 

by using the earliest streamflow and water-level records, 

which began in about 1930 and 1940, respectively. Changes 

in the ground-water flow system, following initial groundwater withdrawals for mining use in 1902 and before early 

water levels were measured, were assumed to be minimal. 

Estimated evapotranspiration rates and observed steadystate streamflow and water levels were used as controls to 

calibrate evapotranspiration rates and aquifer-system hydraulic 

properties, hydraulic conductivity and streambed conductance. 

An oscillatory steady-state simulation that included seasonal 

variations in evapotranspiration was done to establish initial 

conditions for simulations of transient conditions during 

1902–2003. Calibration of vertical hydraulic conductivity and 

storage properties was accomplished by matching simulated 

and observed vertical hydraulic gradients and changes in water 

levels and streamflow.

The ground-water flow model approximates observed 

trends in water levels throughout most of the model area 

and streamflow at the Charleston streamflow-gaging station 

on the San Pedro River. Differences between observed and 

simulated water-level and streamflow trends occur near 

Palominas, where observed water-level declines were less than 

those simulated. Observed rapid streamflow declines during 

1940–50 at the streamflow-gaging station at Palominas were 

simulated as a gradual declining trend during the simulation 

period. Much of the difference in observed and simulated 

conditions near Palominas may be related to insufficient 

information on agricultural withdrawal rates and the vertical 

distribution of withdrawals among the unconfined and 

confined parts of the ground-water flow system. Estimated 

increases in evapotranspiration rates in the Sierra Vista 

subwatershed could not be approximated because simulated 

water-level decline was greater than that observed in some 

areas near the San Pedro River, including the Palominas area. 

Simulation of near-stream conditions and evapotranspiration 

rates could be improved by simulating available climaterelated variations in recharge rates and recharge from flood 

flow infiltration by using available improved streamflowrouting methods. 

Summary and Conclusions  45
References Cited

Anderson, T.W., Freethey, G.W., and Tucci, Patrick, 1992, 

Geohydrology and water resources of alluvial basins in 

southcentral Arizona and adjacent states: U.S. Geological 

Survey Professional Paper 1406-D, 74 p.

Arizona Department of Water Resources, 1990, Hydrographic 

survey report (HSR) for the San Pedro River watershed, 

Volume 1: Arizona Department of Water Resources, 548 p.

Arizona Department of Water Resources, 2005, Upper San 

Pedro Basin active management area report, March 2005: 

Phoenix, Arizona Department of Water Resources, 146 p., 

appendices A-M.

Barnes, R.L, and Putman, F., 2004, Maps showing 

groundwater conditions in the Upper San Pedro Basin, 

Cochise, Graham, and Santa Cruz Counties, Arizona: 

Phoenix, Arizona Department of Water Resources 

Hydrologic Map Series Report No. 34, 2 sheets.

Brown, S.G., Davidson, E.S., Kister, L.R, and Thomsen, 

B.W., 1966, Water resources of Fort Huachuca Military 

Reservation, southeastern Arizona: U.S. Geological Survey 

Water-Supply Paper 1819–D, p. D1–D57.

Coes, A.L., and Pool, D.R., 2005, Ephemeral-channel 

and basin-floor infiltration in the Sierra Vista 

subwatershed, Arizona: U.S. Geological Survey Open-File 

Report 05–1023, 67 p. 

Condor Consulting, Inc., 2003, Inversion of airborne EM 

data—Fort Huachuca and Sierra Vista areas, Arizona: 

Consultants Report, 7 p., 1 compact disc.

Consultores en Agua Subterranea S.A. por Mexicana 

de Cananea, S.A. de C.V., 2000, Actualizacion del 

estudio geohidrologico de las cuencas del Rio San 

Pedro y norte del Rio Sonora en Cananea, Son., 136 p., 

5 appendices.

Corell, S.W., Putman, F., Lovvik, Daryl, and Corkhill, Frank, 

1996, A groundwater flow model of the Upper San Pedro 

Basin, southeastern Arizona: Phoenix, Arizona Department 

of Water Resources Modeling Report no. 10, 85 p.

Dickinson, J.E., Hanson, R.T., Ferré, T.P.A., and Leake, 

S.A., 2004, Inferring time-varying recharge from inverse 

analysis of long-term water levels: Water Resources 

Research, v. 40, W07403, doi:10.1029/Water Resources 

Research.2003.002650, 15 p.

Esparza, Jose Guillermo De Aguinaga Ruiz, 2002, 

Modelacion geohydrologica del aquifero del Rio San Pedro: 

Hermosillo, Sonora, Mexico, Universidad De Sonora, 

professional thesis in Geology.

Fleming, John, and Pool, D.R., 2002, Geophysical surveys 

for delineation of shallow structure and lithology near 

the San Pedro River, Southeast Arizona: Proceedings of 

the Society of Environment and Engineering Geophysical 

Society, February 2002, accessed March 28, 2007, at 

http://www.eegs.org/sageep/proceedings.cfm

Freethey, G.W., 1982, Hydrologic analysis of the 

upper San Pedro Basin from the Mexico–United 

States international boundary to Fairbank, Arizona: 

U.S. Geological Survey Open-File Report 82–752, 52 p.

Gettings, M.E., and Houser, B.B., 1995, Preliminary results 

of modeling the gravity anomaly field in the Upper 

San Pedro Basin, southeastern Arizona: U.S. Geological 

Survey Open-File Report 95–76, 12 p. 

Goode, Tomas, and Maddock, Thomas III, 2000, Simulation 

of groundwater conditions in the Upper San Pedro 

Basin for the evaluation of alternative futures: Tucson, 

Department of Hydrology and Water Resources 

and University of Arizona Research Laboratory for 

Riparian Studies, HWR no. 00-30, 113 p. 

Goodrich, D.C., Williams, D.G., Unkrich, K.L., Hogan, 

J.F., Scott, R.L., Hultine, K.R., Pool, D., Coes, A.L., and 

Miller, S., 2004, Comparison of methods to estimate 

ephemeral channel recharge, Walnut Gulch, San Pedro 

River Basin, Arizona: in J.F. Hogan, F.M. Phillips, and 

B.R. Scanlon, eds., Groundwater recharge in a desert 

environment, the southwestern United States, Water 

Science and Applications Series: American Geophysical 

Union, v. 9, p. 77–99.

Gray, R.S., 1965, Late Cenozoic sediments in the 

San Pedro Valley near St. David, Arizona: Tucson, 

University of Arizona, Ph.D. dissertation, 198 p.

Halverson, P.H., 1984, An exploratory gravity survey in the 

upper San Pedro Valley, southeastern Arizona: Tucson, 

University of Arizona, masters thesis, 1984, 85 p.

Hanson, R.T, Newhouse, M.W., and Dettinger, M.D., 2004, 

A methodology to assess relations between climate 

variability and variations in hydrologic time-series in the 

southwestern United States: Journal of Hydrology, v. 287, 

issues 1–4, p. 252–269.

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

M.G., 2000a, 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.

46  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Harbaugh, A.W., Banta, E.R., Hill, M.C., and Anderman, E.R., 

2000b, MODFLOW-2000, The U.S. Geological Survey 

modular ground-water model—User guide to observation, 

sensitivity, and parameter-estimation processes and three 

post-processing programs: U.S. Geological Survey 

Open-File Report 00–184, 209 p.

Hereford, Richard, 1993, Entrenchment and widening of the 

upper San Pedro River, Arizona: Geological Society of 

America Special Paper 282, 46 p. 

Hollyday, E.F., 1963, A geohydrologic analysis of mine 

dewatering and water development, Tombstone, Cochise 

County, Arizona: Tucson, University of Arizona, master’s 

thesis, 90 p.

Kepner, W.G., and Edmonds, C.M., 2002, Remote sensing and 

geographic information systems for decision analysis in the 

public resource administration: A case study of 25 years 

of landscape change in a southwestern watershed: U.S. 

Environmental Protection Agency, EPA/600/R-02/039, 31 p.

Leenhouts, J.M., Stromberg, J.C., and Scott, R.L., 2005, 

Hydrologic requirements of and consumptive groundwater use by riparian vegetation along the San Pedro River, 

Arizona: U.S. Geological Survey Scientific Investigations 

Report, 2005–5163, 211 p.

Phelps Dodge Corporation, 1998, Hydrologic assessment 

for the tailing impoundments-CTSA APP Project Area, 

Bisbee, Arizona: Golden, Colo., Consultant report, SAVCI 

Environmental Technologies, LLC, 102 p., 47 figs., 

1 appendix.

Pool, D.R., 2005, Variations in climate and natural recharge 

in southeast Arizona: Water Resources, 41, W11403, 

doi:10.1029/2004WR003255, 24 p.

Pool, D.R., and Coes, A.L., 1999, Hydrogeologic 

investigations of the Sierra Vista subwatershed of the Upper 

San Pedro Basin, Cochise County, southeast Arizona: 

U.S. Geological Survey Water-Resources Investigations 

Report 99−4197, 41 p.

Pool, D.R., and Leenhouts, J.M., 2002, A multiparameter 

approach for measuring flood-induced aquifer- and bankstorage changes along the San Pedro River, Arizona 

[abs.], in Program & Abstracts, AGU fall 2002 meeting, 

December 6-10, 2002, San Francisco, accessed March 28, 

2007, at http://www.agu.org/meetings/fm02/program.shtml

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

simulate stream-aquifer relations using a modular, finitedifference, ground-water flow model: U.S. Geological 

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

Prudic, D.E., Konikow, L.F., and Banta, E.R., 2004, A new 

streamflow-routing (SFR1) package to simulate streamaquifer interaction with MODFLOW-2000: U.S. Geological 

Survey Open-File Report 2004–1042, 104 p.

Reichardt, K.L., Schladweiler, Brenda, and Stelling, J.L., 

1978, An inventory of riparian habitats along the San Pedro 

River: Tucson, University of Arizona, The Applied Remote 

Sensing Program, Office of Arid Lands Studies, 22 p.

Reynolds, S.J., 1988, Geologic map of Arizona: Arizona 

Geological Survey, Map 26, 1 sheet. 

Schwartzman, P.N., 1990, A hydrologic assessment of the 

lower Babocomari watershed, Arizona: University of 

Arizona, master’s thesis, 3 pls., 212 p.

Southwest Ground-Water Consultants, 2004, Water supply 

potential Phelps Dodge Copper Queen Mine: Phoenix, 

Ariz., Consultant report, Southwest Ground-water 

Consultants, Inc., 24 p.

Thomas, Blakemore, and Pool, D.R., 2006, Seasonal 

precipitation and streamflow trends in southeastern Arizona 

and southwestern New Mexico: U.S. Geological Survey 

Professional Paper 2005–1712, 79 p.

University of Arizona Geophysics Field Camp, 2001, 

Geophysical surveys near Sierra Vista, Arizona: 

Laboratory for Advanced Subsurface Imaging (LASI) 

Report LASI-01-01, May 4, 2001, 34 p.

University of Arizona Geophysics Field Camp, 2002, 

Geophysical surveys near Fort Huachuca, Arizona: 

Laboratory for Advanced Subsurface Imaging (LASI) 

Report LASI-02-01, May 4, 2002, 34 p.

University of Arizona Geophysics Field Camp, 2004, 

Geophysical surveys near Sierra Vista, Arizona: 

Laboratory for Advanced Subsurface Imaging (LASI) 

Report LASI-04-01, June 7, 2004, 109 p.

U.S. Army Corps of Engineers (Topographic Engineering 

Center), 2001, Vegetation map of the San Pedro Riparian 

National Conservation Area and Babocomari River: 

Fort Huachuca, Ariz., Final report submitted to U.S. Army 

Garrison, 63 p.

U.S. Department of Defense, 2002, Fort Huachuca 

programmatic biological assessment for ongoing and 

programmed future military operations and activities: 

Fort Huachuca, Ariz., Environmental and Natural Resources 

Division, Directorate of Installation Support, U.S. Army 

Garrison, 468 p.

References Cited  47
Vionnet, L.B., and Maddock, T., 1992, Modeling of groundwater flow and surface water/groundwater interactions 

in the San Pedro River Basin—Part I—Cananea, Mexico 

to Fairbank, Arizona: Tucson, University of Arizona, 

Department of Hydrology and Water Resources, 

HWR No. 92-010, 96 p.

Webb, R.H., and Betancourt, J.L., 1992, Climatic 

variability and flood frequency of the Santa Cruz River, 

Pima County, Arizona: U.S. Geological Survey WaterSupply Paper 2379, 40 p.

48  Ground-Water Flow Model, Sierra Vista and Sonoran Portions of the Upper San Pedro Basin, Arizona and Mexico
Manuscript approved for publication, 2007.

Prepared by the U.S. Geological Survey Reports Section, 

Tucson, Arizona.

USGS Publishing staff

Tracey L. Suzuki, Technical Editor

John Callahan, Illustrator

For more information concerning the research in this report, 

contact the Arizona Water Science Center Director,

U.S. Geological Survey, 520 N. Park Ave., Suite 221, 

Tucson, AZ 85719

http://az.water.usgs.gov
This page left blank intentionally.
Inside Back Cover

BLANK
Pool and Dickinson—GROUND-WATER FLOW MODEL OF THE SIERRA VISTA SUBWATERSHED AND SONORA PORTIONS 

OF THE UPPER SAN PEDRO BASIN, SOUTHEASTERN ARIZONA, UNITED STATES, AND NORTHERN 

SONORA, MEXICO—U.S. Geological Survey Scientific Investigations Report 2006–5228