Simulated Ground-Water Flow in the Hueco

Bolson, an Alluvial-Basin Aquifer System

near El Paso, Texas

Prepared in cooperation with

EL PASO WATER UTILITIES

and the

U.S. ARMY- FORT BLISS

U.S. DEPARTMENT OF THE INTERIOR

U.S. GEOLOGICAL SURVEY

Water-Resources Investigations Report 02-4108
U.S. GEOLOGICAL SURVEY

Water-Resources Investigations Report 02-4108

Simulated Ground-Water Flow in the Hueco

Bolson, an Alluvial-Basin Aquifer System

near El Paso, Texas

Prepared in cooperation with

EL PASO WATER UTILITIES

and the

U.S. ARMY - FORT BLISS

By Charles E. Heywood and Richard M. Yager

Albuquerque, New Mexico

2003
U.S. DEPARTMENT OF THE INTERIOR

GALE A. NORTON, Secretary

U.S. GEOLOGICAL SURVEY

Charles G. Groat, Director

The use of firm, trade, and brand names in this report is for identification purposes only and does

not constitute endorsement by the U.S. Geological Survey.

For additional information write to: Copies of this report can be purchased

U.S. Geological Survey

Information Services

Box 25286

Denver, CO 80225-0286

District Chief

U.S. Geological Survey

Water Resources Division

5338 Montgomery Blvd. NE, Suite 400

Albuquerque, NM 87109-1311

from:

Information regarding research and data-collection programs of the U.S. Geological Survey is available

on the Internet via the World Wide Web. You may connect to the home page for the New Mexico District

Office using the URL http://nm.water.usgs.gov.
iii

CONTENTS

Page

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

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

Purpose and scope ....................................................................................................................................................... 2

Previous investigations ................................................................................................................................................ 2

Hydrogeology of the Hueco Bolson ...................................................................................................................................... 2

Extent and thickness of hydrogeologic facies ............................................................................................................. 5

Hydraulic conductivities.............................................................................................................................................. 5

Ground-water levels..................................................................................................................................................... 5

Rio Grande Valley surface-water system..................................................................................................................... 8

Agricultural canals and drains ..................................................................................................................................... 12

Recharge ...................................................................................................................................................................... 12

Mountain-front recharge.................................................................................................................................... 12

Underflow from Tularosa and Mesilla Basins ................................................................................................... 13

Human-induced recharge................................................................................................................................... 13

Discharge ..................................................................................................................................................................... 13

Ground-water withdrawals ................................................................................................................................ 13

Evapotranspiration............................................................................................................................................. 13

Seepage to Rio Grande and agricultural drains ................................................................................................. 13

Land subsidence .......................................................................................................................................................... 14

Steady-state and transient ground-water flow model............................................................................................................. 14

Numerical method ....................................................................................................................................................... 15

Aquifer dewatering ............................................................................................................................................ 15

Multi-aquifer wells ............................................................................................................................................ 15

Spatial discretization.................................................................................................................................................... 18

Vertical datums ............................................................................................................................................................ 18

Temporal discretization ............................................................................................................................................... 18

Boundary conditions.................................................................................................................................................... 19

Specified-flow boundaries ................................................................................................................................. 19

Head-dependent flow boundaries ...................................................................................................................... 19

Streamflow routing ............................................................................................................................................ 20

Faults............................................................................................................................................................................ 23

Storage properties ........................................................................................................................................................ 23

Model calibration................................................................................................................................................................... 23

Hydraulic-head measurements .................................................................................................................................... 23

Flow measurements ..................................................................................................................................................... 23

Parameter estimation ................................................................................................................................................... 23

Model evaluation and simulation results ............................................................................................................................... 26

Estimates of aquifer properties .................................................................................................................................... 26

Parameter sensitivities ................................................................................................................................................. 26

Water-level drawdowns................................................................................................................................................ 29

Ground-water budget ................................................................................................................................................... 29

Seepage from the Rio Grande...................................................................................................................................... 29

Summary and conclusions..................................................................................................................................................... 51

Selected references ................................................................................................................................................................ 53

Appendix 1: Modifications to MODFLOWP ........................................................................................................................ 56

Appendix 2: Modifications to MODFLOW .......................................................................................................................... 63

Appendix 3: Multi-aquifer well package............................................................................................................................... 68
iv

FIGURES

Page

1. Shaded-relief map showing principal physiographic features and thickness of alluvial deposits in the

Hueco Bolson..................................................................................................................................................... 3

2. Shaded-relief map showing ground-water flow-model domain, location of wells, and trace of geologic

sections shown in figure 3.................................................................................................................................. 4

3. Generalized geologic sections A-A′ along the western margin of the Hueco Bolson and B-B′ along an

east-west transect north of El Paso, Texas ......................................................................................................... 6

4. Cutaway perspective view showing generalized horizontal hydraulic conductivity of hydrogeologic facies......... 7

5. Graph showing ground-water withdrawals from the Hueco Bolson, 1903-96........................................................ 8

6. Hydrographs showing water-level declines in selected wells in and near El Paso and Ciudad Juarez, 1935-96.... 9

7. Graphs showing increases in chloride concentration in selected wells, 1951-96.................................................... 10

8. Schematic showing stream segments used for routing surface-water flows............................................................ 11

9. Schematic diagrams showing modifications to MODFLOWP and MODFLOW:

A. Computation of head in partially dewatered, multilayer observation wells ....................................................... 16

B. Stream leakage through dewatered model cell ................................................................................................... 16

C. Distribution of hydraulic head and flow near multilayer pumped well .............................................................. 17

10. Graphs showing mean annual and mean monthly flows in Rio Grande and diversions, 1903-96 .......................... 22

11. Map showing locations of head measurements, stream segments, agricultural drains, and model faults............... 24

12. Graph showing number of annual head measurements in wells included in the nonlinear regression used

in the ground-water flow model......................................................................................................................... 25

13. Graphs showing simulated and measured heads and weighted residuals computed in transient-state

simulation of the Hueco Bolson......................................................................................................................... 27

14. Maps showing simulated water table in shallow aquifer (model layers 1 and 2) in:

A. 1902 (steady-state conditions)............................................................................................................................ 30

B. 1958 .................................................................................................................................................................... 31

C. 1973 .................................................................................................................................................................... 32

D. 1980 .................................................................................................................................................................... 33

15. Maps showing simulated potentiometric surface in model layer 5 in:

A. 1958 .................................................................................................................................................................... 34

B. 1973 .................................................................................................................................................................... 35

C. 1980 .................................................................................................................................................................... 36

D. 1996 .................................................................................................................................................................... 37

16. Perspective views showing water-level drawdown in the Hueco Bolson aquifer in:

A. 1973 .................................................................................................................................................................... 38

B. 1996 .................................................................................................................................................................... 39

17. Hydrographs showing measured and simulated water levels from 1930 through 1995 in selected wells in the:

A. Mesa-Nevins well field....................................................................................................................................... 40

B. Airport well field ................................................................................................................................................ 42

C. Cielo Vista and Town well fields ........................................................................................................................ 43

D. Water-plant well field ......................................................................................................................................... 44

E. Juarez well field.................................................................................................................................................. 45

F. East observation well field .................................................................................................................................. 47

G. Lower valley and north observation well fields.................................................................................................. 48

18. Graphs showing simulated annual inflows and outflows to Hueco Bolson aquifer:

A. Ground-water withdrawals. ................................................................................................................................ 49

B. Storage................................................................................................................................................................ 49

C. Stream leakage.................................................................................................................................................... 49

D. Agricultural drains and evapotranspiration ........................................................................................................ 49

19. Graphs showing Rio Grande seepage between Chamizal zone and Riverside Dam ...................................... 50

20. Graph showing simulated and measured flows at Riverside Dam ................................................................ 52
v

TABLES

1. Measured and simulated flow loss from unlined section of the Rio Grande above Riverside Dam

and Franklin Canal............................................................................................................................................. 12

2. Conversions between vertical datums in the El Paso, Texas, area........................................................................... 18

3. Optimum parameter values estimated for Hueco Bolson aquifer, through nonlinear regression in transientstate simulation, and their approximate 95-percent confidence intervals .......................................................... 28

CONVERSION FACTORS, ABBREVIATIONS, AND DATUMS

Altitude, as used in this report, refers to distance above sea level.

Sea level: In this report, "sea level" refers to the National Geodetic Vertical Datum of 1929—a

geodetic datum derived from a general adjustment of the first-order level nets of the United States

and Canada, formerly called Sea Level Datum of 1929.

Multiply By To obtain

centimeter (cm) 0.06102 inch (in.)

millimeter 0.03937 inch (in.)

meter (m) 3.281 foot (ft)

kilometer (km) 0.6214 mile (mi)

square kilometer (km2

) 0.3861 square mile (mi2)

square kilometer (km2

) 247.1 acre

cubic meter (m3) 0.0008107 acre-foot (acre-ft)

cubic meter (m3

) 35.31 cubic foot (ft3)

cubic meter (m3) 264.2 gallon (gal)

cubic meter per day (m3

/d) 0.0004087 cubic foot per second (ft3/s)

