IV    Report Title

Scientific Investigations Report 2010–5147

U.S. Department of the Interior

U.S. Geological Survey

Simulated Effects of Groundwater Pumping and Artificial

Recharge on Surface-Water Resources and Riparian

Vegetation in the Verde Valley Sub-Basin, Central Arizona

Prepared in Cooperation with The Nature Conservancy
COVER: Verde River Greenway below Deadhorse Ranch State Park, Cottonwood, Arizona. Courtesy of Jeanmarie Haney, The Nature Conservancy.
Simulated Effects of Groundwater Pumping 

and Artificial Recharge on Surface-Water 

Resources and Riparian Vegetation in the 

Verde Valley Sub-Basin, Central Arizona 

Prepared in Cooperation with The Nature Conservancy

By Stanley A. Leake and Donald R. Pool

Scientific Investigations Report 2010–5147

U.S. Department of the Interior

U.S. Geological Survey
Cataloging-in-publication data are on file with the Library of Congress (URL http://www.loc.gov/).

U.S. Department of the Interior

KEN SALAZAR, Secretary

U.S. Geological Survey

Marcia K. McNutt, Director

U.S. Geological Survey, Reston, Virginia 2010

This report and any updates to it are available online at:

http://pubs.usgs.gov/sir/2010/5147/

For more information on the USGS—the Federal source for science about the Earth,

its natural and living resources, natural hazards, and the environment, 

visit http://www.usgs.gov or call 1-888-ASK-USGS (1-888-275-8747)

For an overview of USGS information products, including maps, imagery, and publications, visit 

http://www.usgs.gov/pubprod 

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

Any use of trade, 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 material contained within this report.

Suggested citation:

Leake, S.A., and Pool, D.R., 2010, Simulated effects of groundwater pumping and artificial recharge on surface-water 

resources and riparian vegetation in the Verde Valley sub-basin, Central Arizona: U.S. Geological Survey Scientific 

Investigations Report 2010-5147, 18 p.
iii

Contents

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

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

Background............................................................................................................................................3

Effects of Groundwater Pumping on Connected Surface Water and Riparian 

Vegetation........................................................................................................................3

The Northern Arizona Regional Groundwater-Flow Model..................................................6

Scope ....................................................................................................................................................8

Method.............................................................................................................................................................8

Results ...........................................................................................................................................................10

Applicability and Limitations......................................................................................................................15

Summary........................................................................................................................................................15

Acknowledgments.......................................................................................................................................18

References Cited..........................................................................................................................................18

Figures

1. Verde Valley sub-basin and relation to the extent of the Northern Arizona Regional 

Groundwater-Flow Model. ..........................................................................................................2

2. Groundwater withdrawal rates in the Verde Valley sub-basin from 1950–2006 simulated 

in the model by Pool and others (in press). ..............................................................................3

3. Stages of effects of a pumping well on a nearby stream. (A) A cone of depression has 

developed around the well, but groundwater still discharges to the stream. (B) After 

some time, the cone of depression has deepened and spread to the nearby stream, 

resulting in elimination of groundwater discharge to the stream and further loss of 

water in the stream to the aquifer. (C) When the stream cannot supply the quantity 

of water pumped, the stream may eventually lose all of its water to the aquifer and 

become ephemeral and nearby riparian vegetation also may be affected. Modified 

from Webb and others (2007)......................................................................................................4

4. Sources of water to a pumped well through time...................................................................5

5. Conceptualized relations among hydrogeologic units and Northern Arizona Regional 

Groundwater-Flow Model layers in the Verde Valley sub-basin. Volcanic rocks and the 

lower Supai Group are confining units. Modified from Pool and others (in press). ..........6

6. Active unconfined and confined model cells in the Northern Arizona Regional Groundwater-Flow Model within the Verde Valley sub-basin for model layers 1, 2, and 3...........7

7. Decrease in surface-water flow and evapotranspiration for the example shown in 

figure 4 for a well that pumps 30 acre-feet per year for years 1–5; 

then for years 6–10 the pumping rate continues (case 1), stops (case 2), and doubles 

(case 3). ........................................................................................................................................11

8. Computed reduction in flow of streams and springs and reduction in riparian evapotranspiration as a fraction of pumping rate that would result from pumping groundwater from model layer 1 at a constant rate for 10 years. ........................................................13
iv

9. Computed reduction in flow of streams and springs and reduction in riparian evapotranspiration as a fraction of pumping rate that would result from pumping groundwater from model layer 1 at a constant rate for 50 years. ........................................................14

10. Computed reduction in flow of streams and springs and reduction in riparian evapotranspiration as a fraction of pumping rate that would result from pumping groundwater from model layer 2 at a constant rate for 10 years. ........................................................16

11. Computed reduction in flow of streams and springs and reduction in riparian evapotranspiration as a fraction of pumping rate that would result from pumping groundwater from model layer 2 at a constant rate for 50 years. ........................................................17

Tables

1. Statistics of aquifer properties and saturated thickness simulated in the part of the 

Northern Arizona Regional Groundwater-Flow Model that underlies the Verde Valley 

sub-basin. Values are for active parts of each model layer. .................................................9

2. Example calculations of absolute depletion from depletion fractions for a well that 

pumps 30 acre-feet per year for years 1–5, and then continues pumping at that rate for 

years 6–10 (case 1), shuts off for years 6–0 (case 2), and doubles in pumping rate for 

years 6–10 (case 3).....................................................................................................................12

3. Average depletion after 100 years of pumping expressed as a fraction of the pumping 

rate for features simulated with the Drain, Evapotranspiration, and Stream Packages, 

for pumping locations in layers 1, 2, and 3. ............................................................................15
v

Conversion Factors 

Multiply By To obtain

Length

foot (ft) 0.3048 meter (m)

mile (mi) 1.609 kilometer (km)

Area

square mile (mi2) 2.590 square kilometer (km2) 

Flow rate

acre-foot per year (acre-ft/yr) 1,233 cubic meter per year (m3/yr)

gallon per minute (gal/min) 0.06309 liter per second (L/s)

Hydraulic conductivity

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

Specific storage

per foot (1/ft) 3.281 per meter (1/m)

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

(NAVD 88).”

Horizontal coordinate information is referenced to the “North American Datum of 1927 (NAD 27).”

Altitude, as used in this report, refers to distance above the vertical datum.
Simulated Effects of Groundwater Pumping and Artificial 

Recharge on Surface-Water Resources and Riparian 

Vegetation in the Verde Valley Sub-Basin, Central Arizona

By Stanley A. Leake and Donald R. Pool

Abstract 

In the Verde Valley sub-basin, groundwater use has 

increased in recent decades. Residents and stakeholders in the 

area have established several groups to help in planning for 

sustainability of water and other resources of the area. One 

of the issues of concern is the effect of groundwater pumping in the sub-basin on surface water and on groundwaterdependent riparian vegetation. The Northern Arizona Regional 

Groundwater-Flow Model by Pool and others (in press) is the 

most comprehensive and up-to-date tool available to understand the effects of groundwater pumping in the sub-basin. 

Using a procedure by Leake and others (2008), this model was 

modified and used to calculate effects of groundwater pumping on surface-water flow and evapotranspiration for areas 

in the sub-basin. This report presents results for the upper 

two model layers for pumping durations of 10 and 50 years. 

Results are in the form of maps that indicate the fraction of the 

well pumping rate that can be accounted for as the combined 

effect of reduced surface-water flow and evapotranspiration. 

In general, the highest and most rapid responses to pumping 

were computed to occur near surface-water features simulated 

in the modified model, but results are not uniform along these 

features. The results are intended to indicate general patterns 

of model-computed response over large areas. For site-specific 

projects, improved results may require detailed studies of the 

local hydrologic conditions and a refinement of the modified 

model in the area of interest.

Introduction

The study area includes the “Verde Valley sub-basin,” 

as designated by the Arizona Department of Water Resources 

(http://www.azwater.gov/azdwr/StatewidePlanning/

RuralPrograms/OutsideAMAs_PDFs_for_web/documents/

verde_river.pdf, accessed April 25, 2010). The area lies mostly 

within the transition zone between the Basin and Range 

lowlands to the south and the Colorado Plateau to the north 

(fig. 1). Land-surface altitudes range from about 2,900 ft 

where the Verde River leaves the sub-basin, to about 12,600 ft 

near the northernmost point on the boundary of the sub-basin. 

Perennial flow of surface water occurs in the Verde River and 

reaches of Sycamore Creek, Oak Creek, Wet Beaver Creek, 

and West Clear Creek. Major vegetation types are desert scrub 

at lower altitudes, piñon-juniper woodlands at intermediate 

altitudes, and conifer forests at higher altitudes. Along the 

surface-water channels in the Verde Valley sub-basin, some of 

the more common native species of woody riparian vegetation are Frémont cottonwood, mesquite, sycamore, ash, black 

willow, seepwillow, coyote willow, and Arizona alder (Webb 

and others, 2007). These and other plants that use groundwater 

generally occur in narrow bands, commonly within tens to 

hundreds of feet of stream banks. 

Population estimates for July 1, 2008 for incorporated cities and towns in the sub-basin (http://www.census.gov/popest/

cities/tables/SUB-EST2008-04-04.csv, accessed March 19, 

2010) include 11,599 in Sedona, 11,412 in Cottonwood, 

10,849 in Camp Verde, 4,263 in Clarkdale, and 353 in Jerome. 

These data suggest that average annual population growth 

rates from 2000 to 2008 for these cities and towns range from 

about 1 to more than 3 percent. Unincorporated areas of the 

sub-basin, including communities of Cornville and Rimrock 

and privately owned rural areas, also are populated. Average annual population growth for the entire sub-basin likely 

has been 1.5 percent or higher since 2000. If that growth rate 

were maintained, population would double in about 28 years. 

The larger cities of Flagstaff and Prescott with estimated 2008 

populations of 60,222, and 42,697, respectively, are in areas 

adjacent to the sub-basin.

Although surface water is diverted from the Verde River 

for irrigation, nearly all water for public, domestic, and 

industrial supplies is groundwater (Blasch and others, 2006). 

Groundwater pumping occurs beneath more populated parts 

of the sub-basin in hydrogeologic units that comprise shallow 

alluvium along streams, fluviolacustrine basin fill, sandstone 

and carbonate sediments, and fractured crystalline rocks. 

Decadal estimates of groundwater withdrawals used in the 

groundwater model by Pool and others (in press) suggest a 
2 Simulated Effects of Groundwater Pumping and Artificial Recharge in the Verde Valley Sub-Basin, Central Arizona

Colorado Plateau

Basin and Range Lowlands

Transition Zone

ARIZONA

Extent of

Northern Arizona

Regional GroundwaterFlow Model

Verde

Valley

sub-basin

EXPLANATION UTM

ZONE 12

NAD 27

0 5 10 20 Miles

Stream reach represented in the model:

Reach with flow to or from the aquifer in steady-state model

Reach without flow to or from the aquifer in steady-state model

Major Roads

Verde Valley sub-basin

Spring in Verde Valley simulated with the Drain Package

COTTONWOOD

JEROME

CAMP VERDE

RIMROCK

CORNVILLE

CLARKDALE

SEDONA

FLAGSTAFF

PRESCOTT

400000 450000

38500 3900000

Verde

Oak

Clear

West

Creek

Wet Beaver Creek

Creek

Creek

River

Sycamore

Figure 1. Verde Valley sub-basin and relation to the extent of the Northern Arizona Regional Groundwater-Flow Model.
Introduction 3

greater than ten-fold increase in groundwater use in the subbasin in the latter half of the 20th century (fig. 2). 

With increasing population and water use, sustainability 

of water resources and long-term health of the Verde River is 

a concern of many residents of the sub-basin. Several groups 

of stakeholders are considering strategies for sustainability 

of water supplies in the Verde Valley sub-basin and adjacent 

areas. This report provides analyses that can help water managers and others understand long-term effects of groundwater 

withdrawals on the Verde River and adjacent groundwaterdependent vegetation.

Background

Winter and others (1998) characterize groundwater and 

surface water as being a single resource. Development of 

groundwater affects connected surface-water resources and 

development of surface-water resources affects connected 

groundwater resources. This report addresses general aspects 

of how groundwater pumping affects surface water and 

evapotranspiration, and presents an overview of the Northern 

Arizona Regional Groundwater-Flow Model (Pool and others, 

in press) as a tool for studying interaction of groundwater and 

surface water in the Verde Valley sub-basin. 

Effects of Groundwater Pumping on Connected 

Surface Water and Riparian Vegetation

In an aquifer system such as the one that underlies the 

Verde Valley sub-basin, the principal source of groundwater 

is recharge from precipitation and snowmelt that percolates 

Figure 2. Groundwater withdrawal rates in the Verde Valley sub-basin from 

1950–2006 simulated in the model by Pool and others (in press).

to the water table. In much of the southwestern United States, 

recharge tends to be focused in stream channels where permeable sediments allow water to infiltrate (Stonestrom and 

others, 2007). Water moves through the aquifer system and 

eventually discharges at lower elevations to springs, streams, 

rivers, and by evapotranspiration in areas of vegetation that 

use groundwater. The general term “evapotranspiration” 

includes evaporation and transpiration of water in shallow soil, 

as well as from the saturated groundwater system. In the context of this report, however, “evapotranspiration” refers to use 

of water from the saturated groundwater system by vegetation.

The potential energy in an aquifer system is characterized by the groundwater level, sometimes referred to as the 

“head.” In an unconfined aquifer system, the potential energy 

is characterized by the elevation of the water table. The slope 

of the potential-energy surface at any point is the groundwater 

“gradient.” Groundwater levels are determined through measurements of water levels in wells, and groundwater gradients 

and directions of movement are estimated by using measurements of water levels from a number of wells in different parts 

of the aquifer system. Groundwater moves down-gradient 

from areas of higher water levels to areas of lower water 

levels. The presence of gradients indicates that groundwater is 

moving from areas of recharge or inflow to areas of discharge. 

The rates of groundwater movement vary greatly among aquifer systems. Under natural or “predevelopment” conditions, 

rates of recharge to and discharge from an aquifer system 

vary seasonally and interannually with varying climate, but 

when viewed as long-term averages, commonly are assumed 

to be equal. The equilibrium condition of equal recharge and 

discharge with no change in groundwater levels is referred to 

as “steady state.”

GROUNDWATER WITHDRAWALS, IN

THOUSANDS OF ACRE FEET PER YEAR

1950 1960 1970 1980 1990 2000 2010

0

2

4

6

8

10

12

14

16

18

20
4 Simulated Effects of Groundwater Pumping and Artificial Recharge in the Verde Valley Sub-Basin, Central Arizona

Removal of groundwater by wells initiates changes to 

an aquifer system with effects that are variable through time. 

When water is initially pumped from a well, all of the water 

comes from storage around the well. A cone of depression 

develops around the well and the gradient of the groundwater 

head is the driving force for movement of water into the well. 

If the aquifer is unconfined, nearly all of the storage change is 

from draining of pore spaces at the water table. This resulting cone of depression can be thought of as a lowering of the 

water table around the well, with the greatest decline in the 

water table at the location of the well (fig. 3A). If the aquifer 

is confined, the change in storage is the combined effect of 

reduction of the sizes of pore spaces and expansion of water 

that results from decreasing the fluid pressure in the aquifer. 

Water table

Water table

Water table

Direction of groundwater 

movement

Direction of groundwater 

movement

Direction of groundwater 

movement

Aquifer

Aquifer

Aquifer

Stream

Stream

Stream

C

B

A

Cone of

depression

Cone of

depression

Cone of

depression

Figure 3. Stages of effects of a pumping well on a nearby stream. (A ) A cone of depression has developed around the well, but 

groundwater still discharges to the stream. (B ) After some time, the cone of depression has deepened and spread to the nearby 

stream, resulting in elimination of groundwater discharge to the stream and further loss of water in the stream to the aquifer. (C ) When 

the stream cannot supply the quantity of water pumped, the stream may eventually lose all of its water to the aquifer and become 

ephemeral and nearby riparian vegetation also may be affected. Modified from Webb and others (2007).

This resulting cone of depression is an area of decreased fluid 

pressure around the well, with the greatest decrease in fluid 

pressure at the well. 

As pumping continues, the cone of depression expands 

to increasing distances from the well. When the cone of 

depression expands into areas with streams (fig. 3B), wetlands, rivers, and lakes that are in hydraulic connection with 

the aquifer, and areas of vegetation that use groundwater, 

the natural gradients that drive the exchange of groundwater between the aquifer and these features are altered. For 

surface-water features that under natural conditions lose water 

to the aquifer, the cone of depression from a well can increase 

the gradients from those features; and for surface-water features that gain water from the aquifer, the cone of depression 
Introduction 5

can decrease the gradients to those features. Whether the case 

is a reduction of groundwater flow to a surface-water feature 

or an increase in surface-water flow to the aquifer system, 

there is a net loss in water in the feature. In areas where plants 

use groundwater (fig. 3), the cone of depression can lower the 

water table and reduce the uptake of water by these plants. The 

change in availability of water to these plants and the surfacewater features caused by a pumping well is generally referred 

to as “depletion” or “capture,” or as “streamflow depletion” 

where dominant changes are decreased outflow to or increased 

inflow from streams and rivers. In this report, surface-water 

bodies including streams, springs, and wetlands and areas of 

vegetation that use groundwater are referred to collectively as 

“connected features.”

The decreased outflow to and increased inflow from connected features is a source of water to the pumped well that 

tends to stabilize the cone of depression over time. If a well 

pumps for long enough at a constant rate and sufficient water 

can be captured from the connected features to supply the rate 

of water being pumped, the storage change in the aquifer from 

the pumping will diminish to zero and all of the water pumped 

by a well can be accounted for as changes in flow to or from 

connected features.

 The time over which all or nearly all of the groundwater 

pumping is supplied by connected features is dependent on 

(1) the proximity of the pumping wells to connected features 

that can supply water and (2) the hydraulic properties of the 

aquifer. In general, wells that are close to connected features 

will receive water from these features much faster than more 

distant wells. In cases where wells are many tens of miles 

from connected surface water, the time at which depletion is 

the dominant source of water to the well can be decades or 

even centuries after pumping begins. A graphical example of 

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

0 10 20 30 40 50 60 70 80 90 100

Decrease in storage

Decrease in surface-water flow plus 

decrease in evapotranspiration

TIME SINCE START OF GROUNDWATER PUMPING, IN YEARS

FRACTION OF WELL PUMPING RATE

this process is shown in figure 4. At time zero, when pumping begins, all of the water pumped is coming from storage 

in the aquifer. Over a period of 100 years in this example, the 

fraction of pumped water that comes from storage diminishes 

to 1 percent. The fraction of the pumping that is depletion 

(decrease in surface-water flow plus decrease in evapotranspiration) is the complement of the change in storage for any 

given time. For this example, depletion increases from zero 

initially to 99 percent at 100 years.

The aquifer properties that affect this timing are aquifer 

thickness, hydraulic conductivity, and storage coefficient. If 

an aquifer is dominated by horizontal flow, the relevant property that controls the rate at which pumped water is supplied 

from connected features is the ratio of the product of hydraulic 

conductivity and aquifer thickness to the storage coefficient—

the hydraulic diffusivity of the aquifer. In general, higher 

hydraulic conductivity and aquifer thicknesses and lower 

storage coefficients tend to accelerate the process of depletion. 

The timing of effects of pumping on surface water also may 

be affected by complexities including heterogeneous aquifer 

properties, faults or other structural features that act as flow 

barriers or conduits, and any intervening hydrogeologic units 

of low vertical hydraulic conductivity that inhibit propagation 

of pressure changes from pumping locations at depth to connected features. These hydrogeologic complexities may slow 

down or speed up the process of depletion and may determine 

the extent to which depletion will affect certain connected 

features. In a groundwater system with such complexities, 

however, it is unreasonable to expect that groundwater pumping would never deplete connected features.

Note that effects of pumping shown in figure 4 are given 

as a fraction of the pumping rate. In systems that respond 

linearly to pumping, the fractional response at any given time 

Figure 4. Sources of water to a pumped well through time.
6 Simulated Effects of Groundwater Pumping and Artificial Recharge in the Verde Valley Sub-Basin, Central Arizona

is independent of the pumping rate. The decrease in surfacewater flow plus reduction in evapotranspiration (uptake of 

groundwater by vegetation) for this example is a fraction of 

0.72 or 72 percent of the pumping rate at a time of 20 years 

since the start of groundwater pumping. If a system responds 

linearly to pumping, depletion at 20 years would be 7.2 gal/

min for a well pumping 10 gal/min and would be 72 gal/min 

for a well pumping 100 gal/min. Furthermore, in such a system the effects of pumping can be reversed to determine the 

effects of injecting water. If a well was injecting 100 gal/min 

for 20 years, for example, the increase in surface-water flow 

plus the increase in evapotranspiration would be 72 gal/min. 

Nonlinearities occur when pumping (or injection of water) 

change the properties of the system. The most common pumping- or recharge-induced change is a large decrease or increase 

in saturated thickness in unconfined aquifers. Another nonlinearity occurs when pumping causes surface water or plants 

that use groundwater to become disconnected from the aquifer 

(fig 3C). Most aquifers respond linearly over some range of 

pumping, but become nonlinear with some level of increased 

pumping, especially if that pumping is concentrated in a few 

areas (Leake and others, 2010). 

The Northern Arizona Regional GroundwaterFlow Model

The timing of effects of groundwater pumping in the 

Verde Valley sub-basin on connected features in and around 

the sub-basin can be evaluated using a groundwater-flow 

model. The Northern Arizona Regional Groundwater-Flow 

Model (NARGFM) by Pool and others (in press) is the most 

comprehensive and up-to-date groundwater-flow model that 

includes the sub-basin. With an active area of more than 

51,000 mi2, NARGFM simulates groundwater flow in an area 

much larger than the sub-basin, including parts of the Colorado Plateau, the Verde River Basin, and adjacent areas (fig. 

1). Despite the large extent of the model, it is well suited to 

simulate conditions in the Verde Valley sub-basin because it 

does not include artificial head-dependent or specified-flow 

boundary conditions and because the uniform horizontal grid 

spacing of 3,281 ft (1,000 m) is appropriate for evaluating system responses to pumping in an area the size of the sub-basin.

The model was developed using the U.S. Geological Survey (USGS) groundwater model MODFLOW-2000 (Harbaugh 

and others, 2000) and MODFLOW-2005 (Harbaugh, 2005). 

In the horizontal dimension, the model grid includes 600 rows 

and 400 columns of finite-difference cells. The geographic 

origin of the grid, at the outer corner of the cell at row 600 and 

column 1 is at Universal Transverse Mercator zone 12 easting 660,000 m and northing 3,580,000 m. The model grid is 

rotated 60° in the counter-clockwise direction about this point 

so that model columns are incremented from the cell at row 1, 

column 1, to column 400 in a northeasterly direction. Similarly, rows are incremented from the cell at row 1, column 1 to 

row 600 in a southeasterly direction. The 60° rotation is made 

so that rows and columns are aligned with the assumed principal directions of the hydraulic conductivity tensor. The model 

was constructed using units of meters for length and days for 

time. For this report, descriptions of values in the model are 

given in units of feet for length and days for time.

NARGFM uses three model layers to simulate different 

hydrogeologic units. The hydrogeologic units simulated within 

each of the three model layers vary across the model domain. 

For details, see Pool and others (in press). Within the Verde 

Valley sub-basin, the major hydrogeologic unit included in the 

upper model layer (layer 1) is the fluviolucustrine facies of the 

Verde Formation (fig. 5). The major units in the middle layer 

(layer 2) are the sand and gravel facies of the Verde Formation, volcanic rocks, and the upper, middle, and lower Supai 

Formations. The major units of the lower model layer (layer 3) 

are the Redwall Limestone and crystalline rock.

The active unconfined and confined model areas differ for each of the three model layers—layer 1 is active in a 

small area adjacent to the Verde River below Clarkdale, layer 

2 is active over much of the central and eastern parts of the 

sub-basin, and layer 3 is active over nearly all of the sub-basin 

(fig. 6). Elevations of tops and bottoms of model layers were 

estimated from land-surface topography and geology, and any 

available well-log and geophysical data. Information for layer 

elevations is limited in undeveloped areas and at depths not 

penetrated by wells.

Important properties of the aquifer included in the model 

are hydraulic conductivity, specific storage, and specific yield. 

Figure 5. Conceptualized relations among hydrogeologic units 

and Northern Arizona Regional Groundwater-Flow Model layers 

in the Verde Valley sub-basin. Volcanic rocks and the lower Supai 

Group are confining units. Modified from Pool and others (in press).

QUATERNARY ALLUVIUM

FLUVIOLACUSTRINE FACIES OF

THE VERDE FORMATION

SAND AND GRAVEL

FACIES

OF THE

VERDE

FORMATION

VOLCANIC

ROCKS

UPPER AND MIDDLE

SUPAI FORMATIONS

LOWER SUPAI

FORMATION

REDWALL LIMESTONE AND OTHER

CARBONATE ROCKS

CRYSTALLINE ROCK

NOT SIMULATED

LAYER 1

LAYER 2

LAYER 3

HYDROGEOLOGIC UNITS MODEL LAYERS

BASALT
Introduction 7

Locations of active unconfined model cells

Locations of active confined model cells

Boundary of Verde Valley sub-basin

EXPLANATION

UTM

Zone 12

NAD27

EASTING, IN METERS

EASTING, IN METERS

NORTHING, IN METERS NORTHING, IN METERS

Layer 3

Layer 2

3800000

3840000

3880000

3920000

340000 380000 420000 460000

340000 380000 420000 460000

3800000

3840000

3880000

3920000

Layer 1

These are represented in the Layer-Property Flow Package 

of MODFLOW. Hydraulic conductivity is a tensor quantity 

that can be different in principal directions. In the area of the 

Verde Valley sub-basin in NARGFM, hydraulic conductivity 

is the same in both horizontal directions in layer 1, but it is 

higher along columns than along rows in layers 2 and 3 (table 

1). Having different hydraulic conductivity values in principal 

horizontal directions is referred to as “horizontal anisotropy.” 

In this case, horizontal anisotropy represented in the model 

results in more resistance to groundwater movement in the 

northeast-southwest direction than in the northwest-southeast 

direction. Mean vertical hydraulic conductivity in all layers 

is lower than mean hydraulic conductivity along rows and 

columns, and is lower in layer 1 than in layers 2 and 3 (table 1). 

In MODFLOW, specific storage is applied where the aquifer is 

confined, as indicated by head or water level above the top of 

the aquifer. The product of specific storage and saturated thickness is the aquifer storage coefficient. The USGS has defined 

the storage coefficient of an aquifer as “the volume of water it 

releases from or takes into storage per unit surface area of the 

aquifer per unit change in the component of head normal to 

that surface” (http://water.usgs.gov/admin/memo/GW/gw55.28.

html, accessed February 7, 2010). The mechanisms for these 

storage changes are compression and decompression of the 

aquifer skeleton and water. Similarly, specific yield is applied 

where the aquifer is unconfined as indicated by head or 

water level below the top of the aquifer. Specific yield is the 

quantity of water released per unit unit volume of the aquifer 

by gravity drainage from lowering the water table. In parts of 

NARGFM within the Verde Valley sub-basin layer 3 that are 

overlain by saturated parts of layer 2, the aquifer is confined 

and storage changes are calculated in the model using the product of specific storage and saturated thickness. In the remaining 

active part of layer 3, the aquifer is unconfined and specific 

yield is used. Similarly for layer 2, the part overlain by the 

saturated area of layer 1 is confined and the rest of the active 

Figure 6. Active unconfined and confined model cells in the Northern Arizona Regional Groundwater-Flow Model 

within the Verde Valley sub-basin for model layers 1, 2, and 3.
8 Simulated Effects of Groundwater Pumping and Artificial Recharge in the Verde Valley Sub-Basin, Central Arizona

part is unconfined. The entire active area of layer 1 is unconfined. Mean specific yield in the sub-basin is similar in layers 1 

and 2 but is an order of magnitude lower in layer 3 (table 1).

Saturated thickness also is an important factor in calculating groundwater movement in the model. Where heads are 

below the top of the model layer, saturated thickness is the difference between the model-calculated head or water level and 

the specified bottom of the aquifer. Mean saturated thickness 

is much larger in layers 2 and 3 than in layer 1 (table 1).

Discharge of water in NARGFM occurs through 

streams, simulated with the Stream Package (Prudic, 1989), 

and the Evapotranspiration and Drain Packages (Harbaugh, 

2005). Within the Verde Valley sub-basin, two active features (springs) were simulated with the Drain Package 

(fig. 1) and all streams and other springs are simulated with 

the Stream Package. Although most connected features 

simulated with the Drain Package are outside of the subbasin, it is possible that some depletion could occur in these 

distant features from pumping in the sub-basin. Depletion 

from pumping in the sub-basin, however, is more likely 

to occur from much closer connected features within the 

sub-basin represented with the Stream and Evapotranspiration Packages. The Stream Package keeps track of available 

flow and calculates stream stage in each section or reach of 

stream that crosses a model cell. If leakage of water from 

the stream exceeds available streamflow, no further leakage is allowed from downstream reaches until the point 

where there is additional water available to the downstream 

reaches from tributary inflow or flow of groundwater into 

the stream. Input for the Stream Package includes specification of quantities for each of the stream reaches such as the 

elevations of the top and bottom of the streambed, streambed conductance, stream width, and Manning’s roughness 

coefficient. For details on representation of features with 

the Stream, Drain, and Evapotranspiration Packages in the 

NARGFM, see Pool and others (in press).

Scope

The scope of this study includes modification of 

NARGFM for the purpose of studying possible effects of 

groundwater pumping on connected surface water and riparian 

plants that use groundwater. The main goal is to provide data 

sets and maps to support the development of sustainable water 

management in the Verde Valley that addresses both human 

needs and desired streamflow conditions. The data sets are 

comprehensive, including results of pumping groundwater in 

the Verde Valley sub-basin in all three model layers for steadystate pumping conditions, and for transient conditions during 

100 years of pumping. Maps presented in this report show the 

effects of pumping 10 and 50 years in active parts of models 

layers 1 and 2. Maps showing the effects of pumping any of 

the model layers for any time during 100 years of pumping or 

for steady-state pumping conditions can be made using data 

sets from this study.

Method

Methods and considerations for using groundwater 

models to show effects of groundwater pumping or recharge 

on surface water and riparian vegetation are given by Leake 

and Reeves (2008) and Leake and others (2008, 2010). Similar 

methods were used for this study. The procedure for assessing 

effects of recharge on enhancement of surface flow and evapotranspiration is the same as the procedure for assessing depletion of surface-water flow and evapotranspiration, except for 

the directions of the forcing action (recharge to or discharge 

from groundwater) and response (increased or decreased 

surface-water flow and evapotranspiration). For simplicity, this 

discussion will reference the more common case of discharge 

(pumping) from groundwater. 

The procedure used here calculates change from an 

assumed steady-state (predevelopment) condition. Inherent in 

this approach is the assumption that natural variability, such 

as changes in recharge, and development of groundwater 

resources do not greatly change the aquifer transmissivity and 

configuration of connected features (such as connected surface 

water and evapotranspiration areas). If this assumption is 

reasonable, then results from the method may be used to look 

at changes from transient conditions accompanying natural 

system variations and groundwater-resources development.

To use a groundwater model to assess effects of groundwater pumping on connected surface water and evapotranspiration at a location of interest, the model must first be run 

without the added pumping. For that run, all components of 

flow to and from surface-water features and evapotranspiration 

areas must be recorded. If the analysis is transient, water-budget components must be recorded for points in time at which 

effects of the added well are to be calculated. 

A brief outline of the step-wise procedure used to 

construct the maps showing the effect of well location on the 

reduction of surface-water flow and (or) evapotranspiration 

from groundwater pumping as applied in this study follows:

1. Run the steady-state model without added 

pumping, and save values of groundwater flow 

to or from groups of stream reaches and overall 

groundwater discharge by riparian vegetation 

within the Verde Valley sub-basin.

2. For an active cell location in the sub-basin, add 

a well to the MODFLOW Well Package. The 

pumping rate is set at 18.3 gal/min (100 m3/d).

Run a transient model with the added well for a 

100-year period.

3. At 1-year intervals from 1 to 100 years, compute and save changes in flow components from 

results calculated without the well in the steadystate model (step 1). Changes in flow to and 

from groups of stream reaches and in discharge 

by riparian vegetation are divided by the well 
Method 9

pumping rate to get the effect of pumping as 

a fraction of the pumping rate. In addition to 

changes in flow rates at these features, the location of the well is saved.

4. When the model has been run for all desired 

locations of added wells, map fractions of total 

reduction in surface water plus reduction in 

evapotranspiration for the area over which calculations were made and saved.

In this procedure, steps 2–4 are repeated until results 

are obtained for each desired well location. In this application, wells are placed sequentially at every active model cell 

in the area to achieve the maximum resolution possible in 

mapping the effects of pumping. The horizontal grid spacing 

of the model is 3,281 ft (1,000 m); therefore, sequences of 

steps 2–4 resulted in computed values on a horizontal grid 

with this spacing along the orientation of model rows and 

columns. Mapping for this area in step 5 includes plotting 

pumping-induced changes from wells in active model cells in 

model layers 1 and 2 in the Verde Valley sub-basin at pumping times of 10 and 50 years. From results computed and 

saved, different maps may be generated in the future, including pumping-induced changes from wells in any of the model 

layers, at any pumping time from 1 to 100 years.

Property Units Statistic Value in layer 1 Value in layer 2 Value in layer 3

Hydraulic conductivity along rows feet per day Mean 5.39×100 2.03×100 1.16×100

Minimum 3.28×10-1 1.31×10-2 3.28×10-3

Maximum 1.31×101 3.94×101 1.31×101

Hydraulic conductivity along 

columns feet per day Mean 5.39×100 4.54×100 1.60×100

Minimum 3.28×10-1 5.91×10-3 3.28×10-3

Maximum 1.31×101 1.18×102 6.56×101

Vertical hydraulic conductivity feet per day Mean 2.77×10-3 2.91×10-2 1.16×10-1

Minimum 1.31×10-4 1.31×10-5 3.28×10-4

Maximum 6.56×10-3 1.31×100 1.31×100

Specific storage1 per foot Mean 3.05×10-4 3.04×10-5 3.05×10-7

Minimum 3.05×10-4 3.05×10-6 3.05×10-7

Maximum 3.05×10-4 3.05×10-5 3.05×10-7

Specific yield1 — Mean 1.00×10-1 1.13×10-1 1.00×10-2

Minimum 1.00×10-1 5.00×10-2 1.00×10-2

Maximum 1.00×10-1 2.50×10-1 1.00×10-2

Saturated thickness Feet Mean 2.06×102 1.32×103 1.87×103

Minimum 9.15×100 1.20×101 4.49×102

Maximum 7.69×102 4.27×103 4.86×103 

Table 1. Statistics of aquifer properties and saturated thickness simulated in the part of the Northern Arizona Regional 

Groundwater-Flow Model that underlies the Verde Valley sub-basin. Values are for active parts of each model layer.

The model documented by Pool and others (in press) 

included one steady-state stress period and nine transient 

stress periods. The model was not directly usable for making 

the steady-state model run in step 1 and the many transient 

model runs in step 3. Modifications made for the purposes of 

this study follow:

A. A steady-state model was constructed by using 

the first steady-state stress period in NARGFM 

followed by a dummy transient stress period of 1 

day in length. The dummy stress period allowed 

direct use of the Layer-Property Flow Package 

input data sets from NARGFM. To maximize the 

mass balance of the steady-state components, 

closure criteria for head change and residual 

mass balance (input parameters HCLOSE and 

RCLOSE) of the Preconditioned Conjugate-Gradient solver were lowered, and the model was 

run with a double-precision version of MODFLOW-2005 (Harbaugh, 2005). NARGFM originally was set up with the capability of allowing 

rewetting of model cells during the iterative 

solution. The rewetting option, however, limited 

the accuracy of the mass balance because a number of cells would change from wet to dry and 

1

In computing statistics for this storage property, values at locations where the property is not used by the model were excluded.
10 Simulated Effects of Groundwater Pumping and Artificial Recharge in the Verde Valley Sub-Basin, Central Arizona

dry to wet on alternate iterations. For this study, 

the computed head and IBOUND arrays were 

saved from an initial solution with the rewetting option turned on. For that solution, the final 

status of the alternating cells was dry, or inactive. 

A subsequent steady-state run was made using 

these head and IBOUND arrays with the rewetting option turned off. In doing so, computed 

inflow and outflow for the revised steady-state 

model of about 1.27 ×107 ft3/d were in agreement to within about 3 ft3

/d (< 3.0 ×10-5 percent). 

Computed head values from this steady-state 

model were saved for use as the starting head in 

the transient model. The implication of removing 

the rewetting option for the steady-state model is 

that certain cells with limited saturated thickness 

are removed from the active flow domain. These 

cells occur around the perimeter of the active 

model domain and are not important to include 

in the model analysis completed for this project. 

The transient analysis of depletion is appropriate 

to use without the rewetting option because there 

are no stresses being simulated that would cause 

head to rise and rewet previously dry cells.

B. NARGFM input specified locations of 1,957 

drain cells. Some of these, however, were inactive in the steady-state solution because computed head was below the specified drain elevation. For this study, inactive drain cells were 

eliminated, resulting in 1,245 active drain cells, 2 

of which are in the Verde Valley sub-basin. These 

remaining active cells were classified in 9 groups 

for which drain observations were set up to 

record pumping-induced change in flow for each 

time step using the Drain Observation Package. 

C. NARGFM input included 377 segments (groups 

of reaches) for the Stream Package, but only 

312 segments simulated flow to or from the 

aquifer in the steady-state model. For those 

312 segments, the MODFLOW Stream Observation Package (http://water.usgs.gov/nrp/

gwsoftware/modflow2000/MFDOC/index.

html, accessed Feb. 4, 2010) was used to record 

flows exchanged with the groundwater system. 

Computed flow to and from each of these stream 

segments was saved for the final precise steadystate model run. Components from the overall 

MODFLOW mass balance also were saved.

D. The transient model used in step 3 was constructed from the steady-state model by (a) using 

the heads from the precise steady-state run as 

starting heads, (b) changing the single steadystate stress period to a single transient stress 

period with 100 1-year time steps, (c) adding a 

single well as outlined in step 2, (d) computing 

flow values from the steady-state run for the 312 

stream segments that were used as observations 

of interchange between the segments and the 

groundwater system at 1-year intervals, and (e) 

computing flow from the steady-state run for each 

of the 9 groups of drain cells that were used as 

observations of discharge from the groundwater 

system to drains at 1-year intervals. For (d) and 

(e), Stream and Drain Observation Packages were 

used. The transient model also was run with a 

double-precision version of MODFLOW-2005. 

In using the steady-state values of flow to or from 

stream reaches as observations, calculated differences between computed and observed flows are 

the direct effect of the added pumping well.

In addition to the procedure described above for calculating the effects of pumping through time since pumping began, 

the effects of the new steady-state condition (when storage 

change is zero and all of the water pumped is from reduced 

surface water and reduced evapotranspiration) were calculated, 

regardless of the time required to reach this condition. This 

steady-state approach is useful in evaluating effects on multiple 

connected features. For example, if a system included stream 

A and stream B, pumping at some locations (presumably near 

stream A) would tend to ultimately affect only stream A by the 

amount of the well pumping rate. For other locations, pumping 

would affect both streams, with the sum of the effect on each 

stream equal to the total pumping rate (where depletion from 

features simulated with the Evapotranspiration and Drain Packages are zero). The procedure for carrying out the steady-state 

analysis is similar to that for the transient analysis outlined 

above, except that the transient model in step 3 is replaced with 

a steady-state model. Instead of effects at 1-year intervals up to 

100 years, this model calculates a single ultimate effect when 

steady-state conditions under the imposed pumping stress are 

reached, without any indications of how many years, decades, 

or centuries it might take to reach that condition. 

Results

In model layers 1, 2, and 3, pumping from a single 

hypothetical well was simulated sequentially in each of 10,789 

active model cells. This was achieved by running the model 

10,789 times. The repeated model runs produced a variety of 

results, only a small subset of which is presented in this report. 

All results are depletion values, expressed as a fraction of the 

pumping rate for each of 10,789 active model cells in layers 1–3 in the Verde Valley sub-basin. Most transient results 

include depletion values at 1-year intervals from 1–100 years 

for each cell and steady-state results include a single ultimate 

depletion value for each cell.The complete results of model 

runs include the following:
Results 11

1. Total reduction in streamflow (simulated with the 

Stream and Drain Packages) plus total reduction 

in evapotranspiration, simulated in the transient 

model.

2. (a) Rate of change in groundwater storage, (b) 

reduction in outflow to features simulated with 

the Drain Package, (c) reduction in outflow to 

features simulated with the Evapotranspiration 

Package, and (d) reduction of net flow between 

the aquifer and features simulated with the 

Stream Package at 100 years, simulated in the 

transient model.

3. Total reduction in streamflow for each of 312 

stream segments, simulated in the transient 

model.

4. Total reduction in streamflow for each of 312 

stream segments, simulated in the steady-state 

model.

5. Total reduction in groundwater outflow for each 

of 9 groups of drain cells, simulated in the transient model.

6. Total reduction in groundwater outflow for 

each of 9 groups of drain cells, simulated in the 

steady-state model.

These results are in ASCII files archived with the models used 

for this study and are available upon request.

By assuming linearity in response to pumping, results 

in item 1 can be used to compute effects of more complex 

scenarios than a single well pumping at a constant rate for a 

period of time. To compute additive effects, fractions in item 

1 must be first multiplied by pumping rates to calculate actual 

depletion. For given times, actual depletions from different 

pumping locations may be added to calculate the total depletion from all of the pumping locations. It is also possible to 

compute the effects of pumping at a variable rate at a single 

well location, including the effects of a well that is pumped 

intermittently. For example, data from the curve of reduced 

surface-water flow and evapotranspiration shown in figure 4 

may be used to compute absolute effects for a well that pumps 

30 acre-ft per year for years 1–5 and then for years 6–10 the 

pumping rate continues (case 1), stops (case 2), and doubles 

(case 3). The fractions that make up the curve for times 0–10 

years in figure 4 are designated as line 1 in table 2. For case 1, 

fractions in line 1 are multiplied by 30 acre-ft/year to compute 

values in line 2, the absolute reduction in surface-water flow 

and evapotranspiration, in acre-ft/yr, that would result from 

this well pumping at a constant rate for 1–10 years. To obtain 

the effects of a change in pumping rate at 5 years, the pumping rates for times 1–5 years must be put into the positions of 

times 6–10 years as shown in line 3 in table 2. These values 

are equivalent to responses from a hypothetical well at this 

location that counteracts the pumping of the existing well for 

case 2 or adds to the pumping rate of the existing well for 

case 3. For case 2, values in line 3 are subtracted from values 

in line 2 to calculate the response for the well ceasing to be 

pumped, and for case 3, values in line 3 are added to values in 

line 2 to calculate the response for a doubling of the pumping 

rate (table 2, fig. 7).

In addition to these results, some information was generated to indicate possible nonlinearity for each of the 10,789 

pumping locations. This included a code to denote where 

pumping caused the cell being pumped to go dry and the number of stream reaches (sections of streams that cross a model 

cell), if any, for which the connection to the aquifer became 

nonlinear because of the added well pumping 18.3 gal/min. 

Dry pumping cells occurred mostly around the margins of the 

active part of model layer 1 and in some locations in model 

layer 2 where initial saturated thickness was low. Nonlinear 

0

5

0 1 2 3 4 5 6 7 8 9 10

10

15

20

25

DECRESE, IN ACRE-FEET PER YEAR

TIME, IN YEARS

Case 3

Case 1

Case 2

Figure 7. Decrease in surface-water flow and evapotranspiration for the example shown 

in figure 4 for a well that pumps 30 acre-feet per year for years 1–5; then for years 6–10 the 

pumping rate continues (case 1), stops (case 2), and doubles (case 3).
12 Simulated Effects of Groundwater Pumping and Artificial Recharge in the Verde Valley Sub-Basin, Central Arizona

connection between stream reaches and the aquifer can occur 

for two conditions. The first condition occurs when leakage 

of water to the aquifer is limited by available streamflow and 

the second condition occurs when aquifer head falls below 

the specified streambed bottom elevation and leakage to the 

aquifer becomes steady. At a simulation time of 100 years, 

one or both of these conditions occurred at 1 stream reach for 

22 pumping locations and at 2 stream reaches for 33 pumping locations. Further analysis of nonlinearity would require 

recomputing effects at all locations with different pumping 

rates. Such analyses are beyond the scope of this study.

Results presented here for pumping locations in model 

layers 1 and 2 show the fraction of the 18.3 gal/min pumping 

rate that is computed to be a reduction in surface-water flow 

and evapotranspiration after 10 and 50 years since the onset of 

pumping. Results for model layer 1 at 10 years (fig. 8) indicate 

that the greatest effects occur along the Verde River between 

Clarkdale and Cottonwood and from several miles upstream of 

Camp Verde to the sub-basin boundary. Values also are higher 

along Oak Creek upstream and downstream from Cornville. 

Values are lower along the Verde River along an 11-mile reach 

of the river below Cottonwood. Results for model layer 1 at 

50 years (fig. 9) indicate similar patterns of higher and lower 

values but with higher fractions of decrease in surface-water 

flow and evapotranspiration because of the longer pumping 

time. Nonetheless, model results show fractions of reduction 

in flow and evapotranspiration of less than 0.1 (10 percent of 

Line 

number

Item Result

Time, in years 0 1 2 3 4 5 6 7 8 9 10

1 Depletion fractions 

from first 10 years 

of graph shown in 

figure 4

0.00 0.06 0.12 0.18 0.23 0.28 0.32 0.36 0.40 0.44 0.47

2 Values in line 1 times 

pumping rate of 30 

acre-feet per year

0.00 1.94 3.68 5.31 6.84 8.26 9.64 10.88 12.04 13.18 14.19

3 Values from years 

0–5 in line 2 put 

into columns for 

years 5–10

0.00 0.00 0.00 0.00 0.00 0.00 1.94 3.68 5.31 6.84 8.26

4 Values in line 2 minus 

values in line 3

0.00 1.94 3.68 5.31 6.84 8.26 7.69 7.20 6.72 6.34 5.93

5 Values in line 2 minus 

values in line 3

0.00 1.94 3.68 5.31 6.84 8.26 11.58 14.55 17.35 20.02 22.45

Table 2. Example calculations of absolute depletion from depletion fractions for a well that pumps 30 acre-feet per 

year for years 1–5, and then continues pumping at that rate for years 6–10 (case 1), shuts off for years 6–10 (case 2), 

and doubles in pumping rate for years 6–10 (case 3).

[Results for case 1, years 0–10, are shown in line 2; results for case 2, years 6–10 , are shown in the last five columns of line 4; results for 

case 3, years 6–10 , are shown in the last five columns of line 5.]

the pumping rate) at 50 years for some areas along the north 

and west edges of the active part of model layer 1 (fig. 9).

Fractions of reduced surface-water flow and evapotranspiration for pumping locations in model layer 2 at 10 years 

(fig. 10) indicate highest responses along the Verde River 

upstream of Cottonwood and downstream of Camp Verde, 

along reaches of Oak Creek, and along Wet Beaver Creek 

upstream of Rimrock. Only the most upper and lower parts 

of West Clear Creek are connected to model layer 2. These 

results, therefore, do not indicate effects of groundwater 

pumping on most of West Clear Creek. Results for model layer 

2 at 50 years (fig. 11) indicate general increases in effects 

along the surface-water features and at greater distances from 

these features, compared to results at 10 years. A complex 

pattern of lower response exists in part of this area overlain by 

model layer 1, east and west of Cornville. Patterns of response 

generally are correlated with properties such as hydraulic conductivity. For example, some of the highest areas of depletion 

from pumping in model layer 2 are in areas were hydraulic 

conductivity along rows is greater than 5 ft/d (figs. 10–11). 

Other high areas are adjacent to streams.

In addition to evaluating possible effects of groundwater pumping, figures 8–11 also can be used to understand the 

effects of managed recharge on surface-water flow and evapotranspiration. For recharge at a constant rate, areas with higher 

fractions (orange and red) for a given time indicate a greater 

enhancement of surface-water flow and evapotranspiration 
Results    13

400000 420000

3820000 3840000

UTM

ZONE 12

NAD27

0 2.5 5 10 Miles

Verde

River

Oak

Creek

West Clear Creek

Wet

Beaver

Creek

Total reduced flow as a fraction of pumping rate at 10 years

Verde River and major tributary streams

Major Roads

Verde Valley sub-basin

EXPLANATION

0-0.1

0.1-0.20.2-0.30.3-0.40.4-0.50.5-0.60.6-0.7

0.7-0.80.8-0.90.9-1.0

CLARKDALE

JEROME COTTONWOOD

CORNVILLE

RIMROCK

CAMP

VERDE

Figure 8. Computed reduction in flow of streams and springs and reduction in riparian evapotranspiration as a fraction of

pumping rate that would result from pumping groundwater from model layer 1 at a constant rate for 10 years. The color at

any location represents the fraction of the pumping rate by a well at that location that can be accounted for as changes in

outflow from or inflow to the aquifer for model boundaries representing steams, riparian vegetation and springs.
14    Simulated Effects of Groundwater Pumping and Artificial Recharge in the Verde Valley Sub-Basin, Central Arizona

Figure 9.  Computed reduction in flow of streams and springs and reduction in riparian evapotranspiration as a fraction of

pumping rate that would result from pumping groundwater from model layer 1 at a constant rate for 50 years. The color at

any location represents the fraction of the pumping rate by a well at that location that can be accounted for as changes in

outflow from or inflow to the aquifer for model boundaries representing streams, riparian vegetation and springs.

Verde

River

400000 420000

3820000 3840000

UTM

ZONE 12

NAD27

0 2.5 5 10 Miles

Oak

Creek

West Clear Creek

Wet

Beaver

Creek

Total reduced flow as a fraction of pumping rate at 50 years

Verde River and major tributary streams

Major Roads

Verde Valley sub-basin

EXPLANATION

0-0.1

0.1-0.2

0.2-0.3

0.3-0.40.4-0.50.5-0.60.6-0.7

0.7-0.80.8-0.9

0.9-1.0

CLARKDALE

JEROME

COTTONWOOD

CORNVILLE

RIMROCK

CAMP

VERDE
Summary 15

than areas with lower fractions (blues). If recharge is by gravity drainage to the water table, results from confined parts of 

model layers 2 and 3 should not be used.

The relative contributions to total depletion from features 

simulated with the Drain, Evapotranspiration, and Stream 

Packages after 100 years can be seen in results 2b–d listed 

above. Average values of depletion by type of feature and by 

model layer (table 3) indicate that most of the depletion occurs 

in features represented with the Stream Package for pumping 

in all model layers.

Applicability and Limitations

As indicated by Leake and others (2008), the procedure 

used for the analysis described in this report was designed to 

depict responses to pumping or recharge over large areas of 

an aquifer. The results are meant to help resource managers 

and the general public understand patterns of responses that 

are based on features represented in the part of the NARGFM 

underlying the Verde Valley sub-basin. The distributions 

shown on figures 8–11 do not necessarily mean that groundwater pumping and (or) recharge to the aquifer at all locations 

is appropriate or even technically feasible. In the generation 

of maps presented in this report, an automated procedure 

computed effects at all active model cells regardless of properties such as hydraulic conductivity and storage coefficient. 

For site-specific projects involving groundwater pumping or 

recharge to the aquifer, improved results may require detailed 

studies of the local hydrologic conditions and a refinement of 

the modified NARGFM in the area of interest.

As discussed previously, results are independent of the 

pumping or recharge rate if the system responds linearly to 

these changes. Although comprehensive testing of linearity was 

not done for this study, other withdrawal rates should produce 

patterns of fractions of decrease in surface-water flow and 

evapotranspiration that are similar to those shown in figures 

8–11, particularly for pumping or recharge rates that are smaller 

than the rate of 18.3 gal/min used in this analysis. For example, 

with higher pumping rates, there is an increased likelihood of 

nonlinearity from reducing saturated thicknesses, from drying stream reaches, and from lowering water levels below the 

evapotranspiration extinction depth. Changes are most likely 

to occur in areas where saturated thickness is lowest and near 

streams where base flow is sufficiently low to be entirely 

depleted by the pumping. If nonlinearity exists, it will affect 

the timing of effects on surface-water flow and evapotranspiration, but nonlinearities will not change the tendency towards 

the ultimate effect of the entire pumping rate being accounted 

for as reduction in surface-water flow and evapotranspiration. 

If nonlinearities exist, responses to pumping rates higher than 

18.3 gal/min will be overestimated by results in this report, and 

responses to higher recharge rates will be underestimated.

The analyses use a constant pumping rate for a period 

of 100 years. In reality, wells are not managed to pump at a 

Layer Drain 

Package

Evapotranspiration 

Package

Stream 

Package

Number 

of 

solutions

1 2.41×10-4 8.34×10-2 6.93×10-1 384

2 6.17×10-4 1.12×10-2 3.44×10-1 3,959

3 2.39×10-2 1.11×10-2 5.21×10-1 6,037 

Table 3. Average depletion after 100 years of pumping expressed 

as a fraction of the pumping rate for features simulated with the 

Drain, Evapotranspiration, and Stream Packages, for pumping 

locations in layers 1, 2, and 3.

constant rate, but rather are pumped in response to demands 

for water. The responses in this study approximate the response 

for a well location in which a time-average of the variable rate 

of pumping is used as a constant rate of pumping. In areas 

where substantial depletion of surface water and evapotranspiration occurs more quickly (for example in the red areas of figures 8 and 10), approximation of effects with a time-averaged 

rate might not represent the possible maximum depletion of 

surface water and evapotranspiration. If more exact responses 

are needed for cyclical or intermittent pumping, the procedure 

for more complex scenarios described above (table 2, fig. 7) 

may be applied. The use of 1-year time steps for this analysis, 

however, limits operational periods to multiples of 1-year. 

Analysis of shorter periods would require interpolation of 

model results from this study, or a separate study.

Summary

In the Verde Valley sub-basin, groundwater is the primary 

source of domestic and municipal supply. Groundwater use 

has increased in recent decades and a better understanding of 

the effects of groundwater pumping and managed recharge on 

surface water and groundwater-dependent riparian vegetation 

is needed. The Northern Arizona Regional Groundwater-Flow 

Model by Pool and others (in press) is the most comprehensive and up-to-date tool available to understand the effects 

of groundwater pumping in the sub-basin. For this study, 

the model was modified to calculate the combined effect 

of reduced surface-water flow and evapotranspiration from 

groundwater pumping for locations in the sub-basin where 

model cells were active. The procedure by Leake and others 

(2008) was used to run the model repeatedly, each time with 

a single pumping well in a different location for each active 

model cell. For each location, effects of the single well on 

surface-water flow and evapotranspiration were calculated 

and saved. Results for pumping times of 10 and 50 years for 

locations in active parts of model layers 1 and 2 were used 

to map total reduced surface-water flow and evapotranspiration as a fraction of pumping rate after 10 and 50 years (figs. 

8–11). These results also can be used to evaluate response of 

surface-water flow and evapotranspiration to recharge at a 
16 Simulated Effects of Groundwater Pumping and Artificial Recharge in the Verde Valley Sub-Basin, Central Arizona

400000 450000

3850000 3900000

Total reduced flow as a fraction of pumping rate at 10 years

UTM

ZONE 12

NAD 27 0 5 10 20 Miles

Verde

Oak

Sycamore

West

Clear

Creek

Wet Beaver Creek

Creek

Creek

River

Verde River and major tributary streams

Major Roads

Verde Valley sub-basin

EXPLANATION

Area where hydraulic conductivity along

rows in layer 2 is greater than 5 feet per day

JEROME

CAMP

VERDE

RIMROCK

CORNVILLE

CLARKDALE

SEDONA

FLAGSTAFF

0-0.1

0.1-0.2

0.2-0.30.3-0.40.4-0.50.5-0.60.6-0.70.7-0.80.8-0.9

0.9-1.0

Figure 10. Computed reduction in flow of streams and springs and reduction in riparian evapotranspiration as a fraction 

of pumping rate that would result from pumping groundwater from model layer 2 at a constant rate for 10 years. The color 

at any location represents the fraction of the pumping rate by a well at that location that can be accounted for as changes 

in outflow from or inflow to the aquifer for model boundaries representing streams, riparian vegetation and springs.
Applicability and Limitations 17

Figure 11. Computed reduction in flow of streams and springs and reduction in riparian evapotranspiration as a fraction 

of pumping rate that would result from pumping groundwater from model layer 2 at a constant rate for 50 years. The color 

at any location represents the fraction of the pumping rate by a well at that location that can be accounted for as changes 

in outflow from or inflow to the aquifer for model boundaries representing streams, riparian vegetation and springs.

400000 450000

3850000 3900000

Total reduced flow as a fraction of pumping rate at 50 years

UTM

ZONE 12

NAD 27 0 5 10 20 Miles

Verde

Oak

Sycamore

West

Clear

Creek

Wet Beaver Creek

Creek

Creek

River

Verde River and major tributary streams

Major Roads

Verde Valley sub-basin

EXPLANATION

Area where hydraulic conductivity along

rows in layer 2 is greater than 5 feet per day

JEROME

CAMP

VERDE

RIMROCK

CORNVILLE

CLARKDALE

SEDONA

FLAGSTAFF

0-0.1

0.1-0.20.2-0.30.3-0.40.4-0.50.5-0.60.6-0.70.7-0.80.8-0.9

0.9-1.0
18 Simulated Effects of Groundwater Pumping and Artificial Recharge in the Verde Valley Sub-Basin, Central Arizona

constant rate. The highest and most rapid simulated responses 

to pumping occurred near surface-water features simulated in 

the model, but results were not uniform along these features. 

For example, higher responses for pumping times of 10 and 50 

years for locations in both layers 1 and 2 were calculated adjacent to the Verde River upstream of Cottonwood and from the 

Camp Verde area to the boundary of the sub-basin. Responses 

along the intermediate part of the Verde River at these times 

were lower. The results shown in figures 8–11 indicate general 

patterns over large areas, as computed by the modified model. 

Investigations of pumping or recharge at specific sites may 

require more detailed study.

Acknowledgments

Jeanmarie Haney and Brenda Burman of The Nature 

Conservancy, Arizona, provided valuable insights on the 

types of model results that would support the evaluation of 

water management strategies in the Verde Valley sub-basin. 

Thoughtful reviews by Steve Peterson and Nora Herrera of 

the USGS, and Jeanmarie Haney greatly improved the report.

References Cited

Blasch, K.W., Hoffmann, J.P., Graser, L.F., Bryson, J.R., and Flint, 

A.L., 2006, Hydrogeology of the upper and middle Verde River 

watersheds, central Arizona: U.S. Geological Survey Scientific 

Investigations Report 2005–5198, 102 p., 3 plates, accessed 

July 14, 2010, at http://pubs.usgs.gov/sir/2005/5198/.

Harbaugh, A.W., 2005, MODFLOW-2005, the U.S. Geological Survey modular ground-water model -- the GroundWater Flow Process: U.S. Geological Survey Techniques 

and Methods 6-A16, variously p. accessed July 14, 2010, at 

http://pubs.usgs.gov/tm/2005/tm6A16/.

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

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

modular ground-water model—user guide to modularization 

concepts and the ground-water flow process: U.S. Geological Survey Open-File Report 00–92, 121 p. accessed July 

14, 2010, at http://water.usgs.gov/nrp/gwsoftware/modflow

2000/ofr00-92.pdf.

Leake, S.A., Pool, D.R., and Leenhouts, J.M., 2008, Simulated 

effects of ground-water withdrawals and artificial recharge 

on discharge to streams, springs, and riparian vegetation 

in the Sierra Vista Subwatershed of the Upper San Pedro 

Basin, southeastern Arizona: U.S. Geological Survey Scientific Investigations Report 2008-5207, 14 p. accessed July 

14, 2010, at http://pubs.usgs.gov/sir/2008/5207/.

Leake, S.A., and Reeves, H.W., 2008, Use of models to map 

potential capture of surface water by ground-water withdrawals: Golden, Colorado, International Groundwater 

Modeling Center, Proceedings of MODFLOW and More 

2008, p. 204–208.

Leake, S.A., Reeves, H.W., and Dickinson, J.E., 2010, A 

new capture fraction method to map how pumpage affects 

surface water flow: Ground Water, doi: 10.1111/j.1745-

6584.2010.00701.x, 11 p.

Pool, D. R., Blasch, K.W., Callegary, J., and Glaser, L., in press, 

Groundwater-flow model of the Redwall-Muav, Coconino, 

and Alluvial Basin aquifer systems of northern and central 

Arizona: U.S. Geological Survey Scientific Investigations 

Report 2010-5180.

Prudic, D.E., 1989, Documentation of a computer program to 

simulate stream-aquifer relations using a modular, finite-difference, ground-water flow model: U.S. Geological Survey 

Open-File Report 88-729, 113 p.

Stonestrom, D.A., Constantz, J., Ferré, T.P.A., and Leake, 

S.A., eds., 2007, Ground-water recharge in the arid and 

semiarid southwestern United States: U.S. Geological Survey Professional Paper 1703, 414 p. accessed July 14, 2010, 

at http://pubs.usgs.gov/pp/pp1703/.

Webb, R.H., Leake, S.A., and Turner, R.M, 2007, The ribbon 

of green—change in riparian vegetation in the Southwestern 

United States: Tucson, University of Arizona Press, 480 p.

Winter, T.C., Harvey, J.W., Franke, O.L., and Alley, W.M., 

1998, Ground water and surface water—a single resource: 

U.S. Geological Survey Circular 1139, 79 p. accessed July 

14, 2010, at http://pubs.usgs.gov/circ/circ1139/.
References Cited 19

Produced in the Western Region, Menlo Park, California

Manuscript approved for publication 

Text edited by Jan Zigler

Layout and design by Judy Weathers
20 Simulated Effects of Groundwater Pumping and Artificial Recharge in the Verde Valley Sub-Basin, Central Arizona

Leake and Pool—Simulated Effects of Groundwater Pumping and Artificial Recharge, Verde Valley Sub-Basin, Central Arizona—Scientific Investigations Report 2010–5147