cubic meter per day (m3

/d) 0.2961 acre-foot per year (acre-ft/yr

cubic meter per day (m3

/d) 96,489 gallons per year (gal/yr)

acre-foot (acre-ft) 1,233 cubic meter (m3

)
1

ABSTRACT

The neighboring cities of El Paso, Texas, and

Ciudad Juarez, Chihuahua, Mexico, have historically

relied on ground-water withdrawals from the Hueco

Bolson, an alluvial-aquifer system, to supply water to

their growing populations. By 1996, ground-water

drawdown exceeded 60 meters in some areas under

Ciudad Juarez and El Paso.

A simulation of steady-state and transient

ground-water flow in the Hueco Bolson in westernmost

Texas, south-central New Mexico, and northern

Chihuahua, Mexico, was developed using

MODFLOW-96. The model is needed by El Paso

Water Utilities to evaluate strategies for obtaining the

most beneficial use of the Hueco Bolson aquifer

system. The transient simulation represents a period of

100 years beginning in 1903 and ending in 2002. The

period 1903 through 1968 was represented with 66

annual stress periods, and the period 1969 through

2002 was represented with 408 monthly stress periods.

The ground-water flow model was calibrated

using MODFLOWP and UCODE. Parameter values

representing aquifer properties and boundary

conditions were adjusted through nonlinear regression

in a transient-state simulation with 96 annual time steps

to produce a model that approximated (1) 4,352 water

levels measured in 292 wells from 1912 to 1995, (2)

three seepage-loss rates from a reach of the Rio Grande

during periods from 1979 to 1981, (3) three seepageloss rates from a reach of the Franklin Canal during

periods from 1990 to 1992, and (4) 24 seepage rates

into irrigation drains from 1961 to 1983. Once a

calibrated model was obtained with MODFLOWP and

UCODE, the optimal parameter set was used to create

an equivalent MODFLOW-96 simulation with monthly

temporal discretization to improve computations of

seepage from the Rio Grande and to define the flow

field for a chloride-transport simulation.

Model boundary conditions were modified at

appropriate times during the simulation to represent

changes in well pumpage, drainage of agricultural

fields, and channel modifications of the Rio Grande.

The model input was generated from geographic

information system databases, which facilitated rapid

model construction and enabled testing of several

conceptualizations of hydrogeologic facies boundaries.

Specific yield of unconfined layers and hydraulic

conductance of Quaternary faults in the fluvial facies

were the most sensitive model parameters, suggesting

that ground-water flow is impeded across the fault

planes.

INTRODUCTION

The neighboring cities of El Paso, Texas, and

Ciudad Juarez, Chihuahua, Mexico, have historically

relied on ground-water withdrawals from the Hueco

Bolson, an alluvial-aquifer system, to supply water to

their growing populations. At the end of the 20th

century, about 680,000 people lived in the El Paso

metropolitan area (U.S. Department of Commerce,

2002) and the total population (including Ciudad

Juarez) was about 2 million. The term “bolson,”

literally meaning “handbag” or “purse” in Spanish, is

local terminology that may be considered synonymous

with “basin.” The Hueco Bolson is considered the

southern portion of the Tularosa-Hueco Basin. The

term “Tularosa Basin” is used for the northern portion,

which lies entirely in the State of New Mexico and is

not modeled in this study.

In the United States, diversions from the Rio

Grande and ground-water withdrawals from the

Mesilla Basin (west of the study area) also supply the

freshwater demands of the military, industries, and

public in the El Paso area. In Mexico, diversions from

the Rio Grande are used for agriculture; water needed

by Ciudad Juarez is supplied solely by extraction from

the Hueco Bolson. By 1996, ground-water drawdown

exceeded 60 m in some areas under Ciudad Juarez and

El Paso. Fresh ground water stored in the aquifer

system beneath these cities is bordered by regions of

brackish to saline ground water. As water levels in the

freshwater portions of the aquifer declined, intrusion of

the surrounding brackish water degraded water quality

SIMULATED GROUND-WATER FLOW IN THE HUECO

BOLSON, AN ALLUVIAL-BASIN AQUIFER SYSTEM NEAR

EL PASO, TEXAS

By Charles E. Heywood and Richard M. Yager
2

in public supply wells, which sometimes required well

abandonment.

Prior to human intervention, infiltration of water

from the Rio Grande was the dominant mechanism of

aquifer-system recharge in the El Paso area. During the

20th century, numerous diversions from the Rio

Grande affected the distribution of this potential

recharge water. A section of the Rio Grande between

downtown El Paso and Ciudad Juarez was converted to

a lined canal in 1968, which prevented subsequent

recharge from the river to the ground-water system

along that section. This decreased recharge

exacerbated the effect of increased ground-water

withdrawals, increasing the rate of ground-water-level

declines (Land and Armstrong, 1985).

Purpose and Scope

This report describes the hydrogeology of the

Hueco Bolson and documents a transient ground-water

flow model of the Hueco Bolson. The model,

developed in cooperation with El Paso Water Utilities

(EPWU) and the U.S. Army at Fort Bliss, is needed by

the EPWU to evaluate strategies for obtaining the most

beneficial use of the Hueco Bolson aquifer system.

Included in the report are a (1) description of the

hydrogeologic features relevant to the numerical

simulation, (2) summary of the computer codes used

for the simulations, (3) description of the modelcalibration procedure, and (4) discussion of the results

of the steady-state and transient simulations, including

estimation of seepage loss with a specified design flow

in the Rio Grande and no diversion to the American

Canal Extension (ACE). Modifications to

MODFLOWP and MODFLOW are presented in two

appendixes.

 The active model area encompasses 5,303 km2

(2,408 mi2) in Mexico and the United States, extending

from 37 km north of the New Mexico-Texas State

boundary to a point on the Rio Grande 3 km south of

Fort Hancock, Texas, 79 km southeast of El Paso.

Model boundaries represent the perimeter of Hueco

Bolson deposits throughout most of the modeled area

and coincide with the Franklin and Organ Mountains to

the west and the Hueco Mountains to the east in the

United States, and the Sierra Juarez, Sierra El Presidio,

and Sierra Guadalupe to the west and south in Mexico

(figs. 1 and 2).

Previous Investigations

Sayre and Livingston (1945) first provided a

comprehensive overview of ground-water resources in

the El Paso area. Several ground-water flow models

have been developed to investigate the effects of

pumping on water levels and salinity in the Hueco

Bolson. Leggat and Davis (1966) constructed an

electric analog model to predict ground-water

drawdowns through 1990 resulting from proposed

ground-water withdrawals. A two-layer transient

model by Meyer (1976) represented freshwater with a

dissolved-solids concentration less than 1,000

milligrams per liter in both alluvial and bolson

deposits. The model was used to estimate the total

volume of freshwater in storage and to simulate waterlevel declines resulting from planned ground-water

withdrawals from 1973 to 1991.

Lee Wilson and Associates (1985a,b; 1991)

developed a four-layer model using MODFLOW in

which layer thickness corresponded to the presumed

thickness of water-quality zones. Kernodle (1992) used

the Lee Wilson and Associates (1985a,b) model to

estimate additional elastic aquifer compaction that

might result from diverting flow in a segment of the Rio

Grande into an extension of the American Canal.

Groschen (1994) developed a four-layer model using

MODFLOW and HST3D (Kipp, 1987) to simulate the

movement of saline water in response to ground-water

withdrawal and concluded that increased salinity in

wells screened in the bolson deposits was caused by

leakage from the overlying alluvial aquifer.

HYDROGEOLOGY OF THE HUECO

BOLSON

The Hueco Bolson is a fault-bounded structural

depression associated with the Rio Grande Rift (fig. 1).

At the inception of Rio Grande rifting about 26 million

years ago (Chapin and Seager, 1975), normal faults

accommodated regional extension, resulting in

downdropped structural grabens. Igneous rocks of

Precambrian age and sedimentary rocks of Paleozoic

and Mesozoic age surround and underlie the Hueco

Bolson. Unconsolidated to poorly consolidated

deposits of Tertiary and Quaternary age consisting

primarily of gravel, sand, silt, and clay have filled the

basin. These deposits compose the alluvial-aquifer

system known as the Hueco Bolson.
3

106 ̊30 ́ 106 ̊15 ́

32 ̊

15 ́

32 ̊

00 ́

31 ̊

45 ́

31 ̊

30 ́

Figure 1. Principal physiographic features and thickness of alluvial deposits in the Hueco Bolson.

0

0 10

10

20 MILES

20 KILOMETERS

1200

1800

2400

300

900

2700

1800

1200

900

600

600

900

Rio

Presidio

Guadalupe

Sierra

Mountains

Organ

Sierra

Juarez

Sierra

El

H

ueco Mountains

NM

TEXAS

NEW MEXICO

MEXICO

Franklin Mountains

Grande

El Paso

Ciudad

Juarez

EXPLANATION

300 Line of equal thickness of

alluvium. Interval 300 meters

Shaded-relief map from U.S. Geological Survey

30-meter Digital Elevation Models

Tularosa

Basin

White

Sands

Missile

Range

The

Narrows

New Mexico

Hueco

Bolson

Mexico Texas

Rio Grande Rift

106 ̊00 ́
4

Figure 2. Ground-water flow-model domain, location of wells, and trace of geologic sections shown in figure 3.

0 5 10 15 20 25 30 35 40 MILES

0 5 10 15 20 25 30 35 40 KILOMETERS

106 ̊30 ́ 106 ̊15 ́ 106 ̊00 ́

32 ̊

15 ́

32 ̊

00 ́

31 ̊

45 ́

31 ̊

30 ́

A

B ́

A ́

B

EXPLANATION

Trace of

geologic section

Model boundary

Rivers and canals

Wells

Franklin Canal

Grande

Rio

Acequia Madre

Riverside

Canal

A A ́

El Paso

Shaded-relief map from

U.S. Geological Survey

30-meter Digital

Elevation Models
5

Extent and Thickness of Hydrogeologic

Facies

The alluvial deposits that make up the Hueco

Bolson can be classified into four hydrogeologic facies

(collectively, the bolson-fill facies) on the basis of their

depositional processes and resulting sedimentary

structures:

(1) Fluvial facies. From 3.8 million to 0.67

million years ago, the ancestral Rio Grande meandered

south along the east side of the Franklin Mountains,

depositing a thick sequence of fluvial sediments

consisting of fine- to coarse-grained channel sand

interbedded with silt and clay overbank deposits. The

predominant geologic formation in this facies is the

Camp Rice Formation (Strain, 1969) of Tertiary and

Quaternary age. Electric logs of 101 wells in the El

Paso area indicate that the fraction of clay interbeds

within the freshwater portion of this facies is

approximately one-third.

(2) Alluvial-fan facies. Alluvial fans originating

from the present-day Organ and Franklin Mountains

and Sierra Juarez consist of poorly sorted gravel and

coarse- to fine-grained sand. Deposits of this facies

interfinger with the fluvial deposits of the Rio Grande.

(3) Lacustrine-playa facies. Thick deposits of

clay and silt exist in the east and southeast parts of the

Hueco Bolson and at depth beneath the fluvial and

alluvial-fan facies. These fine-grained sediments were

deposited in a low-energy environment, possibly a lake

of mid-Cenozoic age that formed a terminal depocenter

for the ancestral Rio Grande (Strain, 1969). The

predominant geologic formation of this facies is the

Tertiary Fort Hancock Formation.

(4) Recent alluvial facies. Deposition of the

fluvial, alluvial-fan, and lacustrine-playa facies and

subsequent erosion resulted in formation of the

topographic mesa that today is east and north of El Paso

(Langford, 2001). About 0.67 million years ago, the

Rio Grande breached “The Narrows” (fig. 1), the gap

separating the present-day Franklin Mountains from

the Sierra Juarez. The Rio Grande eroded the

topographic mesa, forming the present-day Rio Grande

Valley. Approximately 30 to 60 m of late Pleistocene to

recent sediments associated with the modern Rio

Grande have been deposited in the Rio Grande Valley.

The cumulative thickness of these alluvial

deposits in the Hueco Bolson is mapped in figure 1.

The distribution of these facies is illustrated in the

generalized geologic sections in figure 3 and in a

cutaway perspective view in figure 4. The horizontal

boundaries between the fluvial and lacustrine-playa

facies are believed to interfinger, resulting in

gradational changes in hydraulic conductivity at the

scale of the ground-water flow model.

Hydraulic Conductivities

EPWU has conducted aquifer pumping tests in

85 production wells. Production wells are generally

installed in known high-permeability areas, such as the

alluvial-fan and fluvial facies in the Hueco Bolson. The

average horizontal hydraulic conductivity estimated

from these tests is 10 m/d, with a standard deviation of

7 m/d. The minimum and maximum measured

horizontal hydraulic conductivities from these tests

were 1 and 50 m/d, respectively.

Laboratory measurements of permeability and

compressibility were made on clay core samples from

five different depth intervals at a site near the Rio

Grande (Harold Olsen, Colorado School of Mines,

written commun., 1995). Vertical hydraulic

conductivities of the undisturbed samples ranged from

6 x 10-3 to 2 x 10-2 m/d.

Ground-Water Levels

Extensive ground-water pumping from the

Hueco Bolson since the 1940’s has resulted in cones of

depression in the water table under El Paso and Ciudad

Juarez. The observed drawdown cones generally

correspond with the extent of the major production

wells shown in figure 2. Total ground-water pumpage

from 1903 through 1996 from the Hueco Bolson by the

United States and Mexico is shown in figure 5. Groundwater withdrawals started to increase substantially

during the 1950’s drought years and have increased in

Mexico since the early 1970’s. Hydrographs of water

levels measured in the United States and Mexico from

1935 to 1996 illustrate that ground-water levels have

declined concurrently with pumpage. Examples of

declining water levels in selected wells in the United

States are shown in figure 6. Water quality, particularly

in regard to chloride concentration, has degraded in

some areas of ground-water development from 1951 to

1996 (fig. 7), such as in well 86 near the El Paso airport

and well 148 in the Rio Grande Valley. In these wells,

ground-water pumpage may have induced intrusion

from adjacent or overlying regions of brackish water.
6

1995 water table

Layers

1 to 9

Figure 3. Generalized geologic sections A-A ́ along the western margin of the Hueco Bolson and

B-B ́ along an east-west transect north of El Paso, Texas (trace of sections shown in

figure 2).

Layers

1 to 9

Layer

10

Bend in

section

Bend in

section

A A'

600

800

1,000

1,200

Distance, in meters

0 10,000 20,000 30,000 40,000 50,000 60,000 Rio Grande

Altitude, in meters

Land surface

SOUTH NORTH

Vertical exaggeration = 33.9:1Bedrock

Bedrock

Model layer 1

2

3

4

5

6

7

8

9

10

Layer

10

Alluvial-fan facies

Fluvial facies

Lacustrine-playa facies

Recent alluvial facies

600

800

1,000

1,200

0 5,000 10,000 15,000 20,000 25,000

Distance, in meters

Bend in

section

Altitude, in meters

B B'

EXPLANATION

Land surface

Datum is sea level

WEST EAST

Vertical exaggeration = 16.6:1

Bedrock

Bedrock

Model layer 12

3

7

8

9

10

4

5

6

1903 water table 1995 water table

1903 water table
7

Figure 4. Cutaway perspective view showing generalized horizontal hydraulic conductivity of hydrogeologic

facies (logarithmic color scale bar). Area of block is identical to that shown in figure 2.

Horizontal hydraulic

conductivity,

in meters per day

6.8

5.4

3.9

2.4

0.92

NORTH
8

Rio Grande Valley Surface-Water System

The Rio Grande is hydraulically connected to the

Hueco Bolson aquifer system. Locally, water from the

river seeps into the shallow part of the aquifer system

in the Rio Grande Valley. Much of this water may be

transpired by agricultural crops and natural vegetation

in the valley. Seepage from the Rio Grande between El

Paso and Ciudad Juarez decreased substantially after

December 1968, when flow was diverted into the

concrete-lined American Canal upstream from the

Chamizal zone.

The Chamizal (named after the desert shrub

Chamiza) is a zone adjacent to the Rio Grande along

the border between El Paso and Ciudad Juarez.

Originally part of Mexico, the Chamizal became part of

the United States when the border-defining Rio Grande

changed course following a flood in 1864. In 1968, the

Chamizal returned to Mexican ownership when the

border was redefined along the newly constructed

American Canal. At that time, the reach of the natural

river channel south of the Chamizal was abandoned,

and flow in the Rio Grande was diverted into the

American Canal. The rate of ground-water-level

decline near the American Canal accelerated following

this diversion (Land and Armstrong, 1985), suggesting

that decreased local aquifer recharge from the Rio

Grande exacerbated these ground-water drawdowns.

The ACE, completed in 1999, conveys water that

formerly flowed in a 15-km (9.3-mi) reach of the Rio

Grande channel between the Chamizal zone and

Riverside Dam (fig. 8). The canal was constructed, in

part, to salvage water “lost” to seepage through the

riverbed to the underlying shallow aquifer. This

seepage was estimated for 1981 through 1983 by the

International Boundary and Water Commission

(IBWC) (Land and Armstrong, 1985; White and

others, 1997). These seepage losses are summarized in

table 1; the average loss for the 3 years was 1.03 x 105

m3/d (83.3 acre-ft/d). The magnitude of this seepage

varied as the hydraulic gradient changed between the

Rio Grande and underlying aquifer. Ground-water

levels generally declined under this reach through the

1980’s and 1990’s; thus, seepage from the Rio Grande

probably increased during this period. The magnitude

of this seepage (and by implication, the quantity of

water “salvaged”) is of interest to water management

parties in the El Paso area.

Figure 5. Ground-water withdrawals from the Hueco Bolson, 1903-96.

Pumpage, in millions of cubic meters

0

50

100

150

200

250

1903

1906190919121915

1918192119241927

1930193319361939

1942194519481951

1954195719601963

1966196919721975

1978198119841987

199019931996

Total pumpage in model

Reported by Mexico
9

W a t e r l e v e l , i n m e t e r s a b o v e s e a l e v e l

JL-49-13-610 (well 86)

1,075

1,085

1,095

1,105

1,115

JMAS 28 (well 207)

1,075

1,085

1,095

1,105

1,115

1935 1945 1955 1965 1975 1985

JL-49-13-804 (well 104)

1,075

1,085

1,095

1,105

1,115

Figure 6. Water-level declines in selected wells in and near El Paso and Ciudad Juarez, 1935-96

(location of wells shown in figure 11).
10

JL-49-13-610 (well 86)

0

100

200

300

400

500

600

700

800

900

1,000

1,100

1950 1960 1970 1980 1990

JL-49-21-301 (well 148)

0

100

200

300

400

500

600

700

800

900

1,000

1,100

1950 1960 1970 1980 1990

C h l o r i d e c o n c e n t r a t i o n , i n p a r t s p e r m i l l i o n

Figure 7. Increases in chloride concentration in selected wells, 1951-96 (location of wells shown in figure 11).
11

Figure 8. Schematic showing stream segments used for routing surface-water flows.

American

Dam

Canal

American

Mexico

International

Dam

Settling

Basin

wasteway

Leon

wasteway

Ascaratewasteway

Rio Grande

Haskell Waste Water

Treatment Plant

Canal

Water

Treatment

Plant

Riverside

Dam

Riverside

Canal

Acequia Madre

Surface-water

accounting components

Chamizal zone

Franklin Canal

Ame

io

rican

Ca x nal e tens

n Flow measuring station

Extensometer siteDirection of flow

NOT TO SCALE
12

Agricultural Canals and Drains

Because the American Canal and its extension

are concrete lined, they do not substantially interact

with the shallow ground-water system. The Franklin

Canal in the United States and the Acequia Madre in

Mexico are major unlined irrigation-supply canals.

These canals supply water through numerous

subsidiary unlined irrigation canals to agricultural

fields in the lower Rio Grande Valley.

To control shallow ground-water levels and

prevent soil salinization, agricultural drains were

installed in the Rio Grande Valley near El Paso

beginning in the 1930’s. These drains originally were

designed and constructed with constant gradient; the

Bureau of Reclamation (BOR) surveyed these drainbed altitudes in the early 1960’s. Since the 1960’s,

some drains have been destroyed by development, such

as the construction of the ACE. From 1960 through

2000, the original gradient of many drains was not

maintained (Al Blair, AWBlair Engineering, oral

commun., June 2001). Infilling from windblown sand

and lateral collapse increased drain-bed altitudes in

segments of the lower valley drains. This may have

decreased drainage efficiency, leading to the observed

higher shallow ground-water levels and increased soil

salinity in some lower valley agricultural fields.

Recharge

Precipitation over the Hueco Bolson is highly

variable, both spatially and temporally. Mean annual

precipitation over the Hueco Bolson is less than 25 cm

(10 in.), most of which falls during the summer months.

Sparse rainfall over the basin floor outside the Rio

Grande Valley probably evaporates or transpires from

the vadose zone before it can infiltrate to water-table

depths and recharge the aquifer system. In the Rio

Grande Valley, where the water table can be within

several meters of land surface, precipitation or applied

irrigation water has a better chance of infiltrating to the

water table. Concentrated surface flows in arroyos

below mountain canyons may infiltrate sufficiently to

penetrate the vadose zone to the water table.

Ground-water age, defined as the time since

water was in contact with the atmosphere, was

estimated using carbon-14 age-dating techniques at

eight different sites in the Hueco Bolson (Anderholm

and Heywood, 2003). The calculated age of water in

these samples ranges from 12,100 to 25,500 years old.

The dates indicate that this water advected from

recharge locations for some distance through the

saturated zone.

Mountain-Front Recharge

By assuming that 25 percent of precipitation that

falls in the catchments of the Organ and Franklin

Mountains may bypass caliche layers at mountain

fronts, Sayre and Livingston (1945) estimated that

recharge to the bolson-fill aquifer from these areas may

be as much as 50,000 m3/d (15,000 acre-ft/yr). Meyer

(1976) estimated total recharge from the Organ and

Franklin Mountains, the Sierra Juarez, and underflow

from the Tularosa Basin to be about 19,000 m3/d

(5,640 acre-ft/yr). Wilkins (1998) estimated recharge

along the western boundary of the Tularosa-Hueco

Basin to be about 245 m3

/d/km (0.161 ft3/s/mi). For the

boundary length corresponding to the Organ and

Franklin Mountains and Sierra Juarez in this study

(approximately 80 km), this recharge estimate equates

to about 19,600 m3

/d. Waltemeyer (2001) used the

basin-climatic characteristics method to estimate

streamflow available for potential recharge at the

mouth of two canyons at the base of the Organ

Mountains, which are in the area of this model. Total

mean annual streamflow from the Oak and Soledad

Canyon drainages was estimated to be 2,300 m3/d;

Table 1. Measured and simulated flow loss from unlined

section of the Rio Grande above Riverside Dam

and Franklin Canal

[IBWC, International Boundary and Water Commission;

BOR, Bureau of Reclamation; USGS, U.S. Geological

Survey; m3/d, cubic meters per day]

Measurement

Measured

flow loss

Simulated

flow loss

Rio Grande 1981

(IBWC)

8.4 x 104 m3/d 9.64 x 104 m3/d

Rio Grande 1982

(IBWC)

1.06 x 105 m3/d 9.82 x 104 m3/d

Rio Grande 1983

(IBWC)

1.18 x 105 m3/d 1.02 x 105 m3/d

Franklin Canal 1984

(BOR)

1.05 x 104 m3/d 9.81 x 103 m3/d

Franklin Canal 1990

(USGS)

5.06 x 104 m3/d 7.79 x 103 m3/d

Franklin Canal 1991

(USGS)

5.42 x 104 m3/d 1.23 x 104 m3/d

Franklin Canal 1992

(USGS)

5.33 x 104 m3/d 1.32 x 104 m3/d
13

some fraction of this flow may contribute to aquifer

recharge.

Underflow from Tularosa and Mesilla Basins

Precipitation falling on the Sacramento and San

Andres Mountains contributes to streamflow that

infiltrates the basin alluvium shortly after flowing from

mountain canyons onto the Tularosa Basin floor. Some

of this infiltration water probably recharges the

regional ground-water system. Ground-water levels in

the Tularosa Basin (McLean, 1970) indicate that

regional flow is to the south into the Hueco Bolson in

Texas and to Lake Lucero, near White Sands Missile

Range (fig. 1). Ground-water flow from the Tularosa

Basin into the study area is probably a major

component of recharge to the Hueco Bolson.

Ground-water underflow from the Mesilla Basin

to the Hueco Bolson may occur adjacent to the Rio

Grande through The Narrows. Slichter (1905)

determined the alluvial thickness to be less than 26 m

in this area and estimated underflow to be less than 270

m3/d (80 acre-ft/yr).

Human-Induced Recharge

Land-use zoning maps indicate that agricultural

acreage in the United States part of the Rio Grande

Valley within the model area is approximately 227 km2

(56,000 acres). About 1.2 m (4 ft) of irrigation water is

applied to this acreage per year, of which 20 percent, or

about 0.25 m (0.8 ft/yr), is estimated to infiltrate to the

water table (Al Blair, oral commun., 2000). Much of

this irrigation-return flow is intercepted by agricultural

drains and returned to the Rio Grande.

Seepage from the Rio Grande is a major recharge

component to the surrounding shallow aquifer system.

Much of this water probably is subsequently

discharged by evapotranspiration (ET) from the

relatively shallow water table in the El Paso Valley.

From 1948 to 1952, EPWU tested injection of

water into the Hueco Bolson aquifer system through a

well in El Paso Valley (Roger Sperka, El Paso Water

Utilities, oral commun., 1998). From 1971 through

1977, several valley wells were injected at a total

average rate of about 700 m3

/d (200 acre-ft/yr).

Beginning in 1981 and continuing through the 1990’s,

a major artificial recharge project north of El Paso

injected 10 recharge wells at rates as high as 18,000

m3/d (5,300 acre-ft/yr).

Discharge

Discharge from the deeper portions of the Hueco

Bolson under valley and mesa areas is principally from

ground-water pumping. Ground-water discharge

mechanisms from recent alluvial deposits in the Rio

Grande Valley are principally ET of infiltrated river

water and applied irrigation water and seepage to

agricultural drains and the Rio Grande.

Ground-Water Withdrawals

Records of historical pumpage from all known

municipal supply, military, industrial, and private wells

were compiled as part of this study. For wells in

Mexico, annual pumpage summaries by well from

1926 through 1995 were obtained from the Junta

Municipal de Agua y Sanimiento through the IBWC.

Monthly summaries of pumpage by well for EPWU

production wells and some military and industrial wells

after 1967 were available.

Historical withdrawals from the Hueco Bolson,

compiled from annual well pumpage records, are

depicted in figure 5. Beginning in the 1950’s, groundwater pumpage accelerated and reached a maximum of

about 230 million m3

/yr (186,000 acre-ft/yr) by the late

1980’s. The rapid growth of Ciudad Juarez since the

1970’s is mirrored by increased ground-water

withdrawals in Mexico. Ground-water pumping in the

United States decreased during the 1990’s as the City

of El Paso began to use more treated water from the Rio

Grande and ground water from the Mesilla Basin.

Evapotranspiration

The measured pan-evaporation rate from 1950 to

1980 at Ysleta Yard in El Paso was 4.9 mm/d (5.8 ft/yr)

(Al Blair, written commun., 2000). Although measured

pan evaporation is useful for estimating potential ET

when the water table is near land surface, actual

maximum ET rates may be higher or lower. Measured

maximum ET rates for phreatophytes near the Pecos

River in New Mexico (Weeks and others, 1987) and the

Gila River in Arizona (Culler and others, 1982) ranged

from 1 to 4 mm/d. Phreatophyte roots may extend to 9

m deep for cottonwood and greater for mesquite

(Robinson, 1958).

Seepage to Rio Grande and Agricultural Drains

When shallow ground-water levels are higher

than the water level in the adjacent Rio Grande, shallow
14

ground water seeps into the Rio Grande where there is

sufficient hydraulic conductance. Such ground-water

inflow may occur southeast of El Paso, especially

where shallow ground-water levels may be high from

applied irrigation water and insufficient drainage from

irrigated fields. The magnitude of seepage to

agricultural drains can be estimated from flow

measured in these drains. The BOR maintained records

of this flow from 1960 through 1983 (Bureau of

Reclamation, El Paso Field Division, written commun.,

2000).

Land Subsidence

Elastic aquifer compression occurs when

associated aquifer pore pressure is reduced by groundwater withdrawals and generally results in relatively

small but measurable land subsidence. If ground-water

levels decline in such a way that effective stresses in

aquifer-system matrix materials exceed previous

maximum magnitudes (the “preconsolidation stress

level”), inelastic aquifer-system compaction may

occur. For comparable incremental ground-water-level

declines, inelastic aquifer-system compaction is

typically one to two orders of magnitude greater than

elastic compression and may result in substantial land

subsidence (Riley, 1998).

First-order, first-class geodetic level lines were

surveyed by the National Geodetic Survey in the Hueco

Bolson in 1953, 1981, and 1993. By assuming stable

benchmark points on bedrock in the Hueco and

Franklin Mountains, the relative change in benchmark

altitudes was determined for points in the Hueco

Bolson for 1953-81 and 1981-93. As of 1993, the

maximum measured altitude change at a benchmark

near downtown El Paso was 0.25 m (0.82 ft) (Emery

Balasz, National Geodetic Survey, written commun.,

1994). This magnitude of altitude change is consistent

with elastic compression of the aquifer matrix resulting

from measured ground-water drawdowns. Evidence

that preconsolidation levels have been exceeded has

not been documented in the Hueco Bolson.

STEADY-STATE AND TRANSIENT

GROUND-WATER FLOW MODEL

Because brackish to saline ground water is

present in the Hueco Bolson, the significance of

density-driven flow effects was evaluated before the

numerical code to be used for model development was

selected. For this purpose, the Groschen (1994) model

was used for preliminary flow-system analysis.

Running several density-dependent simulations using

HST3D tested the magnitude of density-driven groundwater flow from heterogeneous solute concentrations.

These profile models were run with identical boundary

conditions and hydraulic properties using alternate

assumptions of homogeneous or heterogeneous

ground-water density. For the brackish to moderately

saline conditions present in the Hueco aquifer system,

computed differences between the homogeneous and

heterogeneous density assumptions were noticeable for

only a test case of isotropic hydraulic conductivity. For

horizontal to vertical hydraulic-conductivity ratios

greater than 10:1 to 100:1, which are probable for

Hueco Bolson deposits, density-driven flow effects are

negligible. For the purposes of this study, the

assumption of constant ground-water density was

concluded to be an appropriate approximation of the

physics of ground-water flow.

Because of the large size of this ground-water

flow simulation and the need for effective model

calibration, three versions of the numerical model of

the Hueco Bolson were constructed:

(1) A MODFLOW model with 1 steady-state and

96 annual transient stress periods,

(2) A MODFLOWP model with 96 annual stress

periods, and

(3) A MODFLOW model with 66 annual and

408 monthly stress periods.

The MODFLOW model (version 1) was first

constructed to build a reasonable simulation of the

Hueco Bolson. An equivalent MODFLOWP model

(version 2) was subsequently constructed to enable

model calibration by the inverse method (Hill, 1992).

These two versions were frequently run with equivalent

parameter definitions and distributions to ensure model

equivalency and to provide error checking. After an

optimal parameter set had been determined with

MODFLOWP and translated to an equivalent

MODFLOW input, a MODFLOW (version 3)

simulation with refined monthly discretization was run.

The monthly MODFLOW simulation was checked for

general equivalency to the results of the annual

MODFLOW and MODFLOWP simulations. The

minor differences in hydraulic heads and budgets are

attributable to the difference in temporal discretization.
15

Numerical Method

Both MODFLOW and MODFLOWP use the

finite-difference method to solve a form of the groundwater flow equation:

(1)

where

 Kxx, Kyy, and Kzz = principal components of the

hydraulic-conductivity tensor;

h = hydraulic head;

W = source/sink;

t = time; and

Ss = specific storage.

Pumping from the Hueco Bolson aquifer since

1903 has lowered water levels by more than 50 m (197

ft) in some areas of El Paso and Juarez, dewatering part

of the modeled area that corresponds to the top two

model layers. This situation is represented in

MODFLOW by removal of dewatered (dry) cells from

the simulation and conversion of the underlying cells

from confined to unconfined conditions by increasing

the value of the storage coefficient. Simulation of

dewatering in the upper part of the Hueco Bolson

aquifer required changes to MODFLOWP and

MODFLOW, however. Both codes were further

modified by incorporating the multi-aquifer well

package (MAW) (McDonald, 1984), which computes

flow to pumped wells screened in more than one model

layer. These modifications are detailed in the

appendixes and summarized in the next two

subsections.

Aquifer Dewatering

MODFLOWP allows spatial and temporal

interpolation in the computation of hydraulic head at

specific locations in the modeled area and also

computes mean head in observation wells screened in

more than one model layer. However, the interpolation

procedure used to compute heads in these multilayer

observation wells does not support the dewatering

condition described in the previous section. The

interpolation procedure defined by Hill (1992) was

modified by omitting dry layers screened by a

multilayer observation well and using only the

remaining saturated layers in the computation of head

in the well (hw) (fig. 9A).

Neither MODFLOWP nor MODFLOW

simulates stream leakage in model cells that are dry,

although the actual stream channels represented in the

model continue to leak at a constant rate as the water

table declines. The stream package (Prudic, 1989) was

modified in both codes to allow continuing stream

leakage to the aquifer following the procedure used in

the recharge package, in which areally distributed

recharge enters the topmost active cell if the upper

layers in the model are dry (fig. 9B). This modification

also required changes to the computation of stream

leakage reported in both the water budget and flow

observations computed by MODFLOWP.

Multi-Aquifer Wells

The screened intervals of pumped wells in the

Hueco Bolson average more than 100 m; most of these

wells are screened in more than one model layer. The

MAW package (McDonald, 1984) was used to

compute discharge from each model layer (qk) to a well

using the total discharge (Qw) specified for the well.

The method used in the MAW package was described

by Bennett and others (1982) and has been previously

applied in ground-water flow models developed by

Kontis and Mandle (1988) and by Groschen (1994).

The method is based on the Thiem (1906) equation

describing steady-state radial flow to a discharging

well:

(2)

where h = hydraulic head at a distance r from the

well [L];

hw = hydraulic head at the well radius rw

[L];

Qw = well discharge [L3

T]; and

T = transmissivity [L2

T].

Prickett and Lonnquist (1971) and Trescott and

Larson (1976) used the Theim equation in a finitedifference ground-water flow model to estimate hw in a

discharging well from the head computed at a cell

representing the well. In their applications, h was

assumed to be equivalent to the head at an effective

radial distance (ra) from a well located at the center of

∂x

∂ Kxx∂x

∂h

∂y

∂ K

yy∂y

∂h

∂z

∂ Kzz∂z

∂h

+ + – W ∂t

∂hSs =

h hw – Qw

2πT ---------- r

rw

-----     = ln
16

Model

layer:

1

Water table

p2 = 0.25

p3 = 0.375

p4 = 0.375

h2

h3

h4

h3

h4

2

3

4

Before dewatering:

hw = p2h2 + p3h3 +p4h4

After dewatering:

hw = p3h3 +p4h4

Model

layer:

1

Water table

p3 = 0.5

p4 = 0.5

2

3

4

(A) Computation of head in partially dewatered, multilayer observation wells.

EXPLANATION

Water table

hstr

botstr

Stream leakage

Model

layer:

1

2

Before dewatering, leakage to model layer 1:

qstr = C(hstr -botstr)

(B) Stream leakage through dewatered model cell.

After dewatering, leakage to model layer 2:

qstr = C(hstr -botstr)

Water table

hstr

botstr

Stream leakage

Model

layer:

1

2

EXPLANATION

Figure 9. Modifications to MODFLOWP and MODFLOW.

hstr

botstr

C 

qstr

= stream stage

= altitude of streambed bottom

= streambed conductance

= flow in stream

= fraction of the well

 screened in layer k

= head in layer k

= head in well

pk

hw

hk
17

a square cell with length (∆x) (fig. 9C). Bennett and

others (1982) approximated the effective radius as:

(3)

and showed that hw in a multi-aquifer well at cell ij

screened in model layers m to n can be computed from:

(4)

where the subscripts = row i, column j, and layer k and

Qw = the algebraic sum of discharges to the well from

layers m through n. Discharge to the well from each

model layer (qwk) can then be computed from the

Theim equation (eq. 2) by substituting values for hw

and hijk. McDonald (1984) incorporated this method in

the MAW package that was written for use with

MODFLOW. The MAW package was modified for this

study to support dewatering of model layers by

omitting dry layers from the computation of hw in

equation 4 and apportioning flows qwk to or from the

well in the remaining saturated layers.

ra

∆x

4.81 = ----------

hw

Tijkhijk

rak rw ln( ) ⁄ --------------------------

k m =

n

∑

Tijk

rak rw ln( ) ⁄ --------------------------

k m =

n

∑

-------------------------------------- Qw

2π

Tijk

rak rw ln( ) ⁄ --------------------------

k m =

n

∑

= – ---------------------------------------------

(C) Distribution of hydraulic head and flow near multilayer pumped well.

qw3

qw4

qw1

Qw

qw2

Model

layer:

1

2

3

4

Qw + qw3 + qw4 = qw1 + qw2

hw

h1

h2

h3

h4

ra

EXPLANATION

Figure 9. Modifications to MODFLOWP and MODFLOW: (figure 9C modified from Bennett and others,

 1982)--Concluded.

= flow to or from the well in layer k

= effective well radius

= well pumpage

q

Q 

r a

wk

w
18

The assumption of radial flow to a single well

located in the center of a model cell limits the accuracy

of the method’s representation of certain conditions.

The assumption of radial flow is not strictly satisfied at

model cells adjacent to an impermeable boundary, so

wells in these cells are not represented in the MAW

package. Kuniansky and Hillestad (1980) showed,

however, that if a well is located in the center of a cell

adjacent to an impermeable boundary, the error in hw

using the ra value computed with equation 3 is less than

5 percent of the exact value computed with an

analytical solution.

Another potential problem is representing

closely spaced, multilayer wells with combined

pumpage that produces an actual head (hw) lower than

the value computed by the MAW package. The

combined drawdown produced by closely spaced wells

could be computed by the use of smaller model cells to

more accurately represent each well or superposition

(Reilly and others, 1984) by combining the discharge

or recharge from multiple wells into a single well

within the model cell. The head (h) computed for a

model cell in several test cases was determined to be

relatively insensitive to the choice of these alternative

methods, however, so multi-aquifer wells were

specified individually in this model.

Spatial Discretization

To analyze pumping effects of individual wells,

the finite-difference model grid was designed to

maximize spatial detail yet retain reasonable execution

time. Unconsolidated deposits ranging from 138 to 693

m thick above an altitude of 600 m are represented with

10 model layers. Each model layer consists of 165 rows

and 100 columns of cells 500 or 1,000 m on either side.

Model layer 1 contains 10,895 active cells and

represents a larger area (5,099 km2

) than deeper model

layers because the basin narrows with depth (fig. 3);

layer 10 (the bottom layer) contains 8,809 active cells

and represents an area of 3,839 km2

. To maintain an

approximately constant depth below land surface, the

model-grid altitudes of all model layers increase to the

north with a gradient of 1:1,000 to the approximate

location of the Texas-New Mexico State line, north of

which the layers are horizontal.

Recent alluvial deposits in the Rio Grande Valley

are represented in model layer 1, which has an average

thickness of 30 m (fig. 3). Along the western margin of

the basin, the fluvial facies is represented in model

layers 1 through 9, each of which is 30 m thick. The

lacustrine-playa facies is represented along the eastern

margin of the basin in model layers 1 through 9 and

throughout the model in layer 10, which varies in

thickness from 0 to 276 m.

Vertical Datums

The El Paso area has six different vertical

datums, three of which were encountered in data used

to define altitudes of various model boundaries.

Measured hydraulic heads were referenced to the

National Geodetic Vertical Datum of 1929 (NGVD29), which was also used as the vertical datum for this

study. Altitudes of the Rio Grande channel surveyed by

the IBWC were referenced to the IBWC datum; these

were converted to NGVD-29 by adding 0.348 m before

streambed altitudes in the model were defined.

Altitudes of agricultural drains surveyed by the BOR

were referenced to an old datum of the Santa Fe

Railway; these were converted to NGVD-29 by adding

12.891 m before drain-bed altitudes in the model were

defined. The differences among vertical datums used to

perform these conversions are listed in table 2.

Temporal Discretization

A steady-state simulation was used to represent

predevelopment conditions and to obtain starting heads

for the transient simulation. Sixty-six annual stress

periods and time steps were used to simulate 1903

through 1968. Because the model parameterestimation process required thousands of model runs,

using monthly stress periods during the modelcalibration process was not possible. A MODFLOWP

Table 2. Conversions between vertical datums in the

El Paso, Texas, area

Datum

Santa Fe

Railway

International

Boundary

and Water

Commission

(IBWC)

National

Geodetic

Vertical

Datum of

1929

(NGVD-29)

Santa Fe

Railway 0

+ 12.543

meters

+ 12.891

meters

IBWC

- 12.543

meters 0 + 0.348 meter

NGVD-29

-12.891

meters -0.348 meter 0
19

simulation with steady-state and 93 annual stress

periods was constructed to simulate the time period of

calibration data, which extended through 1995. The

optimal parameter set obtained with MODFLOWP was

used to generate an annual stress-period MODFLOW

data set so that both models could be run and checked

for equivalency. The 1968 head distribution from the

annual stress-period MODFLOW simulation was used

as a starting head distribution for the monthly stressperiod simulation, which represented 1969 through

2002.

Seasonal variations in Rio Grande flow, groundwater pumping, irrigation recharge, and ET cause

hydraulic-head variations that cannot be represented

with annual stress periods. To accurately incorporate

these seasonal effects, a MODFLOW simulation

(version 3) was discretized with 408 monthly stress

periods from January 1969 through December 2002.

This simulation was identical to the annual

MODFLOW (version 1) and MODFLOWP

simulations (version 2) for the first 66 annual stress

periods, which simulate 1903-68. Parameter values and

distributions, as well as annual total pumpage stresses,

were identical to those in the annual simulations.

Monthly temporal discretization enabled more

accurate calculations of Rio Grande seepage losses

along the 15-km (9.3-mi) reach between the end of the

lined section of the Chamizal zone and Riverside Dam.

Boundary Conditions

Specified-flow boundaries represent recharge

along the basin margins and beneath irrigated fields and

discharge from pumped wells screened in the fluvial

facies (figs. 3 and 4). Head-dependent boundaries

represent seepage losses from stream channels,

including the Rio Grande and several irrigation canals,

and ground-water discharge to stream channels,

agricultural drains, and ET.

Specified-Flow Boundaries

Specified-flow boundaries in model layer 1

represent infiltration from ephemeral streams draining

the mountains that border the Hueco Bolson and

leakage from irrigation-return flows. In model layers 1

through 9, specified-flow boundaries represent

underflow from upgradient areas in the Tularosa and

Mesilla Basins. Rates of recharge and underflow were

included in the estimated parameter set considered in

the nonlinear regression. Mountain-front recharge is

applied along arroyo channels on the east side of the

basin and at the base of the Organ and Franklin

Mountains in the United States and a 14-km2 (5.4-mi2)

area at the base of the Sierra Juarez in Mexico. The

estimate of mountain-front recharge (800 m3/d) is

smaller than values used in previous modeling studies

in the area (Meyer, 1976; Groschen, 1994). During

simulation stress periods after 1924, recharge also is

applied to model cells that represent irrigated fields in

the United States and Mexico. The maximum recharge

rate represents the volume of water typically applied in

excess of crop requirements in the United States (Al

Blair, oral commun., 2000). During drought and lowflow years, this rate was scaled so that applied

irrigation water would not exceed the total surface

water available. The types and acreage of irrigated

crops have changed during the past 100 years, but little

information is available to document these changes

over the entire simulation period, particularly in

Mexico. The irrigated area (approximately 180 km2

)

and maximum rate of return flow are, therefore,

considered constant in the model simulation.

Specified flows represent underflow from the

Tularosa Basin along the northern model boundary and

from the Mesilla Basin (Slichter, 1905). Ground-water

withdrawals through pumped wells are specified in

model layers 1 through 9 using 1903-96 annual

pumpage rates in 424 wells. Thirty-two wells that are

screened within only one model layer are represented

with the MODFLOW well package, and the remaining

392 wells are represented with the MAW package

described earlier. About 60 percent of the wells are

located in the United States, including 10 wells that are

used to recharge water to the bolson aquifer.

Head-Dependent Flow Boundaries

The exchange of water between stream channels

and recent alluvial sediment in the Rio Grande Valley

is simulated with the stream package in MODFLOW.

Ground-water discharge to agricultural drains is

simulated with the drain package. In both packages the

flow of water between the boundary and the underlying

model cell is a function of the head in the stream or

drain, the head in the model cell, and the hydraulic

conductance of the sediment in the channel bed:

C (5) kA

b = ------
20

where k = vertical hydraulic conductivity of the bed

sediment [LT -1];

A = area of the streambed or drain bed in the

cell [L2]; and

b = thickness of the bed sediment [L].

Geographic information system (GIS) coverages

of the Rio Grande, canals, and drain channels enabled

accurate calculation of the length of these features

within each model cell. Surveyed widths of the Rio

Grande between the Chamizal zone and Riverside Dam

(Dr. Rong Kuo, U.S. International Boundary and Water

Commission, written commun., 1999) were used to

calculate the streambed area (A) of the Rio Grande

channel for corresponding model cells. The average of

these widths (30 m) was used to calculate streambed

area for reaches where surveyed widths were

unavailable. The portion of agricultural drains through

which water seeps (wetted perimeter) is variable but

was estimated during field reconnaissance to average

about 3 m. Altitudes of drain-bed end points surveyed

by the BOR in 1960 were assigned to the associated

drain in a GIS coverage. Altitudes of intermediate

points along each individual drain were determined by

linear interpolation along the arc length between the

end points. The drain-bed altitude assigned to each

flow-model cell was the average altitude of all drains

within that cell. The reasonability of these altitudes was

verified by comparing them with cell-average values

derived from digital elevation models (DEM’s),

corrected for an assumed drain-bed depth of 2-3 m.

Because drains may have filled in with sediment since

1960, drain-bed altitudes in some model runs were

increased by 1 m during the calibration process.

Channel-bed thickness of the Rio Grande,

canals, and drains was assumed to be about 1 m.

Because the bed thickness is somewhat uncertain and

wetted perimeter widths vary, the vertical conductance

per unit length of each of these features, rather than the

vertical hydraulic conductivity, was included in the

estimated parameter set considered in the nonlinear

regression.

In the ground-water flow simulations, essentially

nonirrigated and undrained agricultural conditions

change to irrigated and drained conditions in the Rio

Grande Valley in 1925. (For ground-water seepage to

the 21 major agricultural drains simulated with the

drain package in MODFLOW and MODFLOWP, the

reader is referred to figure 18D.)

ET from nonirrigated riparian land in the Rio

Grande Valley was represented using the ET package

of MODFLOW and MODFLOWP. Land-surface

altitudes throughout the model were obtained from

28-m DEM’s. Maximum depths of ET between 2.5 and

9 m were tested during model calibration. Although the

model was not particularly sensitive to ET extinction

depth, the best model fit had an ET extinction depth of

5 m. The initial estimate (4 mm/d) of the maximum ET

rate was refined through the regression process to a

final value of 4.6 mm/d.

Streamflow Routing

As calculated with the stream package in

MODFLOW, water seepage between the Rio Grande

and the underlying aquifer is dependent on (1) the

difference in hydraulic head between the boundary

reach representing the Rio Grande and the model cell

representing the underlying aquifer and (2) the

hydraulic conductance between the two.

Boundary heads representing water levels of the

Rio Grande, Franklin Canal, Ascarate wasteway, and

Acequia Madre are computed with the Manning

equation in the stream package (Prudic, 1989).

Streambed slope, width, and roughness are spatially

variable properties specified for each boundary reach.

Stream discharge, or flow in a particular boundary

reach, varies both spatially and temporally. The slope

of each stream reach was computed from surveyed

altitudes, if available, or interpolated from digital

elevation data. Manning’s n was specified as 0.03 for

the Rio Grande channel and 0.004 for the Franklin

Canal, Ascarate wasteway, and Acequia Madre. The

time-varying flow in the first reach of stream segments

representing the Rio Grande, Franklin Canal, Ascarate

wasteway, and Acequia Madre was specified in the

stream package in MODFLOW. The stream stage in

each of these starting stream reaches was computed

from this flow using the Manning equation. Discharge

into each subsequent downstream stream reach was

adjusted by seepage loss or gain from the connected

aquifer model cell, and stream stage in each reach was

calculated with the Manning equation.

Surface-water flows in the Rio Grande,

American Canal, Franklin Canal, Ascarate wasteway,

and Acequia Madre in Mexico were simulated with

seven stream segments in the stream package (fig. 8).

Diversions from the Rio Grande flow into the Franklin

Canal and Acequia Madre. A diversion from the

Franklin Canal back to the Rio Grande flows in the

Ascarate wasteway. Flow in the Franklin Canal is

routed back to the Rio Grande through the Tornillo

Canal.

Daily records of flow in the Rio Grande at

American Dam, diversion into the Franklin Canal, and
21

diversion to Mexico at International Dam from June 1,

1938, through December 31, 1996, obtained from the

IBWC (2000) were used to compute mean annual and

mean monthly input flow specifications for stream

segments representing the Rio Grande, American

Canal, Franklin Canal, and Acequia Madre. Water

diverted from the Franklin Canal back to the Rio

Grande was specified according to measurements

obtained from the BOR (El Paso Field Division,

written commun., 2000).

Mean annual flows in the Acequia Madre,

Franklin Canal, and Rio Grande from 1903 through

1996 are shown in figure 10A. Diversions from the Rio

Grande into the Franklin Canal began in 1938. Flow in

the Rio Grande was virtually nonexistent during the

drought years of 1954-57. The late 1970’s were

relatively dry years; flow during this period is shown in

monthly detail in figure 10B. In contrast, the late

1980’s to early 1990’s were relatively wet; flow during

this period is shown in monthly detail in figure 10C.

Flow in the stream segment simulating the Rio

Grande at the end of the Chamizal zone was specified

for the monthly stress periods beginning in January

1969. Because Rio Grande flow was measured

upstream from the lined section of the American Canal,

the flow specified for the simulation required

adjustment to account for tributary inflow and loss

between the streamflow-gaging station and the start of

the unlined section of the Rio Grande. These

adjustments were made according to an analysis of the

El Paso County Water Improvement District canal

system (Al Blair, written commun., January 26, 2000).

Additions to the flow measured downstream from

American Dam were made to account for inflow from

the Settling Basin wasteway, Leon wasteway, and

Haskell Waste Water Treatment Plant. The

accumulated flow was debited to account for the

diversion to Mexico at International Dam and

estimated evaporation in the canal between American

Dam and the downstream end of the Chamizal zone.

This flow is summarized in equation form:

(6)

where Qseg = specified inflow to the stream

segment below the Chamizal

zone;

Qamer = measured flow in the American

Canal below American Dam;

Qmex = measured diversion to Mexico at

International Dam;

Qsb = flow in the Settling Basin

wasteway;

Qleon = flow in the Leon wasteway;

Qhask = return flow from the Haskell Waste

Water Treatment Plant; and

Qevap = estimated total evaporation from

the American Canal between

American Dam and the end of the

Chamizal zone.

Because tributary inflows from the Settling

Basin wasteway and Leon wasteway were not directly

measured, their inflow was derived by subtracting the

measured flow in Franklin Canal and the diversion to

the Canal Water Treatment Plant from the measured

flow in Franklin Canal downstream from American

Dam. In equation form:

(7)

where Qfrank1 = measured flow in the Franklin

Canal below American Dam;

Qwtp = measured diversion to the Canal

Water Treatment Plant; and

Qfrank2 = measured flow in the Franklin

Canal between the Leon and

Ascarate wasteways.

Monthly return flows from the Haskell Waste Water

Treatment Plant (Qhask) were relatively constant

(Roger Sperka, oral commun., June 2000), enabling

use of a constant average flow of 75,708 m3/d. The

monthly evaporative fluxes (Qevap) were generated

from local pan-evaporation rates by Al Blair (written

commun., 2000).

To specify monthly flow in the Rio Grande from

1994 through 2004, a “design flow” was synthesized

from other components of the El Paso surface-water

budget (Al Blair, written commun., 2000). The

specified flow for the stream-package segment

representing the Rio Grande downstream from the

Chamizal zone was set equal to 60 percent of the design

diversion to El Paso plus the design diversion at

Riverside plus the computed seepage loss in the 15-km

reach between the Chamizal zone and Riverside Dam.

(Because the seepage loss in this segment was

computed with the ground-water flow model, which

required this flow specification as input, several

iterations of the model were run to obtain a flow

specification.) This calculation may be summarized as:

(8)

The resulting design-flow specifications for stream

segments representing the Acequia Madre, Franklin

Canal, and Rio Grande are shown in figure 10D.

Qseg Qamer Qmex – Qsb Qleon Qhask Qevap = ++ + –

Qsb Qleon + Q frank1 Qwtp Q frank2 = – –

RioGrande = 0.6( ) diversion + + seepage Riverside ̇
22

Figure 10. Mean annual and mean monthly flows in Rio Grande and diversions, 1903-96.

S t r e a m f l o w, i n t h o u s a n d s o f c u b i c m e t e r s p e r d a y

(A) 1903-96, by year

1900 1910 1920 1930 1940 1950 1960 1970 1980 1990

0

1,000

2,000

3,000

4,000

5,000

6,000

7,000

8,000

9,000

(C) 1988-92, by month

0

1,000

2,000

3,000

4,000

5,000

1988 1989 1990 1991 1992

(D) Design flow, by month

0

1,000

2,000

3,000

4,000

5,000

Jan Feb Mar Apr May June July Aug Sept Oct Nov Dec

(B) 1975-79, by month

0

1,000

2,000

3,000

4,000

5,000

1975 1976 1977 1978 1979

Acequia Madre Franklin Canal Rio Grande

EXPLANATION
23

Faults

Fault gouge in subvertically dipping fault planes

may impede horizontal fluid flow and affect groundwater flow patterns (Haneberg and others, 1999).

Several intrabasin faults deforming Quaternary

alluvium are recognized northeast of El Paso (Barnes,

1983) (fig. 11). To evaluate the influence of these faults

on ground-water flow, they were incorporated into the

ground-water flow simulation with the Horizontal Flow

Barrier package (Hsieh, 1992).

Storage Properties

Water in storage in the aquifer system may be

released by three main processes: (1) drainage of pore

space in unconfined layers, (2) expansion of water, and

(3) compression of the aquifer matrix. Drainage of

water from pore space above a declining water table is

the dominant process yielding water to the groundwater flow system and is quantified by specific yield

(Sy). Specific yield generally ranges from 10 to 20

percent. Its magnitude in the Hueco Bolson was

estimated in the model-calibration process as a uniform

value throughout the alluvial deposits. Specific yield

and the component of specific storage that is due to the

expansion of water (which is very minor) were

simulated with the Block-Centered-Flow package in

MODFLOW and MODFLOWP. The skeletal

components of specific storage, which quantify water

derived from compression of the aquifer-system

matrix, were simulated with the Interbed-Storage

package in MODFLOW and MODFLOWP. The

magnitude of elastic skeletal specific storage (Sske) is

reasonably well constrained by plausible physical

properties of alluvial deposits. In the El Paso area,

onsite estimates of Sske (Heywood, 1995) have been

made from piezometric-extensometric measurements;

this value (Sske = 7 x 10-6 m-1) was specified in the flow

model.

MODEL CALIBRATION

MODFLOWP computes the changing spatial

distribution of ground-water heads and fluxes across

head-dependent boundaries, which were used to

represent the Rio Grande, irrigation canals, and

agricultural drains. During model calibration,

parameter values were adjusted to minimize

differences between computed and measured heads and

fluxes. A total of 4,439 measurements of heads and

flow losses or gains compose the calibration data set

(fig. 12).

Hydraulic-Head Measurements

The nonlinear regression included 4,352 heads

measured between 1912 and 1995 in 310 wells. These

measurements were compiled from databases of

EPWU and the U.S. Geological Survey (USGS). About

83 percent (3,615) of the measurements were made in

226 wells in the United States, and the remainder (737)

were made in 84 wells in Mexico. Most of the heads

have been measured since ground-water withdrawals

increased significantly in the 1960’s (fig. 12). Of the

total set of measurements, 4,184 were made in wells

with screened intervals corresponding to more than one

model layer; these were defined as multilayer head

measurements in MODFLOWP.

Flow Measurements

Measurements of streamflow loss in the Rio

Grande and Franklin Canal and seepage into the Island

Drain were used to constrain the nonlinear regression.

The IBWC computed seepage loss from the Rio

Grande into the underlying aquifer along a 15-km (9.3-

mi) reach between the concrete-lined section in the

Chamizal zone and Riverside Dam for 1981, 1982, and

1983 (Land and Armstrong, 1985; White and others,

1997). The BOR computed flow loss along an 8.5-km

(5.25-mi) reach of the Franklin Canal for 1984, and the

USGS computed flow loss for 1990 through 1992

(Land and Armstrong, 1985; White and others, 1997).

These computations and simulated flow losses are

summarized in table 1.

Although flow measured in some agricultural

drains may include tributary inflow, such inflow to the

Island Drain is minor. Twenty-four measurements of

flow in the Island Drain were used to constrain drainbed hydraulic conductance in the regression.

Parameter Estimation

Parameter values representing aquifer properties

and specified boundaries were estimated through

nonlinear regression that minimized differences

between measured and computed heads and flows in a

94-year transient-state simulation from 1903 through

1996. Hydraulic heads computed in a steady-state
24

Figure 11. Locations of head measurements, stream segments, agricultural drains, and model faults.

141

59

85

14 15

30

43

202

209

212

176

123

207

405

413

416

148

0 2 4 6 8 10 MILES

0 2 4 6 8 10 KILOMETERS

106 ̊30 ́ 106 ̊22 ́30 ́ ́ 106 ̊15 ́

32 ̊

31 ̊

52 ́

30 ́ ́

31 ̊

45 ́

Location of wells in figures 6 and 7

Location of wells with hydrographs in figure 17

Active model-grid boundary

Agricultural drains

Faults as represented in flow model

EXPLANATION

104

86

124 Franklin 

Canal

Grande

Rio

Madre

Acequia

Shaded-relief map from U.S. Geological Survey

30-meter Digital Elevation Models

Location of other water-level measurement sites
25

simulation representing prepumping conditions, when

no large stresses were imposed on the aquifer system,

provide initial conditions for the transient-state

simulation. The nonlinear regression procedure was

implemented using UCODE, in which sensitivities of

model parameters are estimated through a

perturbation technique (Poeter and Hill, 1998).

Although MODFLOWP provides a similar regression

procedure in which sensitivities of model parameters

are computed directly from an analytical expression,

a major limitation prevented its application for model

calibration. Conversion from confined to unconfined

conditions caused by dewatering the upper part of the

modeled area can be represented in MODFLOWP in

the solution for hydraulic head, but this conversion was

not implemented in MODFLOWP for the regression

procedure. The calibration procedure selected to

circumvent this problem entailed (1) using

MODFLOWP to define the parameter values in the

model and compute hydraulic heads and flows and (2)

running UCODE to estimate the optimum parameter

values.

Improvements in model results were identified

by comparing the sum of squared errors (SSE):

, i = 1, n (9)

where ei = difference between measured and

calculated values of

measurement i;

wi

1/2 = square root of the weight

assigned to the error in the

measured value of

measurement i;

wi

1/2ei = weighted residual corresponding

to measurement i; and

n = number of measurements;

and the standard error of

estimate (SEE):

, (10)

where p = number of model parameters estimated by

the regression.

SSE Σ wi

1 2⁄

ei [ ]

2 =

SEE SSE

n p – ------------ 1 2⁄

=

Figure 12. Number of annual head measurements included in the nonlinear regression used in the

ground-water flow model.

in wells

1903

1908

1913

1918

1923

19281933

1938

1943

1948

1953

1958

1963

1968

1973

1978

1983

1988

1993

100

120

0

20

40

60

80

140

160

Number of measurements

United States

Mexico
26

The weights, wi, were chosen according to

procedures in Hill (1992) to account for the different

units associated with head measurements (m) and flow

measurements (m3

/d). The water-level altitudes

reported for wells in the United States used a datum

consistent with other data in the flow model. Heads in

the United States were generally measured several days

after the cessation of pumping in wells to better

represent nonpumping conditions. The time elapsed

after cessation of pumping is not known for wells in

Mexico, however, and the reported water-level altitudes

are less certain. The wi

 values were therefore adjusted

so that heads in wells in Mexico were weighted 33

percent less than heads in wells in the United States.

Heads measured in the United States were weighted

equally. Flow measurements were weighted to reflect

the relative accuracy of Rio Grande and canal seepageloss measurements and drain-flow measurements.

Measurements of Rio Grande seepage loss (table 1)

were weighted equally to heads measured in wells in

the United States. Measurements of Franklin Canal

seepage loss and flow in the Island Drain were

weighted 67 percent less than measurements of Rio

Grande seepage loss.

MODEL EVALUATION AND SIMULATION

RESULTS

The spatial distribution of horizontal and vertical

hydraulic conductivities was defined by assigning

zones to each of the hydrogeologic facies, to which

corresponding parameters representing hydraulic

conductivities were applied. The horizontal extent of

the recent alluvial facies (fig. 3) is known because

deposition was constrained by the topography of the

Rio Grande Valley. Boundaries between the fluvial,

lacustrine-playa, and alluvial-fan facies are less

constrained. Two possible geometries for the boundary

between the fluvial and lacustrine-playa facies were

tested: (1) a “block” model with a vertical facies

boundary and (2) a “wedge” model in which the facies

boundary sloped down to the west so that lacustrineplaya facies underlie the fluvial facies in parts of layers

2 through 9. In both models, the lacustrine-playa facies

accommodated all of layer 10 and the hydraulic

conductivities between the two facies were gradational,

simulating some facies interfingering. For each model,

the parameter-estimation regression was run to

determine the lowest SSE parameter set. The resulting

parameter sets for these two models differed only

slightly. The “block” model was found to have lower

overall error; the geometry of this hydraulicconductivity distribution is depicted in cut-away

perspective in figure 4. Although the horizontal

hydraulic conductivities of the alluvial-fan and fluvial

facies were independently estimated, they converged to

essentially identical values. This suggests that the

location of this facies boundary may not be significant

in the context of the ground-water flow model.

The total SSE computed using equation 7 was

38,814 m2. The SSE for all model heads computed with

equation 8 was 2.99 m2

. The SSE for heads in United

States wells was 2.90 m2 and for heads in Mexican

wells was 3.44 m2

. Scatter plots of simulated heads in

relation to measured heads and of simulated heads in

relation to weighted residuals for wells in Mexico and

the United States are shown in figure 13C-D. The

scatter plots for the United States wells show no

obvious evidence of model bias. The scatter plot of

simulated heads in relation to weighted residuals for

wells in Mexico suggests some bias toward negative

weighted residuals. This bias indicates that heads in

Mexican wells are, in general, slightly lower than those

predicted by the ground-water flow model.

Estimates of Aquifer Properties

The best-fit parameter values along with their

95-percent confidence intervals computed from the

nonlinear least-squares regression of the “block” facies

model are summarized in table 3. Because the lengths

of the Rio Grande, Franklin Canal, Acequia Madre, and

agricultural drains in each model cell are well known

but the thickness of bed material is uncertain, these

parameters are reported as a conductance per unit

length. Quaternary fault-zone hydraulic conductivity

was estimated on the basis of a fault-zone thickness of

1 m. If actual fault-zone thickness is on the order of 10

or 1 cm, the actual fault-zone hydraulic conductivity is

one or two orders of magnitude smaller, respectively.

Parameter Sensitivities

Composite scaled sensitivities (Hill, 1998) are

listed in table 2 for model parameters estimated in the

nonlinear regression. These sensitivities vary

somewhat as a function of parameter values and

distribution. The sensitivities are useful for interpreting
27

Figure 13. Simulated and measured heads and weighted residuals computed in transient-state simulation

of the Hueco Bolson.

United States wells

1,070 1,080 1,090 1,100 1,110 1,120 1,130 1,140

Simulated head, in meters above sea level

1,070

1,080

1,090

1,100

1,110

1,120

1,130

1,140

Measured head, in meters

above sea level

B

-20

-15

-10

-5

0

5

10

15

20

1,070 1,080 1,090 1,100 1,110 1,120 1,130 1,140

Simulated head, in meters above sea level

Mexico wells

Residual head (measured minus simulated)

C

-20

-15

-10

-5

0

5

10

15

20

1,070 1,080 1,090 1,100 1,110 1,120 1,130 1,140

Simulated head, in meters above sea level

United States wells

Residual head (measured minus simulated)

D

Mexico wells

1,070 1,080 1,090 1,100 1,110 1,120 1,130 1,140

Simulated head, in meters above sea level

1,070

1,080

1,090

1,100

1,110

1,120

1,130

1,140

Measured head, in meters

above sea level

A
28

1Maximum for normal irrigation year after 1933.

2Distributed in arroyos below Organ, Franklin, and Juarez Mountains.

Table 3. Optimum parameter values estimated for Hueco Bolson aquifer, through nonlinear regression

in transient-state simulation, and their approximate 95-percent confidence intervals

[m, meters; m/d, meters per day; m3/d, cubic meters per day; N/A, not applicable]

Parameter Value Confidence interval

Scaled sensitivity,

in meters

Recharge:

Irrigation-return flow1 3.8 x 104 m3/d 2.2 x 104 - 6.5 x 104 m3/d 0.3

Underflow from Tularosa Basin 2.0 x 104

 m3/d 1.9 x 104 - 2.2 x 104 m3/d 7

Mountain front on alluvial fans2 8.0 x 102 m3/d 0 - 2 x 103 m3/d 2

Underflow from Mesilla Basin 3.4 x 102

 m3/d Not estimated N/A

Horizontal hydraulic conductivity:

Alluvial-fan facies 6.8 m/d 6.0 - 7.7 m/d 5

Recent fluvial sediments 4.0 m/d 2.8 - 7.2 m/d 2

Fluvial and alluvial facies 6.8 m/d 6.4 - 7.2 m/d 10

Lacustrine-playa facies 0.9 m/d 0.5 - 1.4 m/d 0.3

Quaternary faults 3.8 x 10-3 m/d 1 x 10-3 - 1 x 10-2 m/d 8

Vertical hydraulic conductivity:

Alluvial-fan facies 1.2 x 10-3 m/d 6 x 10-4 - 2 x 10-3 m/d 3

Recent fluvial sediments 1.3 x 10-1 m/d 6 x 10-2 - 6 x 10-1 m/d 2

Fluvial and alluvial facies 1.3 x 10-2 m/d 1 x 10-2 - 1.5 x 10-2 m/d 0.6

Lacustrine-playa facies 2.5 x 10-2 m/d 6 x 10-3 - 1 x 10-1 m/d 0.4

Specific yield 0.178 0.173 - 0.184 10

Specific storage - elastic 7 x 10-6 m-1 2 x 10-6 - 1 x 10-5 m-1 0.5

Specific storage - inelastic 7 x 10-5 m-1 Not estimated N/A

Conductance per unit length:

Rio Grande 1.8 m/d 1.6 - 2.0 m/d 1

Agricultural drains 5.8 m/d 2 - 1.6 x 101

 m/d 2

Irrigation canals 3.0 m/d Not estimated N/A

Evapotranspiration extinction depth 5 m Not estimated N/A

Maximum evapotranspiration rate 4.6 x 10-3 m/d 1 x 10-3 - 7 x 10-3 m/d 2

Manning’s n:

Rio Grande 0.03 Not estimated N/A

Franklin Canal, Acequia Madre 0.004 Not estimated N/A
29

relative sensitivities between model parameters and the

relative importance of a parameter in reducing model

error in the nonlinear regression. Specific yield of

unconfined layers and horizontal hydraulic

conductivity of the fluvial facies were the most

sensitive model parameters, followed by recharge

underflow from the Tularosa Basin and conductance of

Quaternary fault zones.

Water-Level Drawdowns

Simulated water-table altitudes for 1902 (steadystate conditions) and for 1958, 1973, and 1980

(transient-state conditions) are contoured at 2-m

intervals and shown in figure 14A-D. The progressive

growth of cones of depression under El Paso and

Ciudad Juarez is evident from 1958 through 1980.

Shallow ground water flows away from the Rio Grande

toward these cones of depression on either side of the

international border.

Potentiometric surface altitudes computed for

model layer 5, which is approximately 135 m (440 ft)

below the steady-state water table, for 1958, 1973,

1980, and 1996 are shown in figure 15A-D. This model

layer is in the middle of the depth interval in which

many production wells are screened in the Hueco

Bolson. The contours for 1958 through 1996 show the

expanding cones of depression under El Paso and

Ciudad Juarez; the Rio Grande, however, has much less

influence at this depth than at the water table. Since

about 1980 (fig. 15C), water at this depth generally

flows from north to south beneath the river toward the

lower potentiometric heads beneath Ciudad Juarez. By

1973, ground-water drawdown under Ciudad Juarez

had created a ground-water divide to the southeast,

which isolated ground-water flow near the city from the

regional pattern of southeast flow. This divide has

migrated farther to the southeast as the cone of

depression has deepened. The influence of simulated

Quaternary faults (fig. 11) becomes increasingly

evident as deflections of water-level contours in later

transient stress periods, such as depicted for 1980 and

1996 (fig. 15C,D).

The three-dimensional pattern of ground-water

drawdown beneath El Paso and Ciudad Juarez is

further illustrated by perspective views for 1973 and

1996 (fig. 16A,B). The vertical section of these images

was sliced along an azimuth of N. 30° E. to transect the

drawdown cones beneath both cities. Only active

model cells are shown; missing cells in the El Paso and

Ciudad Juarez regions in 1996 (fig. 16B) result from

deactivation by ground-water drawdown beneath the

bottom of model layer 1. Model layer 2 has similarly

deactivated in some areas.

Hydrographs of measured and model-simulated

water levels for 18 wells are presented in figures

17A-G. The locations of these wells, which are

representative of the other 274 wells in the

measurement set, are shown in figure 11. These

hydrographs illustrate the general quality of model fit

and may be useful for showing areas where the model

could be improved. (A spreadsheet containing

hydrograph measurements for the remaining 274 wells

is available on request.)

Ground-Water Budget

Components of the ground-water budget that

change during the course of the transient simulation are

shown in figure 18. Water pumped from wells is

supplied primarily by releases from ground-water

storage, principally by drainage of aquifer pore space;

this process results in large declines of the water table.

The mirror relation between ground-water pumping

and releases from storage is evident in figure 18A,B. A

similar mirrorlike relation can be seen in figure 18C,D,

which illustrates the principal budget components in

the shallow portion of the aquifer system in the Rio

Grande Valley. These two graphs indicate that Rio

Grande water infiltrating into the shallow aquifer

system is consumed primarily by ET and (or) flows to

agricultural drains.

Seepage from the Rio Grande

The simulated seepage losses from the Rio

Grande between the end of the Chamizal zone and

Riverside Dam for two historical and one hypothetical

(design flow) 5-year periods are depicted in figure

19A-D. Although model fluxes and other report figures

have units in meters per day, these seepage losses are

presented in acre-ft/yr for the benefit of interested

parties. Variations in seepage by month from 1975 to

1979 and 1988 to 1992 are shown in figure 19A.

Annual seepage loss was divided into seepage during

the primary and secondary irrigation seasons in figure

19B-D. Seepage losses for 1975-79 and 1988-92 were

during times of relatively low and high flow in the Rio

Grande, respectively, and correspond to the monthly

periods illustrated in figure 10. The greater seepage in

1988-92 principally results from higher average stage

of the Rio Grande.
30

Figure 14A. Simulated water table in shallow aquifer (model layers 1 and 2) in 1902

(steady-state conditions).

Shaded-relief map from

U.S. Geological Survey

30-meter Digital

Elevation Models

106 ̊30 ́ 106 ̊15 ́ 106 ̊00 ́

32 ̊

15 ́

32 ̊

00 ́

31 ̊

45 ́

31 ̊

30 ́

0 5 10 15 20 25 30 35 40 MILES

0 5 10 15 20 25 30 35 40 KILOMETERS

Ciudad Juarez

El Paso

1082

1084

1086

1088

1090

1092

1094

1096

1098

1100

1102

1104

1106

1108

1110

1112

1114

1116

1118

1120

1122

1124

1126

1128

1132

1134

1136

1138

1140

1142

1144

1146

1148

1150

1152

1154

1156

EXPLANATION

Water-table contour--Shows

altitude of water table.

Contour interval 2 meters.

Datum is sea level
31

106 ̊30 ́ 106 ̊15 ́ 106 ̊00 ́

32 ̊

15 ́

32 ̊

00 ́

31 ̊

45 ́

31 ̊

30 ́

Figure 14B. Simulated water table in shallow aquifer (model layers 1 and 2) in 1958.

0 5 10 15 20 25 30 35 40 MILES

0 5 10 15 20 25 30 35 40 KILOMETERS

1082

1084

1086

1088

1090

1092

1094

1096

1098

1100

1102

1104

11

1106

08

1110

1112

1114

1116

1116

1118

1118

1120

1120

1122

11261124

1128

1130

1132

1134

1136

1138

1142

1144

1146

1148

1150

1152

1154

1156

Ciudad Juarez

El Paso

Shaded-relief map from

U.S. Geological Survey

30-meter Digital

Elevation Models

EXPLANATION

Water-table contour--Shows

altitude of water table.

Contour interval 2 meters.

Datum is sea level
32

106 ̊30 ́ 106 ̊15 ́ 106 ̊00 ́

32 ̊

15 ́

32 ̊

00 ́

31 ̊

45 ́

31 ̊

30 ́

0 5 10 15 20 25 30 35 40 MILES

0 5 10 15 20 25 30 35 40 KILOMETERS

1082

1084

1086

1088

1090

1092

1094

1096

1098

1100

1102

1104

1106

1106

1108

1108

1110

1110

1112

1112

1114

1114

1114

1116

1116

1118 1120

1122

112

1126 4

1128

1130

1132

1134

1136

1138

1140

1144

1148

1150

1152

1154

1156

1158

Ciudad Juarez

El Paso

Figure 14C. Simulated water table in shallow aquifer (model layers 1 and 2) in 1973.

Shaded-relief map from

U.S. Geological Survey

30-meter Digital

Elevation Models

EXPLANATION

Water-table contour--Shows

altitude of water table.

Contour interval 2 meters.

Datum is sea level
33

Figure 14D. Simulated water table in shallow aquifer (model layers 1 and 2) in 1980.

106 ̊30 ́ 106 ̊15 ́ 106 ̊00 ́

32 ̊

15 ́

32 ̊

00 ́

31 ̊

45 ́

31 ̊

30 ́

0 5 10 15 20 25 30 35 40 MILES

0 5 10 15 20 25 30 35 40 KILOMETERS

1082

1084

1086

1088

1090

1092

1094

1096

1098

1098

1100

1100

1102

1102

1104

1104

1106

1106

1108

1108

1110

11

1110

12

1112

1114

1114

1116

1116

1118

1120

1122

1124

1126

1128

1130

1132

1134

1136

1138

1140

1142

1144

1146

1148

1150

1152

1154

1156

Ciudad Juarez

El Paso

Shaded-relief map from

U.S. Geological Survey

30-meter Digital

Elevation Models

EXPLANATION

Water-table contour--Shows

altitude of water table.

Contour interval 2 meters.

Datum is sea level
34

Figure 15A. Simulated potentiometric surface in model layer 5 in 1958.

106 ̊30 ́ 106 ̊15 ́ 106 ̊00 ́

32 ̊

15 ́

32 ̊

00 ́

31 ̊

45 ́

31 ̊

30 ́

0 5 10 15 20 25 30 35 40 MILES

0 5 10 15 20 25 30 35 40 KILOMETERS

1082

1084

1086

1090

1088

1092

1094

1096

1098

110

1102

0

1104

1106

1108

1110

1112

1114 1114

1116

1116

1118

1120

1122

1124

1126

1128

1130

1132

1134

1136

1138

1142

1144

1146

1148

1150

1152

1154

Shaded-relief map from

U.S. Geological Survey

30-meter Digital

Elevation Models

Ciudad Juarez

El Paso

EXPLANATION

Potentiometric contour–-

Shows altitude at which

water level would have

stood in tightly cased

wells. Contour interval

2 meters. Datum is sea level
35

Figure 15B. Simulated potentiometric surface in model layer 5 in 1973.

106 ̊30 ́ 106 ̊15 ́ 106 ̊00 ́

32 ̊

15 ́

32 ̊

00 ́

31 ̊

45 ́

31 ̊

30 ́

0 5 10 15 20 25 30 35 40 MILES

0 5 10 15 20 25 30 35 40 KILOMETERS

B'

A'

A

1104

1106

1108 1110

1112

1082

1084

1086

1088

1090

1092

1094

1096

1098

1100

1102

1104

1106

1108

1110

1112

1114

1116

1118

1120

1122

1124

1126

1128

1130

1132

1134

1136

1138

1140

1142

1144

1146

1148

1150

1152

1154

Shaded-relief map from

U.S. Geological Survey

30-meter Digital

Elevation Models

Ciudad Juarez

EXPLANATION

Potentiometric contour–-

Shows altitude at which

water level would have

stood in tightly cased

wells. Contour interval

El Paso 2 meters. Datum is sea level
36

Figure 15C. Simulated potentiometric surface in model layer 5 in 1980.

106 ̊30 ́ 106 ̊15 ́ 106 ̊00 ́

32 ̊

15 ́

32 ̊

00 ́

31 ̊

45 ́

31 ̊

30 ́

0 5 10 15 20 25 30 35 40 MILES

0 5 10 15 20 25 30 35 40 KILOMETERS

1098

1098

1100

 1102 1104 1106

1108

1110

1112

1082

1084

1086

1090

1088

1092

1094

1096

1098

1100

1102

1104

1106

1108

1110

1112

1114

1116

1118

1120

1122

1124

1126

1128

1130

1132

1134

1136

1138

1140

1142

1144

1146

1148

1150

1152

1154

Shaded-relief map from

U.S. Geological Survey

30-meter Digital

Elevation Models

Ciudad Juarez

El Paso

EXPLANATION

Potentiometric contour–-

Shows altitude at which

water level would have

stood in tightly cased

wells. Contour interval

2 meters. Datum is sea level
37

EXPLANATION

Potentiometric contour–-

Shows altitude at which

water level would have

stood in tightly cased

wells. Contour interval

2 meters. Datum is sea level El Paso

Juarez

Figure 15D. Simulated potentiometric surface in model layer 5 in 1996.

106 ̊30 ́ 106 ̊15 ́ 106 ̊00 ́

32 ̊

15 ́

32 ̊

00 ́

31 ̊

45 ́

31 ̊

30 ́

0 5 10 15 20 25 30 35 40 MILES

0 5 10 15 20 25 30 35 40 KILOMETERS

1076

10781080 1082

1084

1086 1088

1090

1092 1094

1096 1098

1100 1102

1104

1106

1108

1084

1086

1090

1088

1092

1094

1096

1098

1100

1102

1104

1106

1110

1112

1114

1116

1118

1120

1122

1124

1126

1128

1130

1132

1134

1136

1138

1140

1142

1144

1146

1148

1150

1152

Shaded-relief map from

U.S. Geological Survey

30-meter Digital

Elevation Models

Ciudad Juarez

El Paso
38

604530161.0

Figure 16A. Perspective view showing water-level drawdown in the Hueco Bolson aquifer in 1973.

Logarithmic color scale.

Water-leveldrawdown,in meters

CIUDAD JUAREZ

EL PASO
39

60

45

30

16

1.0

Figure 16B. Perspective view showing water-level drawdown in the Hueco Bolson aquifer in 1996.

Logarithmic color scale.

Water-leveldrawdown,in meters

EL PASO

CIUDAD JUAREZ
40

W a t e r - l e v e l a l t i t u d e, i n m e t e r s a b o v e s e a l e v e l

JL-49-05-602 (well 14)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

Measured Simulated

Measured Simulated

JL-49-06-401 (well 15)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

1930 1940 1950 1960 1970 1980 1990

Figure 17A. Measured and simulated water levels from 1930 through 1995 in selected wells in the

Mesa-Nevins well field (location of wells shown in figure 11).
41

Figure 17A. Measured and simulated water levels from 1930 through 1995 in selected wells in the

Mesa-Nevins well field (location of wells shown in figure 11)--Concluded.

JL-49-05-901 (well 30)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

JL-49-13-204 (well 43)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

1930 1940 1950 1960 1970 1980 1990

Measured Simulated

Measured Simulated

W a t e r - l e v e l a l t i t u d e, i n m e t e r s a b o v e s e a l e v e l
42

JL-49-13-601 (well 59)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

JL-49-13-610 (well 86)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

1930 1940 1950 1960 1970 1980 1990

W a t e r - l e v e l a l t i t u d e, i n m e t e r s a b o v e s e a l e v e l

Figure 17B. Measured and simulated water levels from 1930 through 1995 in selected wells in the

Airport well field (location of wells shown in figure 11).

Measured Simulated

Measured Simulated
43

Figure 17C. Measured and simulated water levels from 1930 through 1995 in selected wells in the

Cielo Vista and Town well fields (location of wells shown in figure 11).

1,075

1,085

1,095

1,105

1,115

1,125

1,135

JL-49-14-413 (well 85)

JL-49-13-804 (well 104)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

1930 1940 1950 1960 1970 1980 1990

W a t e r - l e v e l a l t i t u d e, i n m e t e r s a b o v e s e a l e v e l

Measured Simulated

Measured Simulated
44

JL-49-13-702 (well 124)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

JL-49-13-710 (well 141)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

1930 1940 1950 1960 1970 1980 1990

W a t e r - l e v e l a l t i t u d e, i n m e t e r s a b o v e s e a l e v e l

Figure 17D. Measured and simulated water levels from 1930 through 1995 in selected wells in the

water-plant well field (location of wells shown in figure 11).

Measured Simulated

Measured Simulated
45

JMAS 10-R (well 176)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

JMAS 17 (well 202)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

1930 1940 1950 1960 1970 1980 1990

W a t e r - l e v e l a l t i t u d e, i n m e t e r s a b o v e s e a l e v e l

Figure 17E. Measured and simulated water levels from 1930 through 1995 in selected wells in the

Juarez well field (location of wells shown in figure 11).

Measured Simulated

Measured Simulated
46

JMAS13-R (well 209)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

JMAS42 (well 212)

1930 1940 1950 1960 1970 1980 1990

1,075

1,085

1,095

1,105

1,115

1,125

1,135

W a t e r - l e v e l a l t i t u d e, i n m e t e r s a b o v e s e a l e v e l Figure 17E. Measured and simulated water levels from 1930 through 1995 in selected wells in the

Juarez well field (location of wells shown in figure 11)--Concluded.

Measured Simulated

Measured Simulated
47

JL-49-14-102 (well 413)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

JL-49-14-301 (well 416)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

1930 1940 1950 1960 1970 1980 1990

W a t e r - l e v e l a l t i t u d e, i n m e t e r s a b o v e s e a l e v e l

Figure 17F. Measured and simulated water levels from 1930 through 1995 in selected wells in the

east observation well field (location of wells shown in figure 11).

Measured Simulated

Measured Simulated
48

JL-49-21-301 (well 148)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

JL-49-05-502 (well 405)

1,075

1,085

1,095

1,105

1,115

1,125

1,135

1930 1940 1950 1960 1970 1980 1990

W a t e r - l e v e l a l t i t u d e, i n m e t e r s a b o v e s e a l e v e l

Figure 17G. Measured and simulated water levels from 1930 through 1995 in selected wells in the lower

valley and north observation well fields (location of wells shown in figure 11).

Measured Simulated

Measured Simulated
49

T h o u s a n d s o f c u b i c m e t e r s p e r d a y

(A) Ground-water withdrawals

-800

-700

-600

-500

-400

-300

-200

-100

0

(C) Stream leakage

0

100

200

300

400

500

600

700

800

(B) Storage

-200

-100

0

100

200

300

400

500

600

-800

-700

-600

-500

-400

-300

-200

-100

0

(D)

1900

19051910

191519201925

193019351940

19451950

195519601965

19701975

198019851990

1995

Figure 18. Simulated annual inflows and outflows (-) to Hueco Bolson aquifer.

Agricultural drains

Evapotranspiration
50

T h o u s a n d s o f a c r e - f e e t p e r y e a r 

(A) Monthly seepage

0

10

20

30

40

50

1975-79 1988-92

Jan Feb Mar Apr May June July Aug Sept Oct Nov Dec

(B) 1975-79

0

5

10

15

20

25

30

35

40

45

Jan 1975 Jan 1976 Jan 1977 Jan 1978 Jan 1979

(C) 1988-92

0

5

10

15

20

25

30

35

40

45

Jan 1988 Jan 1989 Jan 1990 Jan 1991 Jan 1992

Figure 19. Rio Grande seepage between Chamizal zone and Riverside Dam (locations shown

in figure 8).

(D) Rio Grande seepage with design flow

0

5

10

15

20

25

30

35

40

45

Total annual Primary irrigation season Secondary irrigation season

Jan

year 1

Jan

year 2

Jan

year 3

Jan

year 4

Jan

year 5

EXPLANATION
51

The simulated seepage for the 10th year of a

“design flow” in the Rio Grande without the ACE is

illustrated in figure 10D. This design-flow seepage-loss

calculation is for a hypothetical scenario in which Rio

Grande water is not diverted into the ACE. The

simulated aquifer conditions for this scenario

approximate those that would have occurred during

2002 if the monthly ground-water pumpage rates and

patterns in 1992 repeated for 10 years. The computed

seepage loss from the Rio Grande to the underlying

aquifer along the 15-km reach between the end of the

Chamizal zone and Riverside Dam is 39,000

acre-ft/yr in this scenario.

The altitudes of drain beds specified in the model

influence simulated water-table altitudes in portions of

the Rio Grande Valley, which could affect seepage rates

from the Rio Grande to the underlying aquifer.

Although uncertainty exists for the drain-bed altitudes

specified in the simulation (which are probably +1, -2

m), these altitudes were not formally estimated during

the model-calibration process. To test the possible Rio

Grande seepage-loss effect of the specified drain-bed

altitudes, several forward model runs were made with

drain-bed altitudes decreased in increments of 1 m.

Although drain seepage increased, seepage loss from

the Rio Grande between the end of the Chamizal zone

and Riverside Dam remained at 39,000 acre-ft/yr.

Within the tested drain-bed altitude range, which

brackets values deemed reasonable, seepage loss from

this reach of the Rio Grande is insensitive to specified

drain-bed altitude.

The calibration on both the computed Rio

Grande seepage losses and the adjustments to

measured flow required for stream-package flow

specification (see “Streamflow routing” section) was

further checked by comparing evaporation-corrected,

model-simulated flow at Riverside Dam with measured

flow (computed as the sum of measured flow in the

Riverside Canal and measured flow at Coffer Dam).

These computations are shown in figure 20 for

1988-92. In general, the model fit is excellent, with the

exception of 1990. Of the 348 monthly stress periods

not shown in this figure, 346 had an excellent fit

comparable with those of 1988-89 and 1991-92. The 2

years with months of substandard fit were 1975 and

1986. The years of substandard fit may have resulted

from Rio Grande flow-measurement errors.

SUMMARY AND CONCLUSIONS

The neighboring cities of El Paso, Texas, and

Ciudad Juarez, Chihuahua, Mexico, have historically

relied on ground-water withdrawals from the Hueco

Bolson, an alluvial-aquifer system, to supply water to

their growing populations. In the United States,

diversions from the Rio Grande and ground-water

withdrawals from the Mesilla Basin west of the study

area also supply the freshwater demands of the

military, industries, and public in the El Paso area. In

Mexico, diversions from the Rio Grande are used for

agriculture; water needed by Ciudad Juarez is supplied

solely by extraction from the Hueco Bolson. By 1996,

ground-water drawdown exceeded 60 m in some areas

under Ciudad Juarez and El Paso. Fresh ground water

stored in the aquifer system beneath these cities is

bordered by regions of brackish to saline ground water.

As water levels in the freshwater portions of the aquifer

declined, intrusion of the surrounding brackish water

degraded water quality in public supply wells, which

sometimes required well abandonment.

By the end of 1996, ground-water levels had

declined by as much as 60 m (197 ft) from simulated

steady-state levels in the Hueco Bolson. By

incorporating extensive pumpage records with fine

spatial and temporal discretization, these declines have

been simulated in a numerical ground-water flow

model that matches 4,352 head observations with an

SEE of about 3 m. This monthly temporal

discretization was needed to (1) permit accurate

calculations of seepage losses from the Rio Grande and

(2) provide a useful ground-water management model.

The simulation of steady-state and transient

ground-water flow in the Hueco Bolson was developed

using MODFLOW-96. The transient simulation

represents a period of 100 years beginning in 1903 and

ending in 2002. The period 1903 through 1968 was

represented with 66 annual stress periods, and the

period 1969 through 2002 was represented with 408

monthly stress periods. Model boundary conditions

were modified at appropriate times during the

simulation to represent changes in well pumpage,

drainage of agricultural fields, and channel

modifications of the Rio Grande.

The model was calibrated using MODFLOWP

and UCODE. Parameter values representing aquifer

properties and boundary conditions were adjusted

through nonlinear regression in a transient-state

simulation with 96 annual time steps to produce a

model that approximated (1) 4,352 water levels
52

Figure 20. Simulated and measured flows at Riverside Dam (location shown in figure 8).

0

500

1,000

1,500

2,000

2,500

3,000

3,500

4,000

Measured Simulated

Jan

Mar

May

July

Sept

Nov

Jan

Mar

May

July

Sept

Nov

Jan

Mar

May

July

Sept

Nov

Jan

Mar

May

July

Sept

Nov

Jan

Mar

May

July

Sept

Nov

Flow, in thousands of cubic meters per day

1988 1989 1990 1991 1992
53

measured in 292 wells from 1912 to 1995, (2) three

seepage-loss rates from a reach of the Rio Grande

during periods from 1979 to 1981, (3) three seepageloss rates from a reach of the Franklin Canal during

periods from 1990 to 1992, and (4) 24 seepage rates

into irrigation drains from 1961 to 1983. Once a

calibrated model was obtained with MODFLOWP and

UCODE, the optimal parameter set was used to create

an equivalent MODFLOW-96 simulation with monthly

temporal discretization to improve computations of

seepage from the Rio Grande and to define the flow

field for a chloride-transport simulation.

The optimal values for a set of 17 parameters

were obtained using nonlinear regression. Values of

other parameters that were either well constrained or

insensitive to the model were specified. The regression

was constrained by the head and flow-loss

measurements from the Rio Grande and Franklin Canal

and flow observations in agricultural drains. The model

was most sensitive to (1) horizontal hydraulic

conductivity of the fluvial facies, which composes the

principal aquifer material in which production wells

are drilled, (2) specific yield, (3) recharge underflow

from the Tularosa Basin, and (4) hydraulic

conductance of Quaternary fault zones. Parameter

distribution was generated from a simplified

hydrogeologic model. Several more complex geologic

conceptualizations, such as a dipping boundary

between fluvial and lacustrine-playa facies, were

parameterized and optimized as regression runs. These

alternative models, though probably more geologically

realistic, did not provide an overall improved fit to the

hydraulic-head measurements. Future geologic

refinements incorporated into the flow model may

improve model fit in certain areas, however. Seepage

losses for the 15-km reach of the Rio Grande channel

between the Chamizal zone and Riverside Dam were

computed with monthly stress periods. The seepage

loss from the Rio Grande for a hypothetical "design

flow” in the Rio Grande was 39,000 acre-ft/yr and was

not sensitive to specified drain-bed altitudes within a

reasonable range.

The model input was generated from GIS

databases, which facilitated rapid model construction

and enabled testing of several conceptualizations of

hydrogeologic facies boundaries. The simulation

results were sensitive to the hydraulic conductance of

Quaternary faults in the fluvial-aquifer facies,

suggesting that ground-water flow is impeded across

the fault planes.

SELECTED REFERENCES

Anderholm, S.K., and Heywood, C.E., 2003, Chemistry and

age of ground water in southwestern Hueco Basin, New

Mexico and Texas: U.S. Geological Survey WaterResources Investigations Report 02-4237.

Barnes, V.E., 1983, Geologic atlas of Texas, Van Horn - El

Paso Sheet: Austin, University of Texas, Bureau of

Economic Geology, scale 1:250,000.

Bennett, G.D., Kontis, A.L., and Larson, S.P., 1982,

Representation of multiaquifer well effects in threedimensional ground-water flow simulation: Ground

Water, v. 20, no. 3, p. 334-341.

Chapin, C.E., and Seager, W.R., 1975, Evolution of the Rio

Grande rift in the Socorro and Las Cruces areas: New

Mexico Geological Society, 26th Field Conference,

p. 297-321.

Culler, R.C., Hanson, R.L., Myrick, R.M., Turner, R.M., and

Kipple, F.P., 1982, Evapotranspiration before and after

clearing phreatophytes, Gila River flood plain, Graham

County, Arizona, New Mexico: U.S. Geological Survey

Professional Paper 655-P, 67 p.

Groschen, G.E., 1994, Simulation of ground-water flow and

the movement of saline water in the Hueco Bolson

aquifer, El Paso, Texas, and adjacent areas: U.S.

Geological Survey Open-File Report 92-171, 87 p.

Haneberg, W.C., Mozley, P.S., Moore, J.C., and Goodwin,

L.B., 1999, Faults and subsurface fluid flow in the

shallow crust: Geophysical Monograph 113, American

Geophysical Union, Washington, D.C., 222 p.

Harbaugh, A.W., and McDonald, M.G., 1996, User’s

documentation for MODFLOW-96, an update to the

U.S. Geological Survey modular finite-difference

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

Open-File Report 96-485, 56 p.

Heywood, C.E., 1993, Monitoring aquifer compaction and

land subsidence due to ground-water withdrawal in the

El Paso, Texas - Juarez, Chihuahua area, in Prince,

K.R., Galloway, D.L., and Leake, S.A., eds.: U.S.

Geological Survey Subsidence Interest Group

Conference, Edwards Air Force Base, Antelope Valley,

Calif., November 18-19, 1992—Abstracts and

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

94-532, 84 p.

———1995, Investigation of aquifer-system compaction in

the Hueco Basin, El Paso, Texas, USA, in Proceedings

of the Fifth International Symposium on Land

Subsidence: International Association of Hydrological

Sciences Publication 234, p. 35-45.

Hill, M.C., 1992, A computer program (MODFLOWP) for

estimating parameters of a transient, three-dimensional

ground-water flow model using nonlinear regression:

U.S. Geological Survey Open-File Report 91-484,

382 p.
54

Hill, M.C., 1998, Methods and guidelines for effective

model calibration: U.S. Geological Survey WaterResources Investigations Report 98-4005, 90 p.

Holzbecher, E.H., 1998, Modeling density driven flow in

porous media: Berlin, Springer Verlag, 286 p.

Hsieh, P., 1992, Documentation of a computer program to

simulate horizontal-flow barriers using the U.S.

Geological Survey’s modular three-dimensional finitedifference ground-water flow model: U.S. Geological

Survey Open-File Report 92-477, 32 p.

International Boundary and Water Commission Bulletin,

2000.

Kernodle, J.M., 1992, Results of simulations by a

preliminary numerical model of land subsidence in the

El Paso, Texas, area: U.S. Geological Survey WaterResources Investigations Report 92-4037, 35 p.

Kipp, K.L., Jr., 1987, HST3D—A computer code for

simulation of heat and solute transport in threedimensional ground-water flow systems: U.S.

Geological Survey Water-Resources Investigations

Report 86-4095, 517 p.

Kontis, A.L., and Mandle, R.J., 1988, Modifications of a

three-dimensional, ground-water flow model to account

for variable water density and effects of multiaquifer

wells: U.S. Geological Survey Water-Resources

Investigations Report 87-4265, 78 p.

Kuniansky, J., and Hillestad, J.G., 1980, Reservoir

simulation using bottom-hole pressure boundary

conditions: Society of Petroleum Engineers Journal,

v. 20, p. 473-486.

Land, L.F., and Armstrong, C.A., 1985, A preliminary

assessment of land-surface subsidence in the El Paso

area, Texas: U.S. Geological Survey Water-Resources

Investigations Report 85-4155, 96 p.

Langford, R.P., 2001, Segmentation of the Rio Grande

alluvial surface and evolution of the Rio Grande rift

[abs.]: Proceedings of the 53d Annual Meeting (April

29–May 2, 2001) of the Rocky Mountain Section,

Geological Society of America, Boulder, Colo.

Lee Wilson and Associates, 1985a, Report

3—Hydrogeology of the Hueco Basin: Prepared for the

Public Services Board, City of El Paso, Texas.

———1985b, Report 4—Technical framework for

evaluation of proposed Hueco Basin appropriations:

Prepared for the Public Services Board, City of El Paso,

Texas.

———1991, Research Memoranda 3—Solute transport

model for the Hueco Bolson Recharge Project:

Prepared for the Public Services Board, City of El Paso,

Texas.

Leggat, E.R., and Davis, M.E., 1966, Analog model study of

the Hueco Bolson near El Paso, Texas: Texas Water

Development Board Report 28, 26 p.

McDonald, M.G., 1984, Development of a multi-aquifer

well option for a modular ground-water flow model:

Proceedings of the National Water Well Association

Conference, Practical applications of ground-water

models, p. 786-796.

McDonald, M.G., and Harbaugh, A.W., 1988, A modular

three-dimensional finite-difference ground-water flow

model: U.S. Geological Survey Techniques of WaterResources Investigations, book 6, chap. A1.

McLean, J.S., 1970, Saline ground-water resources of the

Tularosa Basin, New Mexico: U.S. Department of the

Interior, Office of Saline Water Research and

Development Progress Report 561, 128 p.

Meyer, W.R., 1976, Digital model for simulated effects of

ground-water pumping in the Hueco Bolson, El Paso

area, Texas, New Mexico, and Mexico: U.S. Geological

Survey Water-Resources Investigations 58-75, 31 p.

Poeter, E.P., and Hill, M.C., 1998, Documentation of

UCODE, a computer code for universal inverse

modeling: U.S. Geological Survey Water-Resources

Investigations Report 98-4080, 116 p.

Pollock, D.W., 1994, User’s guide for MODPATH/

MODPATH-PLOT, version 3—A particle tracking

post-processing package for MODFLOW, the U.S.

Geological Survey finite-difference ground-water flow

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

94-464, 120 p.

Prickett, T.A., and Lonnquist, C.G., 1971, Selected digital

computer techniques for groundwater resource

evaluation: Champaign, Illinois State Water Survey

Bulletin 55, 62 p.

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.

Reilly, T.E., Franke, O.L., and Bennett, G.D., 1984, The

principle of superposition and its application in groundwater hydraulics: U.S. Geological Survey Open-File

Report 84-459, 36 p.

Riley, F.S., 1998, Mechanics of aquifer systems—The

scientific legacy of Joseph F. Poland, in Borchers, J.W.,

ed., Land subsidence—Case studies and current

research: Proceedings of the Dr. Joseph F. Poland

Symposium on Land Subsidence, Association of

Engineering Geologists Special Publication no. 8,

p. 13-27.

Robinson, T.W., 1958, Phreatophytes: U.S. Geological

Survey Water-Supply Paper 1423, 84 p.

Sayre, A.N., and Livingston, P., 1945, Ground-water

resources of the El Paso area, Texas: U.S. Geological

Survey Water-Supply Paper 919, 190 p.

Seager, W.R., 1982, Geology of Organ Mountains and

southern San Andres Mountains, New Mexico:

Socorro, New Mexico Bureau of Mines and Mineral

Resources Memoir 36.
55

Seager, W.R., Hawley, J.W., Kottlowski, F.E., and Kelley,

S.A., 1987, Geology of the east half of Las Cruces and

northeast El Paso 1 x 2 sheets (scale 1:125,000), New

Mexico: Socorro, New Mexico Bureau of Mines and

Mineral Resources, 5 sheets.

Slichter, C.S., 1905, Observations on the ground waters of

the Rio Grande Valley: U.S. Geological Survey WaterSupply Paper 141, 83 p.

Strain, W.S., 1969, Late Cenozoic strata of the El Paso area,

in Kottlowski, F.E., and Lemone, D.V., eds., Border

Stratigraphy Symposium: Socorro, New Mexico

Bureau of Mines and Mineral Resources Circular 104,

123 p.

Thiem, G., 1906, Hydrologische methoden: Gebhardt,

Leipzig, 56 p.

Trescott, P.C., and Larson, S.P., 1976, Finite-difference

model for aquifer simulation in two dimensions with

results of numerical experiments: U.S. Geological

Survey Techniques of Water-Resources Investigations,

book 6, chap. A1, 116 p.

U.S. Department of Commerce, U.S. Census Bureau, 2002,

accessed February 26, 2002, at http://www.census.gov/

prod/2002pubs/01statab/pop.pdf.

Waltemeyer, S.D., 2001, Estimates of mountain-front

streamflow available for potential recharge to the

Tularosa Basin, New Mexico: U.S. Geological Survey

Water-Resources Investigations Report 01-4013, 8 p.

Weeks, E.P., Weaver, H.L., Campbell, G.S., and Tanner,

B.D., 1987, Water use by salt cedar and by replacement

vegetation in the Pecos River flood plain between Acme

and Artesia, New Mexico: U.S. Geological Survey

Professional Paper 491-G, 33 p.

White, D.E., Baker, E.T., and Sperka, Roger, 1997,

Hydrology of the shallow aquifer and uppermost

semiconfined aquifer near El Paso, Texas: U.S.

Geological Survey Water-Resources Investigations

Report 97-4263, 37 p.

Wilkins, D.W., 1998, Summary of the Southwest Alluvial

Basins Regional Aquifer-System Analysis in parts of

Colorado, New Mexico, and Texas: U.S. Geological

Survey Professional Paper 1407-A, 49 p.

Zheng, Chunmiao, and Wang, P.P., 1998, MT3DMS—A

modular three-dimensional multispecies transport

model for simulation of advection, dispersion, and

chemical reactions of contaminants in groundwater

systems, Documentation and user’s guide, contract

report SERD P-99-1: U.S. Army Engineer Research

and Development Center, Vicksburg, Miss.
56

APPENDIX 1: MODIFICATIONS TO MODFLOWP

Changes made to modflowp.f for all modifications

 DATA CUNIT/'BCF ', 'WEL ', 'DRN ', 'RIV ', 'EVT ', 'TLK ', 'GHB ',

 & 'RCH ', 'SIP ', 'DE4 ', 'SOR ', 'OC ', 'PCG ', 'GFD ',

 & 'PAR ', 'HFB ', 'RES ', 'STR ', 'IBS ', 'CHD ', 'MAW ',

 & ' ', ' ', ' ', ' ', ' ', ' ', ' ',

 & ' ', ' ', ' ', ' ', ' ', ' ', ' ',

 & ' ', ' ', ' ', ' ', ' '/

... skip lines...

 IF (IUNIT(2).GT.0) CALL WEL5AL(ISUM,LENX,LCWELL,MXWELL,NWELLS,

 & IUNIT(2),IOUT,IWELCB,NWELVL,IWELAL,

 & IFREFM)

IF(IUNIT(21).GT.0) CALL MAW5AL(ISUM,LENX,MXMAW,NMAWS,

 1 IUNIT(21),IOUT,NLAY,LCRMAW,LCRMAL,IMAWCB)

... skip lines...

 IF (IUNIT(18).GT.0) CALL STR5AL(ISUM,LENX,LCSTRM,ICSTRM,MXSTRM,

 & NSTREM,IUNIT(18),IOUT,ISTCB1,

 & ISTCB2,NSS,NTRIB,NDIV,ICALC,CONST,

 & LCTBAR,LCTRIB,LCIVAR,LCFGAR,NSTROP)

... skip lines...

 IF (IUNIT(2).GT.0) CALL WEL5RP(X(LCWELL),NWELLS,MXWELL,

 & IUNIT(2),IOUT,NWELVL,IWELAL,

 & IFREFM)

 IF(IUNIT(21).GT.0) CALL MAW5RP(X(LCRMAW),

 1 X(LCRMAL),NMAWS,MXMAW,IUNIT(21),IOUT,NLAY)

... skip lines...

 IF (IUNIT(2).GT.0) CALL WEL5FM(NWELLS,MXWELL,X(LCRHS),

 & X(LCWELL),X(LCIBOU),NCOL,

 & NROW,NLAY,NWELVL)

 IF(IUNIT(21).GT.0) CALL MAW5FM(NMAWS,MXMAW,X(LCRHS),X(LCHCOF),

 1 X(LCIBOU),X(LCRMAW),X(LCRMAL),NCOL,NROW,

 2 NLAY,X(LCCR),X(LCDELC),X(LCDELR),X(LCHNEW),IOUT)

... skip lines...

 IF (IUNIT(18).GT.0) CALL STR5FM(NSTREM,X(LCSTRM),X(ICSTRM),

 & X(LCHNEW),X(LCHCOF),X(LCRHS)

 & ,X(LCIBOU),MXSTRM,NCOL,NROW,

 & NLAY,IOUT,NSS,X(LCTBAR),

 & NTRIB,X(LCTRIB),X(LCIVAR),

 & X(LCFGAR),ICALC,CONST,NSTROP,

 & X(LCSHNW),IUNIT(15),IP,ISN)

... skip lines...

 IF (IUNIT(2).GT.0) CALL WEL5BD(NWELLS,MXWELL,VBNM,VBVL,MSUM,

 & X(LCWELL),X(LCIBOU),DELT,

 & NCOL,NROW,NLAY,KKSTP,KKPER,

 & IWELCB,ICBCFL,X(LCBUFF),IOUT,

 & PERTIM,TOTIM,NWELVL,IWELAL)

 IF(IUNIT(21).GT.0) CALL MAW5BD(NMAWS,MXMAW,X(LCIBOU),X(LCRMAW),

 1 X(LCRMAL),NCOL,NROW,NLAY,X(LCHNEW),KSTP,

 2 KPER,IMAWCB,ICBCFL,X(LCBUFF),IOUT,MSUM,DELT,VBNM,VBVL)

... skip lines...
57

 IF (IUNIT(18).GT.0) CALL STR5BD(NSTREM,X(LCSTRM),X(ICSTRM),

 & X(LCIBOU),MXSTRM,X(LCHNEW),

 & NCOL,NROW,NLAY,DELT,VBVL,

 & VBNM,MSUM,KKSTP,KKPER,

 & ISTCB1,ISTCB2,ICBCFL,

 & X(LCBUFF),IOUT,NTRIB,NSS,

 & X(LCTRIB),X(LCTBAR),

 & X(LCIVAR),X(LCFGAR),ICALC,

 & CONST,IPTFLG,NSTROP)

... skip lines...

 CALL SEN1OT(IUHEAD,IOUT,NROW,NCOL,NLAY,NP,ISN,NH,X(LCNDER),

 & X(LCHNEW),PID,DID,IP,KPER,X(LCBUFF),KSTP,PERTIM,TOTIM,

 & X(LCX),X(LCDELR),X(LCDELC),X(LCIBOU),X(LCCOFF),

 & X(LCROFF),X(LCH),X(LCWT),X(LCHOBS),IPRINT,IFO,ITERP,

 & IPAR,X(LCRINT),X(LCJOFF),X(LCIOFF),X(LCMLAY),X(LCPR),

 & MOBS,NPER,X(LCB),X(LCLN),NQ,NQC,NQT,X(LCNQOB),

 & X(LCNQCL),X(LCIQOB),X(LCQCLS),X(LCIBT),MXBND,NBOUND,

 & X(LCBNDS),MXRIVR,NRIVER,X(LCRIVR),X(LCSHNW),LASTX,

 & ISCALS,X(LCTOFF),MXDRN,NDRAIN,X(LCDRAI),MXSTRM,NSTREM,

 & X(LCSTRM),X(ICSTRM),MAXM,KPRINT,JDRY,IDRY,NPR,X(LCWP),

 & MPR,X(LCPRM),X(LCIWPG),X(LCB1),IOUB,RSQ,RSQP,RSQO,

 & RSQOO,NPO,SOSC,SOSR,IPR,X(LCNIPR),X(LCWPF),ND,RSQF,

 & IOUYR,IOUHDS,IOUFLW,IOUPRI,IUNORM,X(LCWTQ),X(LCWTQS),

 & IOWTQ,NDMH,X(LCPV),X(LCRHS),X(LCCR),X(LCCC),X(LCCV),

 & X(LCBOT),X(LCTOP),NPTH,X(LCNPNT),NTT2,KTDIM,KTREV,

 & X(LCICLS),X(LCPRST),X(LCPOFF),X(LCSMAT),X(LCNLL),NSM,

 & X(LCTRPY),NMM,NZM,LZI1,X(LCSFAC),X(LCLZ),X(LCLM),

 & X(LCMATZ),NLLI1,X(LCST),X(LCTT2),IOUTT2,NRCHOP,

 & X(LCIRCH),X(LCRECH),X(LCSV),X(LCBANI),X(LCTHCK),NCLAY,

 & KTFLG,ADVSTP,X(LCLAYC),NSTROP)
58

Compute heads for multi-layer observations with dry cells: ssen1jz.f

Subroutine SSEN1U

This subroutine interpolates heads and accounts for dry cells, if necessary.

 REAL B, COFF, DELC, DELR, FACT, H, HD, HOBS, PR, PROP, RHS, RINT,

 & ROFF, TOFF, W, WT, X, ZERO, PROLD

...skip lines...

C---------IF THE OBSERVATION THIS IS TO BE SUBTRACTED FROM IS DRY, MAKE

C---------THIS ONE DRY, TOO

 N1 = NDER(5,N)

IF (N1.GT.0) THEN

 IF ((WT(N1).LT.ZERO.OR.COFF(N1).GE.5.)) THEN

 IDRY = IDRY + 1

 WT(N) = -ABS(WT(N))

 WRITE (IOUT,500) N, DID(N)

 GOTO 30

 ENDIF

 ENDIF

C---------CHECK FOR DRY OBSERVATIONS OR INTERPOLATIONS AFFECTED BY DRY

C---------CELLS

 DO 20 M = 1, MM

 KK = K

 IF (K.LT.0) KK = MLAY(M,ML)

 IF (KK.EQ.0) GOTO 30

 IF (LAYCON(KK).EQ.1 .OR. LAYCON(KK).EQ.3) THEN

 IF (IBOUND(JJ,II,KK).EQ.0) THEN

C*** dry cell in observation stack:

C*** more than 1 layer in observations

C*** change MLAY (move layer definitions forward in array)

C*** recalculate pr

C*** at least 1 non dry layer?

C*** goto 20

 IF (MLAY(M+1,ML).eq.0) GOTO 13

 PROLD = PR(M,ML)

 DO 11 MTMP = M, MM-1

 IF (MLAY(MTMP,ML).eq.0) GOTO 12

 MLAY(MTMP,ML) = MLAY(MTMP+1,ML)

 PR(MTMP,ML) = PR(MTMP+1,ML)/(1-PROLD)

 IF (PR(MTMP+1,ML).eq.0) PR(MTMP,ML) = 0

 11 CONTINUE

 MLAY(MM,ML) = 0

 PR(MM,ML) = 0

 12 IF (MLAY(1,ML).gt.0) GOTO 20

 13 IDRY = IDRY + 1

 WT(N) = -ABS(WT(N))

 WRITE (IOUT,500) N, DID(N)

 GOTO 30

 ELSEIF ((RINT(2,N).NE.ZERO.AND.IBOUND(JJ+JO,II,KK)

 & .EQ.0) .OR.

 & (RINT(3,N).NE.ZERO.AND.IBOUND(JJ,II+IO,KK)

 & .EQ.0) .OR.

 & (RINT(4,N).NE.ZERO.AND.IBOUND(JJ+JO,II+IO,KK)

 & .EQ.0)) THEN

C*** dry cell in adjacent stack:

C*** IF (MM.GT.1 .OR. TOFF(N).GT.ZERO) THEN

C*** multi-layer observation

 IF (MM.GT.1) THEN

C*** adjust PR and ML arrays as above

 IF (MLAY(M+1,ML).eq.0) GOTO 16

 PROLD = PR(M,ML)

 DO 14 MTMP = M, MM-1
59

 IF (MLAY(MTMP,ML).eq.0) GOTO 15

 MLAY(MTMP,ML) = MLAY(MTMP+1,ML)

 PR(MTMP,ML) = PR(MTMP+1,ML)/(1-PROLD)

 IF (PR(MTMP+1,ML).eq.0) PR(MTMP,ML) = 0

 14 CONTINUE

 MLAY(MM,ML) = 0

 PR(MM,ML) = 0

 15 IF (MLAY(1,ML).gt.0) GOTO 20

 ENDIF

 16 IF (MM.GT.1 .OR. TOFF(N).GT.ZERO) THEN

 IDRY = IDRY + 1

 WT(N) = -ABS(WT(N))

 WRITE (IOUT,500) N, DID(N)

 GOTO 30

 ENDIF

Compute leakage from streams to underlying active cells: str5p.f

Subroutine STR5AL

This subroutine allocates array storage for streams.

 SUBROUTINE STR5AL(ISUM,LENX,LCSTRM,ICSTRM,MXSTRM,NSTREM,IN,IOUT,

 & ISTCB1,ISTCB2,NSS,NTRIB,NDIV,ICALC,CONST,LCTBAR,

 & LCTRIB,LCIVAR,LCFGAR,NSTROP)

... skip lines...

 READ (IN,505) MXSTRM, NSS, NTRIB, NDIV, ICALC, CONST, ISTCB1,

 & ISTCB2,NSTROP

505 FORMAT (5I10,F10.0,3I10)

... skip lines...

 540 FORMAT (' ***X ARRAY MUST BE DIMENSIONED LARGER***')

C

C3------CHECK TO SEE THAT OPTION IS LEGAL.

 IF(NSTROP.GE.1.AND.NSTROP.LE.3) GO TO 560

C

C3A-----IF ILLEGAL PRINT A MESSAGE AND ABORT SIMULATION

 WRITE(IOUT,550)

 550 FORMAT(1X,'ILLEGAL OPTION CODE. SIMULATION ABORTING')

 STOP

C

C4------PRINT OPTION CODE.

 560 IF(NSTROP.EQ.1) WRITE(IOUT,565)

 565 FORMAT(1X,'OPTION 1 -- LEAKAGE TO DEFINED LAYER')

 IF(NSTROP.EQ.3) WRITE(IOUT,570)

 570 FORMAT(1X,'OPTION 3 -- LEAKAGE TO HIGHEST ACTIVE NODE IN EACH',

 1 ' VERTICAL COLUMN')

C

C10-----RETURN.

 RETURN

 END

Subroutine STR5FM

This subroutine adds stream terms to RHS and HCOF if flow occurs in model cell.

 SUBROUTINE STR5FM(NSTREM,STRM,ISTRM,HNEW,HCOF,RHS,IBOUND,MXSTRM,

 & NCOL,NROW,NLAY,IOUT,NSS,ITRBAR,NTRIB,ARTRIB,

 & IDIVAR,NDFGAR,ICALC,CONST,NSTROP,SHNW,IU,IP,ISN)

... skip lines ...
60

C12----DETERMINE LEAKAGE THROUGH STREAMBED.

C IF DRY CELL AND OPTION 3, CHECK FOR HIGHEST ACTIVE CELL

 IF ((IBOUND(IC,IR,IL).LT.1).AND.(NSTROP.EQ.1)) THEN

 FLOBOT = ZERO

 ELSE

IF((IBOUND(IC,IR,IL).LE.0).AND.(NSTROP.EQ.3)) THEN

C12A------IF OPTION IS 3 LEAKAGE IS INTO HIGHEST INTERNAL CELL.

C CANNOT PASS THROUGH CONSTANT HEAD NODE

 DO 42 NIL=IL+1,NLAY

C

C12B-----IF CELL IS: CONSTANT HEAD MOVE ON TO NEXT STREAM REACH

C INACTIVE MOVE DOWN A LAYER

C ACTIVE SET LEAKAGE TO CONSTANT

 IF(IBOUND(IC,IR,NIL)) 43,42,44

 42 CONTINUE

C--------NO ACTIVE CELLS FOUND, NO LEAKAGE

 43 FLOBOT=0.

 GO TO 45

 ENDIF

 44 IF (FLOWIN.LE.ZERO) HSTR = STRM(5,L)

 CSTR = STRM(3,L)

 SBOT = STRM(4,L)

 H = HNEW(IC,IR,IL)

C*** set head to top active layer if dry cell

 IF(IBOUND(IC,IR,IL).EQ.0) H=HNEW(IC,IR,NIL)

C-----ADDED FOR PARAMETER ESTIMATION

 IF (IU.GT.0 .AND. IP.NE.0) H = SHNW(IC,IR,IL)

IF(IU.GT.0.AND.IP.NE.0.AND.IBOUND(IC,IR,IL).EQ.0)

 + H=SHNW(IC,IR,NIL)

... skip lines...

C16-----STREAMFLOW OUT EQUALS STREAMFLOW IN MINUS LEAKAGE.

IF ((IBOUND(IC,IR,IL).LT.1).AND.(NSTROP.EQ.1)) FLOBOT = ZERO

 ENDIF

 ENDIF

45 FLOWOT = FLOWIN - FLOBOT

 IF (ISTSG.GT.1 .AND. NREACH.EQ.1) STRM(9,LL) = ARTRIB(IFLG)

C

C17----STORE STREAM INFLOW, OUTFLOW AND LEAKAGE FOR EACH REACH.

 STRM(9,L) = FLOWOT

 STRM(10,L) = FLOWIN

 STRM(11,L) = FLOBOT

C

C18----RETURN TO STEP 3 IF STREAM INFLOW IS LESS THAN OR EQUAL TO ZERO

C AND LEAKAGE IS GREATER THAN OR EQUAL TO ZERO OR IF CELL

C IS NOT ACTIVE--IBOUND IS LESS THAN OR EQUAL TO ZERO--

C & NSTROP=1.

 IF ((IBOUND(IC,IR,IL).GT.0).OR.(NSTROP.EQ.3)) THEN

C*** reset top active layer if dry cell

 IF(IBOUND(IC,IR,IL).EQ.0) IL=NIL

IF (FLOWIN.GT.ZERO .OR. FLOBOT.LT.ZERO) THEN

C

C19------IF HEAD > BOTTOM THEN ADD TERMS TO RHS AND HCOF.

 IF (IQFLG.LT.1) THEN

C--------FOR ADJOINT STATES, ONLY CALCULATE CONTRIBUTION TO HCOF

 IF (IU.EQ.0 .OR. IP.EQ.0 .OR. ISN.LT.0)

 & RHS(IC,IR,IL) = RHS(IC,IR,IL) - CSTR*HSTR

 HCOF(IC,IR,IL) = HCOF(IC,IR,IL) - CSTR

 ELSE

C

C20------IF HEAD < BOTTOM THEN ADD TERM ONLY TO RHS.

C-------FOR ADJOINT STATES, ONLY CALCULATE CONTRIBUTION TO HCOF

IF (IU.EQ.0 .OR. IP.EQ.0 .OR. ISN.LT.0) THEN

 RHS(IC,IR,IL)=RHS(IC,IR,IL) - FLOBOT
61

 ENDIF

 ENDIF

 ENDIF

 ENDIF

 50 CONTINUE

Subroutine STR5BD

This subroutine calculates volumetric budget for streams.

 SUBROUTINE STR5BD(NSTREM,STRM,ISTRM,IBOUND,MXSTRM,HNEW,NCOL,NROW,

 & NLAY,DELT,VBVL,VBNM,MSUM,KSTP,KPER,ISTCB1,

 & ISTCB2,ICBCFL,BUFF,IOUT,NTRIB,NSS,ARTRIB,ITRBAR,

 & IDIVAR,NDFGAR,ICALC,CONST,IPTFLG,NSTROP)

... skip lines...

C14----DETERMINE LEAKAGE THROUGH STREAMBED.

IF ((IBOUND(IC,IR,IL).LT.1).AND.(NSTROP.EQ.1)) THEN

 FLOBOT = ZERO

 ELSE

IF((IBOUND(IC,IR,IL).LE.0).AND.(NSTROP.EQ.3)) THEN

 DO 72 NIL=IL+1,NLAY

C

C18A-----IF CELL IS: CONSTANT HEAD MOVE ON TO NEXT STREAM REACH

C INACTIVE MOVE DOWN A LAYER

C ACTIVE SET LEAKAGE TO CONSTANT

 IF(IBOUND(IC,IR,NIL)) 73,72,75

 72 CONTINUE

C--------NO ACTIVE CELLS FOUND, NO LEAKAGE

 73 FLOBOT=0.

 GO TO 77

 75 IF (STRM(3,L).ne.0)

 & WRITE (IOUT,400) NIL, ISTRM(4,L), ISTRM(5,L), IR, IC

 400 FORMAT (/,5X,'LEAKAGE TO LAYER',I5,' FROM STREAM SEG',

 & I6,' STREAM REACH',I6,' IN ROW',I5,' COLUMN',I5)

 ENDIF

 IF (FLOWIN.LE.ZERO) HSTR = STRM(5,L)

CSTR = STRM(3,L)

 SBOT = STRM(4,L)

H = HNEW(IC,IR,IL)

C*** set head to top active layer if dry cell

 IF(IBOUND(IC,IR,IL).EQ.0) H=HNEW(IC,IR,NIL)

... skip lines...

C18----STREAMFLOW OUT EQUALS STREAMFLOW IN MINUS LEAKAGE.

 IF ((IBOUND(IC,IR,IL).LT.1).AND.(NSTROP.EQ.1))

 & FLOBOT = ZERO

 ENDIF

 ENDIF

 77 FLOWOT = FLOWIN - FLOBOT

 IF (ISTSG.GT.1 .AND. NREACH.EQ.1) STRM(9,LL) = ARTRIB(IFLG)

C

C19----STORE STREAM INFLOW, OUTFLOW AND LEAKAGE FOR EACH REACH.

 STRM(9,L) = FLOWOT

 STRM(10,L) = FLOWIN

 STRM(11,L) = FLOBOT

C

C20----IF LEAKAGE FROM STREAMS IS TO BE SAVED THEN ADD RATE TO BUFFER.

C------OPTION=3; LEAKAGE IS INTO HIGHEST CELL IN A VERTICAL COLUMN

C------THAT IS NOT NO FLOW. IF NO ACTIVE CELLS EXIST THEN ZERO STREAM

C------LEAKAGE

 IF ((IBD.EQ.1).AND.(FLOBOT.NE.0.0)) THEN

 IF(IBOUND(IC,IR,IL).GT.0) THEN
62

 BUFF(IC,IR,IL)=BUFF(IC,IR,IL) + FLOBOT

 ELSE

 BUFF(IC,IR,NIL)=BUFF(IC,IR,NIL) + FLOBOT

 ENDIF

 ENDIF

C

C21----DETERMINE IF FLOW IS INTO OR OUT OF MODEL CELL.

... skip lines...

C29-----SAVE STREAMFLOWS OUT OF EACH REACH ON DISK.

 DO 120 L = 1, NSTREM

 IC = ISTRM(3,L)

 IR = ISTRM(2,L)

 IL = ISTRM(1,L)

 IF (IBOUND(IC,IR,IL).GT.0)

 & BUFF(IC,IR,IL) = BUFF(IC,IR,IL) + STRM(9,L)

C*** reset top active layer if dry cell

 IF ((IBOUND(IC,IR,IL).LE.0).AND.(NSTROP.EQ.3)) THEN

 DO 118 NIL=IL+1,NLAY

C

C29A-----IF CELL IS: CONSTANT HEAD MOVE ON TO NEXT STREAM REACH

C INACTIVE MOVE DOWN A LAYER

 IF(IBOUND(IC,IR,NIL)) 118,119,118

 118 CONTINUE

 GO TO 120

 119 BUFF(IC,IR,NIL) = BUFF(IC,IR,NIL) + STRM(9,L)

 ENDIF

 120 CONTINUE

 CALL UBUDSV(KSTP,KPER,STRTXT,ISTCB2,BUFF,NCOL,NROW,NLAY,IOUT)

Subroutine PAR1AQ

This subroutine calculates final iteration parameters and updates parameters.

 CALL SSEN1V(NQ,NQC,NQT,NQOB,NQCL,IQOB,QCLS,IBT,MXBND,NBOUND,

 & BNDS,MXRIVR,NRIVER,RIVR,SHNW,IP,HNEW,NCOL,NROW,NLAY,

 & IOUT,IBOUND,NPER,KPER,NH,X,DID,NP,H,B,LN,TOFF,

 & MXDRN,NDRAIN,DRAI,MXSTRM,NSTREM,STRM,ISTRM,ISN,

 & WTQ,NDMH,NSTROP)

Subroutine SEN1OT

This subroutine prints data for observed heads and flows.

 SUBROUTINE SEN1OT(IUHEAD,IOUT,NROW,NCOL,NLAY,NP,ISN,NH,NDER,HNEW,

 & PID,DID,IP,KPER,BUFF,KSTP,PERTIM,TOTIM,X,DELR,

 & DELC,IBOUND,COFF,ROFF,H,WT,HOBS,IPRINT,IFO,

 & ITERP,IPAR,RINT,JOFF,IOFF,MLAY,PR,MOBS,NPER,B,

 & LN,NQ,NQC,NQT,NQOB,NQCL,IQOB,QCLS,IBT,MXBND,

 & NBOUND,BNDS,MXRIVR,NRIVER,RIVR,SHNW,LASTX,

 & ISCALS,TOFF,MXDRN,NDRAIN,DRAI,MXSTRM,NSTREM,

 & STRM,ISTRM,MAXM,KPRINT,JDRY,IDRY,NPR,WP,MPR,PRM,

 & IWPG,B1,IOUB,RSQ,RSQP,RSQO,RSQOO,NPO,SOSC,SOSR,

 & IPR,NIPR,WPF,ND,RSQF,IOUYR,IOUHDS,IOUFLW,IOUPRI,

 & IUNORM,WTQ,WTQS,IOWTQ,NDMH,PV,RHS,CR,CC,CV,BOT,

 & TOP,NPTH,NPNT,NTT2,KTDIM,KTREV,ICLS,PRST,POFF,

 & SMAT,NLL,NSM,TRPY,NMM,NZM,LZI1,SFAC,LZ,LM,MATZ,

 & NLLI1,ST,TT2,IOUTT2,NRCHOP,IRCH,RECH,SV,BANIV,

 & THCK,NCLAY,KTFLG,ADVSTP,LAYC,NSTROP)

... skip lines...

 IF (NQ.GT.0) CALL SSEN1V(NQ,NQC,NQT,NQOB,NQCL,IQOB,QCLS,IBT,
63

 & MXBND,NBOUND,BNDS,MXRIVR,NRIVER,RIVR,

 & SHNW,IP,HNEW,NCOL,NROW,NLAY,IOUT,

 & IBOUND,NPER,KPER,NH,X,DID,NP,H,B,LN,

 & TOFF,MXDRN,NDRAIN,DRAI,MXSTRM,NSTREM,

 & STRM,ISTRM,ISN,WTQ,NDMH,NSTROP)

Subroutine SSEN1V

This subroutine saves simulated flows and calculates sensitivities.

CHANGE 28.02.96:ARGUMENTS ADDED SUBROUTINE AND DIMENSION STATEMENT

SUBROUTINE SSEN1V(NQ,NQC,NQT,NQOB,NQCL,IQOB,QCLS,IBT,MXBND,

 & NBOUND,BNDS,MXRIVR,NRIVER,RIVR,SHNW,IP,HNEW,NCOL,NROW,

 & NLAY,IOUT,IBOUND,NPER,KPER,NH,X,DID,NP,H,B,LN,

 & TOFF,MXDRN,NDRAIN,DRAI,MXSTRM,NSTREM,STRM,ISTRM,

 & ISN,WTQ,NDMH,NSTROP)

... skip lines...

C-------------ASSIGN VARIABLE VALUES

 IF (IBT1.EQ.3) THEN

 K = ISTRM(1,NB)

 I = ISTRM(2,NB)

 J = ISTRM(3,NB)

C IF DRY CELL AND OPTION 3, CHECK FOR HIGHEST ACTIVE CELL

 IF ((IBOUND(J,I,K).LT.1).AND.(NSTROP.EQ.1)) GOTO 30

 IF ((IBOUND(J,I,K).LE.0).AND.(NSTROP.EQ.3)) THEN

C------IF OPTION IS 3 LEAKAGE IS INTO HIGHEST INTERNAL CELL.

C CANNOT PASS THROUGH CONSTANT HEAD NODE

 DO 5 NIL=K+1,NLAY

C

C-----IF CELL IS: CONSTANT HEAD MOVE ON TO NEXT STREAM REACH

C INACTIVE MOVE DOWN A LAYER

C ACTIVE SET LEAKAGE TO CONSTANT

 IF(IBOUND(J,I,NIL)) 6,5,7

 5 CONTINUE

C--------NO ACTIVE CELLS FOUND, NO LEAKAGE

 6 GO TO 30

 ENDIF

 ENDIF

7 IF (IP.EQ.0) HHNEW = HNEW(J,I,K)

 IF (IP.GT.0) HHNEW = SHNW(J,I,K)

APPENDIX 2: MODIFICATIONS TO MODFLOW

Changes made to modflw96.f for all modifications

 DATA CUNIT/'BCF ','WEL ','DRN ','RIV ','EVT ','TLK ','GHB ',

 1 'RCH ','SIP ','DE4 ','SOR ','OC ','PCG ','GFD ',

 2 'MAW ','HFB ','RES ','STR ','IBS ','CHD ','FHB ',

 3 ' ',' ',' ',' ',' ',' ',' ',

 4 ' ',' ',' ',' ',' ',' ',' ',

 5 ' ',' ',' ',' ',' '/

... skip lines...

 IF(IUNIT(2).GT.0) CALL WEL5AL(ISUM,LENX,LCWELL,MXWELL,NWELLS,

 1 IUNIT(2),IOUT,IWELCB,NWELVL,IWELAL,IFREFM)

C*** multi-layer wells

 IF(IUNIT(15).GT.0) CALL MAW5AL(ISUM,LENX,MXMAW,NMAWS,
64

 1 IUNIT(15),IOUT,NLAY,LCRMAW,LCRMAL,IMAWCB)

 ... skip lines...

 IF(IUNIT(18).GT.0) CALL STR1AL(ISUM,LENX,LCSTRM,ICSTRM,MXSTRM,

 1 NSTREM,IUNIT(18),IOUT,ISTCB1,ISTCB2,NSS,NTRIB,

 2 NDIV,ICALC,CONST,LCTBAR,LCTRIB,LCIVAR,LCFGAR,

3 NSTROP)

 ... skip lines...

 IF(IUNIT(2).GT.0) CALL WEL5RP(X(LCWELL),NWELLS,MXWELL,IUNIT(2),

 1 IOUT,NWELVL,IWELAL,IFREFM)

C*** multi-layer wells

 IF(IUNIT(15).GT.0) CALL MAW5RP(X(LCRMAW),

 1 X(LCRMAL),NMAWS,MXMAW,IUNIT(15),IOUT,NLAY)

... skip lines...

 IF(IUNIT(2).GT.0) CALL WEL5FM(NWELLS,MXWELL,X(LCRHS),X(LCWELL),

 1 X(LCIBOU),NCOL,NROW,NLAY,NWELVL)

C*** multi-layer wells

 IF(IUNIT(15).GT.0) CALL MAW5FM(NMAWS,MXMAW,X(LCRHS),X(LCHCOF),

 1 X(LCIBOU),X(LCRMAW),X(LCRMAL),NCOL,NROW,

 2 NLAY,X(LCCR),X(LCDELC),X(LCDELR),X(LCHNEW),IOUT)

... skip lines...

 IF(IUNIT(18).GT.0) CALL STR1FM(NSTREM,X(LCSTRM),X(ICSTRM),

 1 X(LCHNEW),X(LCHCOF),X(LCRHS),X(LCIBOU),

 2 MXSTRM,NCOL,NROW,NLAY,IOUT,NSS,X(LCTBAR),

 3 NTRIB,X(LCTRIB),X(LCIVAR),X(LCFGAR),ICALC,CONST,

4 NSTROP)

... skip lines...

 IF(IUNIT(2).GT.0) CALL WEL5BD(NWELLS,MXWELL,VBNM,VBVL,MSUM,

 1 X(LCWELL),X(LCIBOU),DELT,NCOL,NROW,NLAY,KKSTP,KKPER,IWELCB,

 1 ICBCFL,X(LCBUFF),IOUT,PERTIM,TOTIM,NWELVL,IWELAL)

C*** multi-layer wells

 IF(IUNIT(15).GT.0) CALL MAW5BD(NMAWS,MXMAW,X(LCIBOU),X(LCRMAW),

 1 X(LCRMAL),NCOL,NROW,NLAY,X(LCHNEW),KSTP,

 2 KPER,IMAWCB,ICBCFL,X(LCBUFF),IOUT,MSUM,DELT,VBNM,VBVL)

... skip lines...

 IF(IUNIT(18).GT.0) CALL STR1BD(NSTREM,X(LCSTRM),X(ICSTRM), STR1

 1 X(LCIBOU),MXSTRM,X(LCHNEW),NCOL,NROW,NLAY,DELT,VBVL,VBNM,MSUM, STR1

 2 KKSTP,KKPER,ISTCB1,ISTCB2,ICBCFL,X(LCBUFF),IOUT,NTRIB,NSS, STR1

 3 X(LCTRIB),X(LCTBAR),X(LCIVAR),X(LCFGAR),ICALC,CONST,IPTFLG,

4 NSTROP)

Compute leakage from streams to underlying active cells: str1.f

Subroutine STR1AL

This subroutine allocates array storage for streams.

 SUBROUTINE STR1AL(ISUM,LENX,LCSTRM,ICSTRM,MXSTRM,NSTREM,IN,

 1 IOUT,ISTCB1,ISTCB2,NSS,NTRIB,NDIV,ICALC,CONST,

 2 LCTBAR,LCTRIB,LCIVAR,LCFGAR,NSTROP)

 ... skip lines...

C2------ READ MXSTRM, NSS, NTRIB, ISTCB1, AND ISTCB2.

 100 READ(IN,3)MXSTRM,NSS,NTRIB,NDIV,ICALC,CONST,ISTCB1,ISTCB2,NSTROP
65

 3 FORMAT(5I10,F10.0,3I10)

 ... skip lines...

 10 FORMAT(1X,' ***X ARRAY MUST BE DIMENSIONED LARGER***')

C

C3------CHECK TO SEE THAT OPTION IS LEGAL.

 IF(NSTROP.GE.1.AND.NSTROP.LE.3) GO TO 250

C

C3A-----IF ILLEGAL PRINT A MESSAGE AND ABORT SIMULATION

 WRITE(IOUT,11)

 11 FORMAT(1X,'ILLEGAL OPTION CODE. SIMULATION ABORTING')

 STOP

C

C4------PRINT OPTION CODE.

 250 IF(NSTROP.EQ.1) WRITE(IOUT,12)

 12 FORMAT(1X,'OPTION 1 -- LEAKAGE TO DEFINED LAYER')

 IF(NSTROP.EQ.3) WRITE(IOUT,13)

 13 FORMAT(1X,'OPTION 3 -- LEAKAGE TO HIGHEST ACTIVE NODE IN EACH',

 1 ' VERTICAL COLUMN')

C10-----RETURN.

 RETURN

 END

Subroutine STR1FM

This subroutine adds stream terms to RHS and HCOF if flow occurs in model cell.

 SUBROUTINE STR1FM(NSTREM,STRM,ISTRM,HNEW,HCOF,RHS,IBOUND,MXSTRM,

 1 NCOL,NROW,NLAY,IOUT,NSS,ITRBAR,NTRIB,ARTRIB,

 2 IDIVAR,NDFGAR,ICALC,CONST,NSTROP)

... skip lines...

C12----DETERMINE LEAKAGE THROUGH STREAMBED. C

C IF DRY CELL AND OPTION 3, CHECK FOR HIGHEST ACTIVE CELL

 IF((IBOUND(IC,IR,IL).LE.0).AND.(NSTROP.EQ.1)) GO TO 315

 IF((IBOUND(IC,IR,IL).LE.0).AND.(NSTROP.EQ.3)) THEN

C12A------IF OPTION IS 3 LEAKAGE IS INTO HIGHEST INTERNAL CELL.

C CANNOT PASS THROUGH CONSTANT HEAD NODE

 DO 3101 NIL=IL+1,NLAY

C

C12B-----IF CELL IS: CONSTANT HEAD MOVE ON TO NEXT STREAM REACH

C INACTIVE MOVE DOWN A LAYER

C ACTIVE SET LEAKAGE TO CONSTANT

 IF(IBOUND(IC,IR,NIL)) 3102,3101,311

 3101 CONTINUE

C--------NO ACTIVE CELLS FOUND, NO LEAKAGE

 3102 FLOBOT=0.

 GO TO 320

 ENDIF

 311 IF(FLOWIN.LE.0.) HSTR=STRM(5,L)

 CSTR=STRM(3,L)

 SBOT=STRM(4,L)

 H=HNEW(IC,IR,IL)

C*** set head to top active layer if dry cell

 IF(IBOUND(IC,IR,IL).EQ.0) H=HNEW(IC,IR,NIL)

 T=HSTR-SBOT

... skip lines...

C16-----STREAMFLOW OUT EQUALS STREAMFLOW IN MINUS LEAKAGE.

315 IF((IBOUND(IC,IR,IL).LE.0).AND.(NSTROP.EQ.1)) FLOBOT=0.

 320 FLOWOT=FLOWIN-FLOBOT
66

 IF((ISTSG.GT.1).AND.(NREACH.EQ.1)) STRM(9,LL)=ARTRIB(IFLG)

C

C17----STORE STREAM INFLOW, OUTFLOW AND LEAKAGE FOR EACH REACH.

 STRM(9,L)=FLOWOT

 STRM(10,L)=FLOWIN

 STRM(11,L)=FLOBOT

C

C18----RETURN TO STEP 3 IF STREAM INFLOW IS LESS THAN OR EQUAL TO ZERO

C AND LEAKAGE IS GREATER THAN OR EQUAL TO ZERO, OR IF CELL

C IS NOT ACTIVE--IBOUND IS LESS THAN OR EQUAL TO ZERO--

C & NSTROP=1.

 IF((IBOUND(IC,IR,IL).LE.0).AND.(NSTROP.EQ.1)) GO TO 500

 IF((FLOWIN.LE.0.0).AND.(FLOBOT.GE.0.0)) GO TO 500

C

C19------IF HEAD > BOTTOM THEN ADD TERMS TO RHS AND HCOF.

C*** reset top active layer if dry cell

 IF(IBOUND(IC,IR,IL).EQ.0) IL=NIL

 IF(IQFLG.GT.0) GO TO 400

Subroutine STR1BD

This subroutine calculates volumetric budget for streams.

 SUBROUTINE STR1BD(NSTREM,STRM,ISTRM,IBOUND,MXSTRM,HNEW,NCOL,NROW,

 1 NLAY,DELT,VBVL,VBNM,MSUM,KSTP,KPER,ISTCB1,ISTCB2,ICBCFL,BUFF,

 2 IOUT,NTRIB,NSS,ARTRIB,ITRBAR,IDIVAR,NDFGAR,ICALC,CONST,IPTFLG,

 3 NSTROP)

... skip lines...

C14----DETERMINE LEAKAGE THROUGH STREAMBED.

 IF((IBOUND(IC,IR,IL).LE.0).AND.(NSTROP.EQ.1)) GO TO 315

 IF((IBOUND(IC,IR,IL).LE.0).AND.(NSTROP.EQ.3)) THEN

 DO 311 NIL=IL+1,NLAY

C

C18A-----IF CELL IS: CONSTANT HEAD MOVE ON TO NEXT STREAM REACH

C INACTIVE MOVE DOWN A LAYER

C ACTIVE SET LEAKAGE TO CONSTANT

 IF(IBOUND(IC,IR,NIL)) 3111,311,3112

 311 CONTINUE

C--------NO ACTIVE CELLS FOUND, NO LEAKAGE

 3111 FLOBOT=0.

 GO TO 320

 3112 IF (STRM(3,L).ne.0)

 & WRITE (IOUT,900) NIL, ISTRM(4,L), ISTRM(5,L), IR, IC

 900 FORMAT (/,5X,'LEAKAGE TO LAYER',I5,' FROM STREAM SEG',

 & I6,' STREAM REACH',I6,' IN ROW',I5,' COLUMN',I5)

 ENDIF

 IF(FLOWIN.LE.0.0) HSTR=STRM(5,L)

 CSTR=STRM(3,L)

 SBOT=STRM(4,L)

 H=HNEW(IC,IR,IL)

C*** set head to top active layer if dry cell

 IF(IBOUND(IC,IR,IL).EQ.0) H=HNEW(IC,IR,NIL)

 T=HSTR-SBOT

... skip lines...

C18----STREAMFLOW OUT EQUALS STREAMFLOW IN MINUS LEAKAGE.

315 IF((IBOUND(IC,IR,IL).LE.0).AND.(NSTROP.EQ.1)) FLOBOT=0.

 320 FLOWOT=FLOWIN-FLOBOT

 IF((ISTSG.GT.1).AND.(NREACH.EQ.1)) STRM(9,LL)=ARTRIB(IFLG)

C

C19----STORE STREAM INFLOW, OUTFLOW AND LEAKAGE FOR EACH REACH.

 STRM(9,L)=FLOWOT
67

 STRM(10,L)=FLOWIN

 STRM(11,L)=FLOBOT

C

C20----IF LEAKAGE FROM STREAMS IS TO BE SAVED THEN ADD RATE TO BUFFER.

C------OPTION=3; LEAKAGE IS INTO HIGHEST CELL IN A VERTICAL COLUMN

C------THAT IS NOT NO FLOW. IF NO ACTIVE CELLS EXIST THEN ZERO STREAM

C------LEAKAGE

 IF ((IBD.EQ.1).AND.(FLOBOT.NE.0.0)) THEN

 IF (IBOUND(IC,IR,IL).GT.0)

 BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+FLOBOT

 ELSE

 BUFF(IC,IR,NIL)=BUFF(IC,IR,NIL)+FLOBOT

 ENDIF

 ENDIF

C

C21----DETERMINE IF FLOW IS INTO OR OUT OF MODEL CELL.

... skip lines...

C29-----SAVE STREAMFLOWS OUT OF EACH REACH ON DISK. C

 DO 615 L=1,NSTREM

 IC=ISTRM(3,L)

 IR=ISTRM(2,L)

 IL=ISTRM(1,L)

 IF((IBOUND(IC,IR,IL).LE.0).AND.(NSTROP.EQ.1)) GO TO 615

 IF (IBOUND(IC,IR,IL).GT.0)

 & BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+STRM(9,L)

C*** set layer to top active cell

 IF (IBOUND(IC,IR,IL).EQ.0.AND.NSTROP.EQ.3)

 & BUFF(IC,IR,NIL)=BUFF(IC,IR,NIL)+STRM(9,L)

 615 CONTINUE
68

APPENDIX 3: MULTI-AQUIFER WELL PACKAGE

SUBROUTINE MAW5AL(ISUM,LENX,MXMAW,NMAWS,IN,IOUT,

 1 NLAY,LCRMAW,LCRMAL,IMAWCB)

C

C-----VERSION 5.0 OCT2001 MAW5AL

C ******************************************************************

C ALLOCATE ARRAY STORAGE FOR MULTI-AQUIFER WELL PACKAGE

C ******************************************************************

C

C1------IDENTIFY PACKAGE AND INITIALIZE NMAWS

 WRITE(IOUT,1)

 1 FORMAT(1H0,’MAW5 -- MULTI-AQUIFER WELL PACKAGE VERSION 5’)

NMAWS=0

C

C2------READ MAX NUMBER OF MAWS AND UNIT OR FLAG FOR

C2------CELL-BY-CELL FLOW TERMS.

 READ(IN,2) MXMAW,IMAWCB

 2 FORMAT(2I10)

 WRITE(IOUT,3) MXMAW

 3 FORMAT(1H ,’MAXIMUM OF’,I5,’ MULTI-AQUIFER WELL CATEGORIES’)

 IF(IMAWCB.GT.0) WRITE(IOUT,9) IMAWCB

 9 FORMAT(1X,’CELL-BY-CELL FLOW WILL BE SAVED ON UNIT’,I3)

 IF(IMAWCB.LT.0) WRITE(IOUT,8)

 8 FORMAT(1X,’CELL-BY-CELL FLOWS WILL BE PRINTED WHEN ICBCFL NOT 0’)

C

C3------ALLOCATE SPACE FOR ARRAYS RMAW AND RMAL.

 ISOLD=ISUM

 LCRMAW=ISUM

 ISUM=ISUM+5*MXMAW

 LCRMAL=ISUM

 ISUM=ISUM+4*NLAY*MXMAW

 ISP=ISUM-ISOLD

C

C4------PRINT NUMBER OF WORDS IN X ARRAY USED BY MAW PACKAGE.

 WRITE(IOUT,4) ISP

 4 FORMAT(1X,I6,’ ELEMENTS IN X ARRAY ARE USED FOR MAWS’)

 ISUM1=ISUM-1

 WRITE(IOUT,5) ISUM1,LENX

 5 FORMAT(1X,I6,’ ELEMENTS OF X ARRAY USED OUT OF ’,I13)

 IF(ISUM1.GT.LENX) WRITE(IOUT,6)

 6 FORMAT(1X,’ ***X ARRAY MUST BE DIMENSIONED LARGER***’)

C

C5------RETURN

 RETURN

 END

C
69

 SUBROUTINE MAW5RP(RMAW,RMAL,NMAWS,MXMAW,IN,IOUT,

 1 NLAY)

C-----VERSION 5.0 OCT2001 MAW1RP

C ******************************************************************

C READ MULTI-AQUIFER WELL LOCATIONS AND STRESS RATES

C modified by RMY to omit categories & read data from single line

C ******************************************************************

C

C SPECIFICATIONS:

C ------------------------------------------------------------------

 DIMENSION RMAW(5,MXMAW),RMAL(4,NLAY,MXMAW)

INTEGER LAYERS(10)

C -----------------------------------------------------------------

C

C1------READ ITMP(# OF MAW WELLS OR FLAG SAYING REUSE MAW DATA)

 READ (IN,1) ITMP

 1 FORMAT(I10)

 IF(ITMP.GE.0) GO TO 50

C

C2------IF ITMP LESS THAN ZERO REUSE DATA. PRINT MESSAGE AND RETURN.

 WRITE(IOUT,6)

 6 FORMAT(1H0,’REUSING MULTI-AQUIFER WELLS FROM LAST STRESS PERIOD’)

 GO TO 260

C

C3------ITMP=>0. SET NMAWS EQUAL TO ITMP.

 50 NMAWS=ITMP

 IF(NMAWS.LE.MXMAW) GO TO 100

C

C4------NMAWS > MXMAW. PRINT MESSAGE. STOP.

 WRITE(IOUT,99) NMAWS,MXMAW

 99 FORMAT(1H0,’ITMP(’,I4,’) IS GREATER THAN MXMAW(’,I4,’)’)

 STOP

C

C5------PRINT NUMBER OF MAW WELLS IN CURRENT STRESS PERIOD.

 100 WRITE (IOUT,2) NMAWS

 2 FORMAT(1H0,/

 + 1X,I5,’ MULTI-AQUIFER WELLS’)

C

C6------IF THERE ARE NO ACTIVE MAWS IN THIS STRESS PERIOD THEN RETURN

 IF(NMAWS.EQ.0) GO TO 260

C

C7------PRINT HEADING FOR MAW INPUT.

 WRITE(IOUT,3)

 3 FORMAT(1X,/

 + 1X,’ ROW COL M.A.W. RATE RADIUS RATIO WELL NO.’,/

 + 1X,’--------------------------------------------------------’)

 DO 250 II=1,NMAWS

C

C8------FOR EACH WELL READ AND PRINT ROW, COLUMN, RATE,

C8------# OF LAYERS SCREENED AND # IN CATEGORY.

READ (IN,4) I,J,Q,RRATIO,NUM,(LAYERS(L),L=1,10)

 4 FORMAT(2I10,2F10.0,I5,10I3)

WRITE (IOUT,7) I,J,Q,RRATIO,II

 7 FORMAT(2X,I8,I7,G16.5,F10.2,I8)

RMAW(1,II)=I

 RMAW(2,II)=J

RMAW(4,II)=RRATIO

C*** Q in Bennett et al (1982) is positive to for discharging well

C*** Q in this code is read as negative for discharging well

 RMAW(5,II)=-Q

 NWLAYS=0

 DO 230 L=1,10

 IF (LAYERS(L).EQ.0) GO TO 230

 NWLAYS=NWLAYS+1

 RMAL(1,NWLAYS,II)= L
70

 230 CONTINUE

 IF ((NWLAYS.ne.num).AND.(Q.ne.0)) THEN

 WRITE (IOUT,9) NUM,II

 9 FORMAT(5x,’NUMBER OF LAYERS NOT EQUAL’,I5,’ FOR WELL #’,I5)

 STOP

 ENDIF

 RMAW(3,II)=NWLAYS

 WRITE(IOUT,8) (RMAL(1,JJ,II), JJ= 1,NWLAYS)

 8 FORMAT(1X,’LAYERS: ’,10f5.0)

 DO 240 JJ=1,NWLAYS

C

C9------FOR EACH SCREENED LAYER READ LAYER # AND RADIUS RATIO.

 RMAL(2,JJ,II)=RRATIO

 240 CONTINUE

 250 CONTINUE

C

C10-----RETURN

 260 RETURN

 END

C

 SUBROUTINE MAW5FM(NMAWS,MXMAW,RHS,HCOF,IBOUND,RMAW,

 1 RMAL,NCOL,NROW,NLAY,CR,DELC,DELR,HNEW,IOUT)

C

C-----VERSION 5.0 OCT2001 MAW5FM

C

C ******************************************************************

C ADD MULTI-AQUIFER WELL FLOW TO RHS AND HCOF

C ******************************************************************

C

C SPECIFICATIONS:

C ------------------------------------------------------------------

 DOUBLE PRECISION HNEW

C

 DIMENSION RHS(NCOL,NROW,NLAY),HCOF(NCOL,NROW,NLAY),

 1 RMAW(5,MXMAW),

 2 RMAL(4,NLAY,MXMAW),CR(NCOL,NROW,NLAY),DELC(NROW),

 3 DELR(NCOL),HNEW(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY)

C ------------------------------------------------------------------

 PI=3.14159

C

C1------PROCESS EACH MAW CATEGORY

 DO 90 I1=1,NMAWS

C

C2------GET THE INFORMATION FOR THE CATEGORY

 NWLAYS=RMAW(3,I1)

 I=RMAW(1,I1)

 J=RMAW(2,I1)

 SUMF=0

 SUMFH=0

C

C3------PROCESS THE CELL IN EACH SCREENED LAYER

 DO 30 I2=1,NWLAYS

 K=RMAL(1,I2,I1)

C

C4------IF THE CELL IS NOT ACTIVE MOVE ON TO THE NEXT LAYER.

 IF(IBOUND(J,I,K).LE.0) GO TO 30

C

C5------CALCULATE THE TRANSMISSIVITY OF THE CELL.

 TRM=CR(J-1,I,K)*((DELR(J-1)+DELR(J))/2)/DELC(I)

 TRP=CR(J,I,K)*((DELR(J)+DELR(J+1))/2)/DELC(I)

 IF ((TRM.NE.0).AND.(TRP.NE.0)) GO TO 50

C

C6----IF EITHER TRANSMISSIVITY IS ZERO THEN CHECK TO SEE WHETHER

C WELL EXTENDS TO LOWER LAYERS, IF NOT, STOP SIMULATION

 IF (I2.LT.NWLAYS) GO TO 30
71

 WRITE(IOUT,888) K,I,J

 888 FORMAT(1H ,’MAW FAILS. WELL IN LAYER’,I5,’ ROW’,I5,’ COLUMN’,

 1 I5,’ IS ADJACENT TO AN INACTIVE CELL’)

 WRITE(IOUT,889)

 889 FORMAT(1H ,’SIMULATION ENDING’)

 STOP

C

C7------CALCULATE F AND FH, STORE THEM IN ARRAY RMAL AND ADD THEM

C7------TO ACCUMULATORS SUMF AND SUMFH.

 50 TR=2*TRP*TRM/(TRP+TRM)

 RRATIO=RMAW(4,I1)

 F=TR/(ALOG(RRATIO))

 SUMF=SUMF+F

 RMAL(3,I2,I1)=F

 HTMP=HNEW(J,I,K)

 FH=F*HTMP

 SUMFH=SUMFH+FH

 RMAL(4,I2,I1)=FH

30 CONTINUE

C*** check for dry wells

 IF (sumf.eq.0) THEN

 WRITE(IOUT,890) I, J

 890 FORMAT(’ DRY WELL IN ROW ’,I5, ’ COLUMN ’,I5)

 GO TO 90

 ENDIF

C

C8------CALCULATE THE HEADS IN THE WELLS IN THIS CATEGORY

 Q=RMAW(5,I1)

 HWELL=(SUMFH/SUMF)-(Q/(2*PI*SUMF))

C

C9------FOR EACH CELL ADD TERMS FOR THIS CATEGORY TO HCOF AND RHS.

C*** categories omitted in this version

 DO 60 I2=1,NWLAYS

 K=RMAL(1,I2,I1)

 IF(IBOUND(J,I,K).LE.0) GO TO 60

 F=RMAL(3,I2,I1)

NINCAT=1

 HCOF(J,I,K)=HCOF(J,I,K)-2*PI*F*NINCAT

 RHS(J,I,K)=RHS(J,I,K)-2*PI*F*HWELL*NINCAT

60 CONTINUE

90 CONTINUE

C

C10-----RETURN

 RETURN

 END

 SUBROUTINE MAW5BD(NMAWS,MXMAW,IBOUND,RMAW,

 1 RMAL,NCOL,NROW,NLAY,HNEW,KSTP,KPER,IMAWCB,

 2 ICBCFL,BUFF,IOUT,MSUM,DELT,VBNM,VBVL)

C-----VERSION 5.0 OCT2001 MAW5BD

C

C ******************************************************************

C CALCULATE CELL-BY-CELL FLOWS AND BUDGET TERMS FOR

C MULTI-AQUIFER WELLS

C ******************************************************************

C

C SPECIFICATIONS:

C ------------------------------------------------------------------

 DOUBLE PRECISION HNEW

C

 DIMENSION HNEW(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),

 1 RMAW(5,MXMAW),

 2 RMAL(4,NLAY,MXMAW),

 3 BUFF(NCOL,NROW,NLAY),VBVL(4,20),VBNM(4,20)

 DIMENSION TEXT(4)

 CHARACTER*4 TEXT,VBNM
72

 DATA TEXT(1),TEXT(2),TEXT(3),TEXT(4) /’MULT’,’I-AQ’,’IFR ’,’WELL’/

C ------------------------------------------------------------------

C

C1------CLEAR RATIN AND RATOUT ACCUMULATORS.

 IBD=0

 RATIN=0.

 RATOUT=0.

C

C2------IF THERE ARE NO MAWS DO NOT ACCUMULATE FLOW

 IF(NMAWS.EQ.0)GO TO 200

C

C3------TEST TO SEE IF CELL-BY-CELL FLOW TERMS WILL BE RECORDED.

 IF(ICBCFL.EQ.0 .OR. IMAWCB.LE.0 ) GO TO 10

C

C4------IF CELL-BY-CELL FLOWS WILL BE SAVED THEN CLEAR THE BUFFER.

 IBD=1

 DO 5 IL=1,NLAY

 DO 5 IR=1,NROW

 DO 5 IC=1,NCOL

 BUFF(IC,IR,IL)=0.

 5 CONTINUE

C

C5------PROCESS MAW CATEGORIES ONE AT A TIME.

 10 PI=3.14159

DO 150 I1=1, NMAWS

 NWLAYS=RMAW(3,I1)

 I=RMAW(1,I1)

 J=RMAW(2,I1)

 SUMF=0

 SUMFH=0

C

C5A-----FOR EACH CELL OPEN TO THE CATEGORY CALCULATE F AND FH

 DO 30 I2=1,NWLAYS

 K=RMAL(1,I2,I1)

 IF(IBOUND(J,I,K).LE.0) GO TO 30

 F=RMAL(3,I2,I1)

 SUMF=SUMF+F

 HTMP=HNEW(J,I,K)

 FH=F*HTMP

 SUMFH=SUMFH+FH

 RMAL(4,I2,I1)=FH

30 CONTINUE

C*** skip over if dry well

IF (sumf.eq.0) GO TO 150

C

C5B-----CALCULATE THE HEAD IN THE WELLS IN THIS CATEGORY

 Q=RMAW(5,I1)

 HWELL=(SUMFH/SUMF)-(Q/(2*PI*SUMF))

C

C5C-----FOR EACH LAYER IN WHICH THE MAW IS SCREENED PROCESS

C5C-----THE CELL WHICH CONTAINS THE MAW.

 DO 100 I2=1,NWLAYS

C

C5C1----CALCULATE RATE OF FLOW FROM THE MAWS INTO THE CELL.

 K=RMAL(1,I2,I1)

 IF(IBOUND(J,I,K).LE.0) GO TO 100

 F=RMAL(3,I2,I1)

 HTMP=HNEW(J,I,K)

NINCAT=1

 RATE=2*PI*F*(HWELL-HTMP)*NINCAT

C

C5C2----IF BUDGET TERMS ARE TO BE SAVED THEN ADD RATE TO BUFFER.

IF(IBD.EQ.1) BUFF(J,I,K)=BUFF(J,I,K)+RATE

C

C5C3----PRINT THE INDIVIDUAL RATES IF REQUESTED(IMAWCB<0).
73

 IF(IMAWCB.LT.0.AND.ICBCFL.NE.0) WRITE(IOUT,900) (TEXT(N),N=1,4),

 1 KPER,KSTP,I1,K,I,J,RATE,HWELL

 900 FORMAT(1H0,4A4,’ PERIOD’,I3,’ STEP’,I3,’ MAW’,I4,

 1 ’ LAYER’,I3,’ ROW ’,I4,’ COL’,I4,’ RATE’,G15.7,

 2 ’ HWELL’,G15.7)

 IF(RATE) 90,100,80

C

C5C4---- RATE IS POSITIVE(RECHARGE). ADD IT TO RATIN.

 80 RATIN=RATIN+RATE

 GO TO 100

C

C5C5----RATE IS NEGATIVE(DISCHARGE). ADD IT TO RATOUT.

 90 RATOUT=RATOUT-RATE

 100 CONTINUE

 150 CONTINUE

C

C6------IF CELL-BY-CELL TERMS WILL BE SAVED THEN CALL UBUDSV TO

C6-----RECORD THEM ON DISK

 IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,IMAWCB,BUFF,NCOL,NROW,

 1 NLAY,IOUT)

C

C7------MOVE RATES INTO VBVL FOR PRINTING BY MODULE BAS1OT.

 200 VBVL(3,MSUM)=RATIN

 VBVL(4,MSUM)=RATOUT

C

C8------MOVE RATES TIMES TIME STEP LENGTH INTO VBVL ACCUMULATORS.

 VBVL(1,MSUM)=VBVL(1,MSUM)+RATIN*DELT

 VBVL(2,MSUM)=VBVL(2,MSUM)+RATOUT*DELT

C

C9------MOVE BUDGET TERM LABELS INTO VBNM FOR PRINTING.

 VBNM(1,MSUM)=TEXT(1)

 VBNM(2,MSUM)=TEXT(2)

 VBNM(3,MSUM)=TEXT(3)

 VBNM(4,MSUM)=TEXT(4)

C

C10-----INCREMENT BUDGET TERM COUNTER(MSUM).

 MSUM=MSUM+1

C

C11-----RETURN

 RETURN

 END
C.E. Heywood and R.M. Yager—SIMULATED GROUND-WATER FLOW IN THE HUECO BOLSON, AN ALLUVIAL-BASIN AQUIFER SYSTEM

NEAR EL PASO, TEXAS—U.S. Geological Survey Water-Resources Investigations Report 02-4108

U.S. Department of the Interior

U.S. Geological Survey, WRD

5338 Montgomery Blvd. NE, Suite 400

Albuquerque, NM 87109-1311

BOOK RATE