GWM—A Ground-Water Management

Process for the U.S. Geological Survey

Modular Ground-Water Model

(MODFLOW-2000)

U.S. Department of the Interior

U.S. Geological Survey

Ground-Water Resources Program

Open-File Report 2005-1072

Ahlfeld and others—GW

M

—

A Ground-Water Management Process for the U.S. Geological Survey Modular Ground-Water Model (MODFLOW-2000)—OFR 2005-1072

COLUMNS

ROWS

LAY

E

R

S

1

2

3

1 2 3 4

1

2

3

4

5

6

Q4

Q3

Q2 Q1

IM1

RESERVOIR
GWM—A Ground-Water Management 

Process for the U.S. Geological Survey 

Modular Ground-Water Model 

(MODFLOW-2000)

By David P. Ahlfeld (University of Massachusetts), Paul M. Barlow, and

Ann E. Mulligan (Woods Hole Oceanographic Institution)

U.S. Department of the Interior

U.S. Geological Survey

Ground-Water Resources Program

Open-File Report 2005-1072
U.S. Department of the Interior

Gale A. Norton, Secretary

U.S. Geological Survey

Charles G. Groat, Director

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

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

Box 25286, Denver Federal Center

Denver, CO 80225

For more information about the USGS and its products:

Telephone: 1-888-ASK-USGS

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

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

endorsement by the U.S. Government.

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

reproduce any copyrighted materials contained within this report.

Suggested citation:

Ahlfeld, D.P., Barlow, P.M., and Mulligan, A.E., 2005, GWM—A ground-water management process for the U.S. 

Geological Survey modular ground-water model (MODFLOW-2000): U.S. Geological Survey Open-File Report 2005-1072, 

124 p.
 iii

Preface

This report describes a Ground-Water Management Process (GWM) for the U.S. Geological 

Survey modular three-dimensional ground-water model, MODFLOW-2000. The performance of 

the program has been tested in a variety of applications. Future applications, however, might 

reveal errors that were not detected in the test simulations. Users are requested to notify the U.S. 

Geological Survey of any errors found in this report or the computer program by using the address 

on the inside of the back cover of the report. Updates might occassionally be made to both the 

report and to the computer program. Users can check for updates on the Internet at URL 

http://water.usgs.gov/software/ground_water.html/.

 v

Contents

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

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

Numerical Modeling of Ground-Water Flow with MODFLOW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

Formulation of Ground-Water Management Problems with GWM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

Decision Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

Flow-Rate Decision Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

External Decision Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

Binary Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

Objective Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

Constraints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

Decision-Variable Constraints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

Linear-Summation Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

Hydraulic-Head Constraints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

Streamflow Constraints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

Complete Statement of Ground-Water Management Formulation Solved by GWM . . . . . . . . . . 18

Solution of Ground-Water Management Problems with GWM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Linear Formulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Calculation of Response Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

Accuracy and Precision of the Response Coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

Simplex Algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Nonlinear Formulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

Mixed-Binary Linear Formulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

Input Instructions and Output Files. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

Name File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

GWM Process Files. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

Description of Selected Conventions, Options, and Variables in GWM . . . . . . . . . . . . . . . . . 30

GWM File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

Decision variables (DECVAR) File. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

Objective Function (OBJFNC) File. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

Constraint Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

Decision-Variable Constraints (VARCON) File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

Linear-Summation Constraints (SUMCON) File. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

Head Constraints (HEDCON) File. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

Streamflow Constraints (STRMCON) File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

Solution and Output-Control Parameters (SOLN) File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

Output Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Sample Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Sample Problem 1: DEWATER. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Linear Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

Mixed-Binary Linear Formulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

Selected Input and Output Files for Sample Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
vi 

Sample Problem 2: SEAWATER. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

Nonlinear Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

Selected Input and Output Files for Sample Problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Sample Problem 3: SUPPLY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

Nonlinear Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

Selected Input and Output Files for Sample Problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

References Cited . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

Appendix 1. Examples of Decision Variables that can be Defined in GWM. . . . . . . . . . . . . . . . . . . . . . . . 99

Appendix 2. GWM Solution Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

Appendix 3. Programmers’ Guide to GWM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

Figures

1. Schematic graph showing seepage rate through a streambed as a function of 

ground-water head in a model cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2–5. Schematic illustrations showing:

2. Types of hydraulic-head constraints supported by GWM. . . . . . . . . . . . . . . . . . . . . . . . . . 16

3. Model grid for DEWATER sample problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4. Model grid for SEAWATER sample problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5. Model grid for SUPPLY sample problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

Conversion Factors and Abbreviations

GWF Ground-Water Flow

GWM Ground-Water Management

LP Linear programming

MPS Mathematical Programming System

RMS Response Matrix Solution

SLP Sequential linear programming

Multiply By To obtain

cubic foot (ft3) 0.02832 cubic meter

cubic foot per day (ft3/d) 0.02832 cubic meter per day

foot (ft) 0.3048 meter

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

square foot per day (ft2/d) 0.09290 square meter per day
GWM—A Ground-Water Management 

Process for the U.S. Geological Survey 

Modular Ground-Water Model 

(MODFLOW-2000)

By David P. Ahlfeld, Paul M. Barlow, and Ann E. Mulligan

Abstract

GWM is a Ground-Water Management Process for the U.S. Geological Survey modular three-dimensional ground-water model, MODFLOW-2000. GWM uses a responsematrix approach to solve several types of linear, nonlinear, and mixed-binary linear 

ground-water management formulations. Each management formulation consists of a set 

of decision variables, an objective function, and a set of constraints. Three types of decision 

variables are supported by GWM: flow-rate decision variables, which are withdrawal or 

injection rates at well sites; external decision variables, which are sources or sinks of water 

that are external to the flow model and do not directly affect the state variables of the simulated ground-water system (heads, streamflows, and so forth); and binary variables, which 

have values of 0 or 1 and are used to define the status of flow-rate or external decision variables. Flow-rate decision variables can represent wells that extend over one or more model 

cells and be active during one or more model stress periods; external variables also can be 

active during one or more stress periods. A single objective function is supported by GWM, 

which can be specified to either minimize or maximize the weighted sum of the three types 

of decision variables. Four types of constraints can be specified in a GWM formulation: 

upper and lower bounds on the flow-rate and external decision variables; linear summations of the three types of decision variables; hydraulic-head based constraints, including 

drawdowns, head differences, and head gradients; and streamflow and streamflowdepletion constraints.

The Response Matrix Solution (RMS) Package of GWM uses the Ground-Water Flow 

Process of MODFLOW to calculate the change in head at each constraint location that 

results from a perturbation of a flow-rate variable; these changes are used to calculate the 

response coefficients. For linear management formulations, the resulting matrix of 

response coefficients is then combined with other components of the linear management 

formulation to form a complete linear formulation; the formulation is then solved by use of 

the simplex algorithm, which is incorporated into the RMS Package. Nonlinear formulations arise for simulated conditions that include water-table (unconfined) aquifers or headdependent boundary conditions (such as streams, drains, or evapotranspiration from the 

water table). Nonlinear formulations are solved by sequential linear programming; that is, 

repeated linearization of the nonlinear features of the management problem. In this 

approach, response coefficients are recalculated for each iteration of the solution process. 

Mixed-binary linear (or mildly nonlinear) formulations are solved by use of the branch and 

bound algorithm, which is also incorporated into the RMS Package.
2 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Three sample problems are provided to demonstrate the use of GWM for typical 

ground-water flow management problems. These sample problems provide examples of 

how GWM input files are constructed to specify the decision variables, objective function, 

constraints, and solution process for a GWM run. The GWM Process runs with the MODFLOW-2000 Global and Ground-Water Flow Processes, but in its current form GWM cannot be used with the Observation, Sensitivity, Parameter-Estimation, or Ground-Water 

Transport Processes. The GWM Process is written with a modular structure so that new 

objective functions, constraint types, and solution algorithms can be added.

Introduction

Since the 1960s, numerical ground-water flow models have become increasingly 

important tools for the analysis of ground-water systems. More recently, ground-water 

flow models have been combined with optimization techniques to determine waterresource management strategies that best meet a particular set of management objectives 

and constraints. Optimization techniques are a set of mathematical programs that seek to 

find the optimal (or best) allocation of resources to competing uses. In the context of 

ground-water management, the resources are typically the ground- and surface-water 

resources of a basin and (or) the financial resources of the communities that depend on the 

water. The management objectives and constraints are stated (or formulated) mathematically in an optimization (management) model. Combined ground-water flow and optimization models have been applied to various ground-water management problems, including 

the control of water-level declines and land subsidence that could result from ground-water 

withdrawals, conjunctive management of ground-water and surface-water systems, capture 

and containment of contaminant plumes, and seawater intrusion. As applied in U.S. Geological Survey (USGS) studies, management agencies and other stakeholders provide 

information on water-resource objectives and constraints. The USGS then provides scientific data, analysis, and expertise in ground-water flow and optimization modeling to help 

decisionmakers understand how the characteristics of the hydrologic (ground-water and 

surface-water) system and the stated objectives and constraints interact to affect options for 

managing the resource. 

A number of computer codes have been developed during the past two decades to 

facilitate linked flow and optimization modeling of ground-water flow systems (Lefkoff 

and Gorelick, 1987; Greenwald, 1998; Zheng and Wang, 2002; Ahlfeld and Riefler, 2003; 

Peralta, 2004). These codes differ in the numerical model used to represent the groundwater flow system, the types of ground-water management problems that can be solved, 

and the approaches used to solve the management problems.

GWM is a new process for the USGS MODFLOW-2000 modular ground-water model 

(Harbaugh and others, 2000). The GWM Process solves several types of linear, nonlinear, 

and mixed-binary linear ground-water management problems. The response-matrix 

approach, which has been used widely in ground-water management modeling, is used to 

transform a ground-water management problem into an optimization formulation that can 

be solved by GWM. GWM uses the simplex and branch and bound optimization algorithms 

to solve the resulting formulations; these algorithms have been coded internally in GWM 

in the FORTRAN-90 computer language. In its current form, the GWM Process can only 

be used with the MODFLOW-2000 Ground-Water Flow (GWF) and Global Processes; it 

cannot be used with the Ground-Water Transport, Observation, Sensitivity, or ParameterEstimation Processes. Currently (2005), MODFLOW-2000 (Harbaugh and others, 2000) is 

the most recent version of the MODFLOW code, which was originally developed in the 
Numerical Modeling of Ground-Water Flow with MODFLOW 3

1980s (McDonald and Harbaugh, 1988; Harbaugh and McDonald, 1996). Unless otherwise 

noted, in the remainder of this report, the term MODFLOW will refer to the MODFLOW2000 version of the code.

The origin of GWM is the MODOFC code developed by Ahlfeld and Riefler (2003). 

MODOFC is based on the MODFLOW-96 version of MODFLOW (Harbaugh and 

McDonald, 1996). Several modifications have been made in the transition from MODOFC 

to GWM. First, the types of decision variables and constraints that can be specified in a 

management-model formulation have been modified and expanded. For example, it is now 

possible to include decision variables that represent a source or sink of water that is external 

to the model domain. Second, the structure of the input and output files has been modified 

substantially. In the new structure, separate input files are used to define the decision variables, the objective function, the constraints, and the solution technique. The new structure 

is intended to facilitate the process of converting each management formulation into a format that can be solved by GWM; the revised structure also should make the process of 

developing new options and associated computer modules for GWM easier. Finally, many 

changes were made to the original MODOFC FORTRAN code; these changes included 

organizing the subroutines into modules and packages that are consistent with the overall 

MODFLOW-2000 structure.

This report describes the formulation of ground-water management problems that can 

be solved with GWM, the approaches that GWM uses for solving the management problems, and the input and output files associated with a GWM run. The report also includes 

three sample problems of the application of GWM to typical ground-water management 

problems. These sample problems provide examples of the input files needed by GWM and 

the output that is generated by GWM. This report, however, is not a guide to the application 

of optimization modeling for aquifer management. For detailed guides to the application of 

management models, refer to textbooks by Willis and Yeh (1987), Gorelick and others 

(1993), and Ahlfeld and Mulligan (2000), and to literature reviews by Gorelick (1983), Yeh 

(1992), Ahlfeld and Heidari (1994), and Wagner (1995). 

The report begins with a brief overview of ground-water flow modeling using the 

MODFLOW Ground-Water Flow Process. The purpose of the overview is to introduce 

concepts that are important to the GWM Process. For more details on the theory and use of 

the Ground-Water Flow Process, the reader is referred to Harbaugh and others (2000).

Numerical Modeling of Ground-Water Flow with MODFLOW

The partial-differential equation of ground-water flow used in MODFLOW is 

(McDonald and Harbaugh, 1988)

 , (1)

where

is the potentiometric head (L);

are values of hydraulic conductivity along the , , and coordinate 

axes, which are assumed to be parallel to the major axes of hydraulic 

conductivity (L/T);

is a volumetric flow rate per unit volume, and represents sources 

and/or sinks of water (T-1);

is the specific storage of the porous material (L-1); and

is time (T).

∂x

∂ Kxx

∂h

∂x

------ ⎝ ⎠ ⎛ ⎞ ∂y

∂ Kyy

∂h

∂y

------ ⎝ ⎠ ⎛ ⎞ ∂

∂z

---- Kzz

∂h

∂z

----- ⎝ ⎠ ⎛ ⎞ W Ss

∂h

∂t

+ + – = -----

h

Kxx Kyy Kzz , , x y z

W

Ss

t
4 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Equation 1 describes three-dimensional ground-water flow of constant density under 

nonequilibrium conditions in a heterogeneous and anisotropic medium, provided the principal axes of hydraulic conductivity are aligned with the coordinate directions (McDonald 

and Harbaugh, 1988). The ground-water flow equation, when combined with a particular 

set of boundary and initial conditions, constitutes a mathematical representation of a 

ground-water flow system. A solution to equation 1 gives the distribution of head (the state 

variable) throughout the flow system as a function of space and time, .

There are three main types of boundary conditions that can be specified for a given 

location and time along the boundary of the model domain. These are (1) specified-head 

boundaries, (2) specified-flow boundaries, and (3) head-dependent flow boundaries. A 

fourth type of boundary condition is required if a free surface, or water table, is present. 

Because the water table is a free-surface boundary whose position is not fixed, the location 

of the boundary is not known. This variable boundary results in a nonlinear relation 

between the position of the water table and conditions along the boundary of, or stresses to, 

a ground-water system. The presence of free-surface or head-dependent boundary conditions has implications for the solution of the ground-water management problem, as 

discussed later in the report.

MODFLOW uses the finite-difference method to solve the partial-differential equation of ground-water flow. In this method, the continuous system described by equation 1 

is converted into a finite set of simultaneous linear algebraic difference equations that 

define the state variable at discrete node points and times. The basic concept of 

the finite-difference method is to approximate the derivatives in equation 1 using the difference in head at adjoining nodes during specified time intervals. To define the locations 

of the node points, the ground-water flow system is discretized into a mesh of blocks 

(called cells) that consists of a set of rows, columns, and layers. Each cell within the mesh 

is identified by its row, column, and layer position by the indexing scheme (that is, 

row , column , and layer ). MODFLOW’s block-centered formulation of the finitedifference equations places the node points at the center of the cells.

A finite-difference approximation also is defined for the time derivative of head on the 

right-hand side of equation 1, in which the continuous time derivative is replaced by a set 

of discretized time intervals called time steps. MODFLOW uses a backward-difference 

approximation to the time derivative because that approach is always numerically 

stable—that is, errors introduced at any time diminish progressively at succeeding times 

(McDonald and Harbaugh, 1988). MODFLOW defines stress periods as time intervals 

during which all of the external stresses to the simulated ground-water flow system are constant. Stress periods are divided into time steps to achieve adequate numerical accuracy. 

Stress periods defined for a particular MODFLOW simulation can be of variable lengths; 

ground-water management problems solved by GWM also can use stress periods that are 

of variable lengths.

Solution of the finite-difference equations gives an approximate value of head, which 

will be notated as , at each of the discrete cells and for each of the discrete times. For 

steady-state problems, the subscript can be ignored and the head is understood to be that 

at steady state. Most of the methods available in MODFLOW to solve the set of finitedifference equations are iterative methods, wherein the equations are solved repeatedly 

(iteratively) for each time step beginning with an initial estimate of the head distribution 

throughout the model domain and stopping when a head distribution is obtained that meets 

one or more convergence (or closure) criteria specified by the user. These convergence criteria most commonly include a requirement that the changes in calculated heads from one 

hx y z t ( ) ,,,

hx y z t ( ) ,,,

ijk , ,

i j k

hijkt ,, ,

t
Numerical Modeling of Ground-Water Flow with MODFLOW 5

iteration to the next must be less than a specified value. Iterative procedures yield only 

approximate solutions to the finite-difference equations for each time step; the accuracy 

of these solutions depends on many factors, including the closure criteria specified by the 

user. Because of the importance of the accuracy of MODFLOW calculations to the solutions of ground-water management problems solved by the GWM Process, the issue of 

model accuracy will be discussed in more detail in the section “Solution of Ground-Water 

Management Problems with GWM.”

The source/sink term, , in equation 1 is used to represent stresses imposed on a 

ground-water flow system. In MODFLOW, these stresses are flow rates applied to the 

simulated ground-water flow system at specified locations and stress periods, and, after 

multiplication by appropriate volumes, have dimensions of volume per unit time (L3/T). 

One type of flow stress that can be represented by is a withdrawal (discharge) or injection (recharge) rate at a simulated well site, , where represents both the location of 

the nth well site and the stress period (or periods) during which the well operates. As used 

in MODFLOW, is positive ( ) for flow into the ground-water system at the 

well (injection) and negative ( ) for flow out of the ground-water system at the well 

(withdrawal).

Withdrawal and injection of water at simulated wells affect the head distribution of the 

simulated ground-water flow system. At a particular location and time in the modeled flow 

system, the relation between head and simulated stresses at the wells can be expressed 

mathematically as

 , (2)

where represents the vector (set) of all withdrawal and injection rates at all well locations and all operative stress periods. Equation 2 simply states that there is a functional relation between the state variable calculated by MODFLOW and the simulated wells. 

Equations similar to equation 2 also could be written for other head-based state variables 

indirectly calculated by MODFLOW, such as head drawdowns and gradients.

Ground-water flow systems commonly interact with streams and other surface-water 

features that are in hydraulic connection with the underlying ground-water system. For 

the case of streams, this interaction occurs as seepage of water across the ground-water/

streambed interface. The rate of seepage to or from a particular stream reach at a particular 

time is dependent upon the conductance of the streambed and the head gradient between 

the ground-water flow system and the stream; when ground-water head is above the stream 

bottom, the seepage rate is described by

 , (3)

where

is the flow rate (seepage rate) between the ground-water flow system at cell

 and stream reach at time (L3/T);

is the streambed conductance in stream reach (L2/T);

is the head in stream reach at time (L); and

is the ground-water head in cell at time (L).

Prudic (1989) developed a Streamflow-Routing (STR1) Package for MODFLOW 

that simulates hydraulic interactions between the simulated ground-water system and 

adjoining streams, and also keeps track of the amount of water within each simulated 

W

W

Qwn n

Qwn Qwn 0.0 >

Qwn 0.0 <

hijkt ,, , h = ijkt ,, , ( ) Qw

Qw

hijkt ,, ,

Qsbr t , Csbr Hsr t , hijkt ,, , = ( ) –

Qsbr t ,

ijk , , r t

Csbr r

Hsr t , r t

hijkt ,, , ijk , , t
6 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

stream reach, , where represents both the location of the rth stream reach and the 

stress period of flow in the reach. Each stream reach corresponds to an individual cell in 

the finite-difference grid used to simulate ground-water flow. The flow in each stream 

reach consists of inflows from adjacent upstream reaches, diversions within the reach, and 

the gain or loss of streambed seepage.

The functional relation between ground-water head in a model cell and the seepage 

rate between the aquifer and stream reach is shown schematically in figure 1. An important 

characteristic of the relation is that stream seepage is a nonlinear function of head. If head 

in the aquifer is above the bottom of the streambed (SBOT in fig. 1), then there is a linear 

relation between head in the model cell and seepage to or from the stream. If head in the 

aquifer falls below the bottom of the streambed, however, then the ground-water head is set 

to the elevation of the bottom of the streambed, and the seepage rate from the stream to the 

model cell depends only on stream stage (as long as there is flow in the stream reach to meet 

the seepage rate).

Qsfr r

Slope = -Csbr

Qsbr,t

Positive Qsbr,t

indicates flow

into aquifer

Negative Qsbr,t

indicates flow

into stream

0

SBOT Hsr,t

hi,j,k,t

Figure 1. Seepage rate, , through a streambed as a function of ground-water head, , in a 

model cell, where is the streambed bottom, is streambed conductance, and is head in 

the stream. Figure modified from Prudic (1989).

Qsbr t, hijkt ,, ,

SBOT Csbr Hsr t ,
Formulation of Ground-Water Management Problems with GWM 7

Because streamflow calculated by the STR1 Package is a function of seepage between 

each simulated stream reach and the modeled ground-water flow system, and seepage, in 

turn, is a function of head in the ground-water system (eq. 3), calculated streamflows are 

also a function of the head in the aquifer. Therefore, streamflows calculated by the STR1 

Package can be considered to be state variables of the flow model, and, as with heads, a 

functional relation between streamflows and simulated wells can be defined:

. (4)

Prudic and others (2004) developed a second Streamflow-Routing (SFR1) Package that can 

be used instead of the original STR1 Package. Because the SFR1 Package had not been 

released during the development of GWM, the functional relations between streamflow and 

simulated wells are calculated in GWM by use of the STR1 Package.

A calibrated MODFLOW ground-water flow model provides the basis by which 

GWM develops specific functional relations between the simulated wells and calculated 

heads and streamflows described in general form by equations 2 and 4. These relations are 

referred to as response coefficients (or response functions), and are described in more detail 

in the solution-procedures section of the report.

Formulation of Ground-Water Management 

Problems with GWM

This section describes the components available in the GWM Process for formulating 

ground-water management problems. A ground-water management formulation consists 

of three components: decision variables, an objective function, and a set of constraints. 

Together, these three components define a mathematical model of the management decision-making (or design) process (Ahlfeld and Mulligan, 2000; Hillier and Lieberman, 

2001).

The decision variables of the management problem are the quantifiable controls (or 

decisions) that are to be determined by the model, such as the withdrawal rates at a set of 

managed wells. The values determined by GWM for these control variables define the solution of the problem. The objective function of the problem, which is stated in terms of one 

or more of the decision variables, is a measure of the performance of the design process; 

the objective function is used to identify the best solution among many possible solutions. 

This function may be maximized or minimized, depending upon the GWM application. The 

third component of the management problem is a set of constraints that impose restrictions 

on the values that can be taken by the decision variables. The solution of a well-defined 

ground-water management formulation consists of values for the decision variables that 

optimize the objective function while satisfying all constraints on decision-variable values 

(Ahlfeld and Mulligan, 2000).

Decision Variables

GWM supports three types of decision variables: flow-rate decision variables, external 

decision variables, and binary variables.

Qsfr Qsf = r ( ) Qw
8 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Flow-Rate Decision Variables

The primary type of decision variable is a withdrawal (discharge) or injection 

(recharge) flow rate at a managed well site, . All flow-rate decision variables are 

treated as positive values in GWM, whether they represent withdrawal or injection. For 

flow-rate variables that represent withdrawal, GWM internally converts to the correct sign 

for consistency with the sign convention of the MODFLOW GWF Process. Each flow-rate 

decision variable can extend over one or more model cells and can be active during one or 

more stress periods; the stress periods for which a decision variable are defined do not need 

to be of equal length. In the simplest case, is the withdrawal or injection rate at a single 

model cell during a single stress period. In more complicated problems, might consist 

of a set of vertical cells that simulate withdrawal from a well screened over several model 

layers, or a set of cells within a single model layer that represent the aggregated withdrawal 

from several wells in a subarea of the model domain. Although each decision variable can 

be defined to extend over multiple cells and multiple stress periods, the flow rate determined for the decision variable is a single value that is constant for each stress period 

during which the decision variable is active, and is apportioned to each cell within the areal 

extent of the decision variable on a percentage basis specified by the user. 

Each flow-rate decision variable represents either withdrawal or injection at the well 

site; a single decision variable cannot be used for both withdrawal and injection. However, 

the user can specify two flow-rate decision variables for a well site that are defined by the 

same set of model cells and stress periods, with one of the decision variables defined as a 

withdrawal well and the other defined as an injection well. Moreover, for transient models, 

a single well site may have more than one flow-rate decision variable associated with it, 

wherein each decision variable for the site has the same location but is active during different stress periods. This is a common situation for wells at which withdrawal rates vary with 

time.

It is important to recognize the difference between flows simulated by the Well (WEL) 

or Multi-Node Well (MNW) Packages of the GWF Process and those defined as decision 

variables for a GWM problem. For the purpose of a GWM problem, flows simulated by the 

WEL (Harbaugh and others, 2000) or MNW (Halford and Hanson, 2002) Packages have 

user-defined withdrawal or injection rates, and are referred to as unmanaged wells. These 

unmanaged withdrawal and injection rates are considered to be background stresses, in the 

sense that they will contribute to the total stress of the ground-water flow system in the 

absence of managed withdrawals and injections. In contrast, for a flow-rate decision variable defined in a GWM problem, the withdrawal or injection rate is unknown at the start of 

the GWM problem, and is determined as part of the solution process. Flow rates defined as 

decision variables are referred to as managed flows, because they are part of the management solution. At the start of a GWM problem, the managed flows are candidates for 

possible selection into the final set of active (that is, nonzero) flow-rate decision variables 

that compose the solution of the problem. It is possible that some of the candidate decision 

variables may not be selected as part of the solution to the problem, in which case the decision variables will be inactive (calculated withdrawal or injection flow rates of zero).

GWM allows for the simultaneous use of both managed and unmanaged wells at 

model cells. For example, the user might specify an unmanaged withdrawal rate (that is, a 

background stress) of 1.0 ft3/s at a particular cell with the WEL Package; the user also could 

define a managed withdrawal at the same cell by use of a flow-rate decision variable in 

GWM. The total withdrawal rate at the cell at the end of the GWM run would then equal 

the sum of the unmanaged withdrawal rate (1.0 ft3/s) and the managed withdrawal rate 

determined by GWM for the decision variable.

Qwn

Qwn

Qwn
Formulation of Ground-Water Management Problems with GWM 9

Although flow-rate decision variables have been described with reference to a well site 

(and are simulated in the GWM Process in a manner that is identical to a well simulated by 

the MODFLOW WEL Package), these variables could represent other types of managed 

stresses to a ground-water system, such as the recharge rate to an artificial recharge basin 

or the discharge rate to an irrigation drain. Several examples of how flow-rate decision variables can be defined in a GWM problem are provided in Appendix 1 of this report. Those 

examples also can be used as a guide during preparation of the input files for GWM.

External Decision Variables

The second type of decision variable that can be specified in GWM is a source or sink 

of water that does not have a direct effect on the state variables of the ground-water flow 

system. All external variables are treated as positive values in GWM, whether they represent a source or sink of water. These decision variables are external to the ground-water 

flow system, and therefore, are called external variables, , where represents the mth 

external variable. In contrast to the flow-rate decision variables, external variables are not 

part of the MODFLOW GWF Process, and response coefficients between the variables and 

MODFLOW state variables are not determined by GWM. Sources of external water to the 

management model are referred to as imports, whereas sinks of water are referred to as 

exports. An example of a source of water to the management model is an interbasin transfer 

from a surface-water reservoir that is external from the ground-water basin; additional 

examples of sources of water are treated wastewater available for artificial recharge or 

desalinated water from a desalination plant. An example of an export of water is a fraction 

of the ground water withdrawn from the simulated ground-water basin that is exported out 

of the basin and, therefore, is unavailable to meet within-basin water-supply demands.

Binary Variables

The third type of decision variable supported by GWM is a binary variable, which is 

defined to indicate the status of associated sets of flow-rate and external decision variables. 

One or more flow-rate or external decision variables, or combinations of flow-rate and 

external decision variables, can be associated with a single binary variable (see example 3 

in Appendix 1). Binary variables are notated as (where represents the lth binary variable), and have values of 0 or 1 (written ). If equals 1, at least one of the flow-rate 

or external decision variables associated with the binary variable is active (that is, the site 

has been constructed or is operational); if equals 0, all of the associated flow-rate and 

external variables are inactive. The binary variables are not directly associated with specific 

stress periods of the model; rather, the stress periods are implied by those assigned to the 

associated flow-rate and external decision variables.

Binary variables often are used to model the construction costs of a water-supply facility (a well site or external facility), but can have other uses in a management model. For 

example, there may be a requirement that the total number of active well sites constructed 

by a water-supply agency must be less than or equal to the total number of developable sites 

owned by the agency; alternatively, the water agency might require that the withdrawal rate 

at each well site either be zero (that is, the site is not constructed) or, if the site is constructed, greater than some minimum, nonzero withdrawal rate.

Binary variables can impose a significant computational burden on the GWM solution 

process, and they should be used with caution when nonlinear responses are present in the 

ground-water model.

Exm m

Il l

Il 0 1 = , Il

Il
10 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Objective Function

GWM supports a single objective function, which is to minimize or maximize the 

weighted sum of the three types of decision variables:

 , (5)

where

is the cost or benefit per unit volume of water withdrawn or injected at well 

site ;

is the cost or benefit per unit volume of water imported or exported at 

external site ;

is the unit cost or benefit associated with the binary variable ;

is the total duration of flow at well site ;

is the total duration of flow at external site ; and

are the total number of flow-rate, external, and binary decision variables, 

respectively.

 and are calculated by GWM by summing the duration of all stress periods during 

which the nth or mth decision variable is active. (Again, note that GWM does not require 

that stress periods specified in a MODFLOW simulation be of equal length.) The coefficients , , and are called the objective-function coefficients.

Equation 5 can be considered in terms of economic costs or benefits of the water withdrawn, injected, imported, or exported in the management model, so that each of the three 

terms of equation 5 has monetary units (such as dollars). The first term in equation 5 is the 

cost (or benefit) of withdrawing or injecting water at the flow-rate decision variables, the 

second term is the cost (or benefit) of importing or exporting water at the external variables, 

and the third term is the cost (or benefit) of making the flow-rate or external variables 

associated with each binary variable active. Note that the coefficients in the first two terms 

in the equation ( and ) imply that the cost (or benefit) of water is linearly proportional 

to the volume of water withdrawn, injected, imported, or exported.

Although equation 5 can be considered in monetary terms, it is often difficult to assign 

economic costs and benefits to all components of a water-resource management problem. 

In cases where economic costs and benefits are unknown or are not required in the management formulation, the user can specify the and coefficients as relative costs or 

benefits among the different flow-rate or external variables. For example, the objective 

might be to maximize the water withdrawn at well sites, in which case the coefficients 

could be set to dimensionless values of 1.0. Equation 5 then would become

Maximize , (6)

and the resulting objective value is in terms of total volume of water withdrawn. Alternatively, the user might want to weight the withdrawal from one well site as being twice as 

beneficial as the withdrawal from a second site, in which case the objective function would 

simplify to

Maximize . (7)

βn QwnTQwn 

n 1=

N

∑ γ m ExmTExm

 

m 1=

M

∑ κl Il

l 1=

L

+ + ∑

βn

n

γ m

m

κl Il

TQwn n

TExm m

NML , ,

TQwn TExm

βn γ m κl

βn γ m

βn γ m

n βn

1.0QwnTQwn

n 1=

N

∑

2.0Qw1TQw1 1.0Qw2TQw2 +
Formulation of Ground-Water Management Problems with GWM 11

Regardless of the meaning assigned to the coefficients , , and in equation 5, the user 

must ensure that the units of these coefficients are consistent. When GWM solves the management problem, it is the relative magnitude of each term in equation 5 that will determine 

the solution obtained.

Constraints

GWM supports four general types of management-model constraints that are described 

in the four subsections that follow. These constraints can be divided broadly into two types: 

those for which response coefficients need not be generated (constraints on the decision 

variables themselves and linear-summation constraints), and those for which response 

coefficients between the decision variables and ground-water flow system state variables 

must be generated (the hydraulic-head and streamflow constraints). Constraints that do not 

require generation of response coefficients are described first.

Decision-Variable Constraints

Lower and upper bounds are commonly placed on the flow-rate and external decision 

variables of a ground-water management model to express bounds on the yield of a well, 

injection rate of a well, and minimum and maximum capacities of an external source or sink 

of water. These lower and upper bounds are written mathematically as

(8)

and

, (9)

where and are the lower and upper bounds on flow-rate decision variable 

, respectively, and and are the lower and upper bounds on external variable 

, respectively. GWM requires that , , , and all be greater than or equal 

to zero, and that and .

If a flow-rate or external decision variable has not been associated with a binary variable in GWM, the lower bound for each variable ( or ) must be defined as 0. 

Equations 8 and 9 then become

(10)

and

. (11)

The user can specify nonzero lower bounds for decision variables that are not associated 

with a binary variable by use of linear-summation constraints described in the next section 

of the report.

βn γ m κl

Qwn

l

 Qwn Qwn

u ≤ ≤

Exm

l

 Exm Exm

u ≤ ≤

Qwn

l Qwnu

n Exm

l Exmu

m Qwn

l Qwnu Exml Exmu

Qwn

l

 Qwn

u ≤ Exml

 Exm

u ≤

Qwn

l Exml

0 Qwn Qwn

u ≤ ≤

0 Exm Exm

u ≤ ≤
12 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

If a flow-rate or external decision variable has been associated with a binary variable, 

then enhanced lower-bound constraints can be placed on the variable. Specifically, each 

variable can be required to be inactive, with a volumetric flow rate of zero, or active, with 

a nonzero lower bound ( or ). This either/or type of constraint can be written with 

binary variables as

(12)

and

, (13)

where is the binary variable associated with decision variables and (or) . Note 

that if , then constraints 12 and 13 imply that the variable lower bound is 0, whereas 

if , then the lower bound is the specified nonzero bound value. Constraints 12 and 13 

are applied separately to each decision variable associated with a given binary variable. 

Additional constraints are automatically added by GWM to force the decision variable to 

be 0 if the binary variable is 0 and to allow the decision variable to take any value in the 

range described by equations 8 and 9 if the binary variable is 1.

Linear-Summation Constraints

Constraints are available in GWM for linear summation of the three types of decision 

variables. These constraints have the general form

, (14)

, (15)

and

, (16)

where and are specified coefficients, is any of the three types of decision variables ( , and ), and is the total number of terms in the summation. Because all 

flow-rate and external decision variables are treated in GWM as positive values whether 

they are withdrawals, injections, imports, or exports, the user must ensure that the sign of 

each coefficient in equations 14–16 is defined to achieve the desired constraint form. The 

user also must ensure that the units used throughout a GWM formulation, including the 

coefficients in equations 14–16, are consistent.

When specifying a linear-summation constraint in GWM, it is not necessary to specify 

the stress period (or periods) for which the constraint is active, because the time period(s) 

for which each constraint applies is implicitly defined by the stress period(s) for which the 

flow-rate or external decision variables have been defined. The generalized nature of these 

constraints requires that the user ensure that the constraint definition makes sense. For 

example, if a withdrawal variable has been defined for a well site as in stress period 1 

Qwn

l Exml

Qwn Qwn

l

I ≥ l

Exm Exm

l

I ≥ l

Il Qwn Exm

Il 0 =

Il 1 =

apGVp b ≤ p

p 1 =

P

∑

apGVp b ≥ p

p 1 =

P

∑

apGVp 

p 1 =

P

∑ b = p

ap bp GVp

Qwn Exm Il P

Qw1
Formulation of Ground-Water Management Problems with GWM 13

and a second withdrawal variable has been defined for the well site as for stress period 

2, there may be a requirement that withdrawal at the site not decrease by more than 25 percent from the first to the second stress period, that is, . This constraint is 

specified in GWM as

, (17)

with the stress periods for which the constraint applies implied by the stress periods 

associated with the two decision variables.

Summation constraints described by equations 14–16 can be used in a variety of ways. 

A few examples of the use of the summation constraints are described in the next paragraphs.

Water-supply demands: Water-management models often include constraints that 

describe the demands on the water-supply system. These constraints can be written in a 

general way as

, (18)

where and are the coefficients (weights) associated with withdrawals and 

sources , respectively, and are the total number of withdrawal and source locations, respectively, and is the water-supply demand during stress period .

As an example of the use of equation 18, an irrigation district might require that the 

total ground water withdrawn from its three wells, in addition to the amount that is available 

from an out-of-basin source (an import), must meet the district’s total irrigation-season 

water-supply demands. If the three wells from which the water is withdrawn during the irrigation season are designated , , and , all sources of water are assigned equal 

weights (that is, all coefficients equal 1.0), the amount available from the out-of-basin 

source is designated as , and the total demand during the irrigation season is 

designated , then the constraint can be written as

.

Net stress constraints: An extension of equation 18 is to require that the difference 

between the total withdrawal and total injection rates from a set of candidate withdrawal 

and injection wells be constrained between upper and lower limits. To impose such a constraint, equations 14 and 15 can be written as

(19)

and

, (20)

where and are the candidate withdrawal and injection wells, and are the 

total numbers of candidate withdrawal and injection wells, and and are the upper 

and lower bounds on the net stress during the stress period (or periods), respectively.

Qw2

Qw2 0.75 ≥ Qw1

 0.75Qw1 – 1.0Qw2 + ≥ 0

aQwnQwn aExmExm 

m 1 =

M

+ ∑ ≥ Dt

n 1 =

N

∑

aQwn aExm Qwn

Exm N M

Dt t

Q1 Q2 Q3

IM1

D

Q1 Q2 Q3 IM1 +++ ≥ D

Qwnw 

nw 1 =

NW

∑ Qwni Qwu ≤

ni 1 =

NI

– ∑

Qwnw 

nw 1 =

NW

∑ Qwni Qwl ≥

ni 1 =

NI

– ∑

Qwnw Qwni NW NI

Qwu Qwl
14 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Stress ratio constraints: Instead of constraining the total net stress from a set of 

candidate withdrawal and injection wells, as in the previous example, there may be a desire 

to define the ratio between total withdrawals from and total injections to a set of candidate 

wells. These constraints can be written as

(21)

and

, (22)

where and are specified stress ratios and other terms are defined as they are for equations 19 and 20. Equations 21 and 22 are rearranged for use in GWM, resulting in

(23)

and

. (24)

Equations 23 and 24 can be used, for example, in the design of a capture and containment system of a ground-water contamination plume, in which a fraction of the water withdrawn to contain the plume is chemically treated and then reinjected into the ground-water 

system.

Constraints on the total number of active well sites: Upper and lower bounds on 

the total number of active (that is, nonzero) well sites can be specified by use of binary variables and summation constraints. These constraints are written in the general form

(25)

and

, (26)

where and are the upper and lower bounds on the number of active well sites and 

 is the total number of candidate well sites. For example, a water district might have 

access to a total of eight possible sites at which wells could be constructed, but funding is 

available for a maximum of five wells only. In addition, a minimum of three wells is desired 

by the water district to guard against the risk of one or two of the wells being taken out of 

service because of damage to the wells or contamination of the ground water. These 

requirements could be specified as

(27)

and

. (28)

Qwnw

nw 1 =

NW

∑

Qwni

ni 1 =

NI

∑

----------------------- ≥ a

Qwnw

nw 1 =

NW

∑

Qwni

ni 1 =

NI

∑

----------------------- 1

b ≤ --

a b

Qwnw 

nw 1 =

NW

∑ a Qwni 

ni 1 =

NI

– ∑ ≥ 0.0

Qwni

ni 1 =

NI

∑ b Qwnw 

nw 1 =

NW

– ∑ ≥ 0.0

Il 

l 1 =

L

∑ N U ≤

Il 

l 1 =

L

∑ NL ≥

N U N L

L

I1 I2 I3 I4 I5 I6 I7 I8 +++++++ ≤ 5

I1 I2 I3 I4 I5 I6 I7 I8 +++++++ ≥ 3
Formulation of Ground-Water Management Problems with GWM 15

Hydraulic-Head Constraints

GWM supports four types of hydraulic-head constraints. The first type of constraint is 

an absolute lower and (or) upper bound placed on a head at a specific location and stress 

period:

(29)

and

, (30)

where and are the specified lower and upper bounds on head at location 

(fig. 2A) at the end of stress period (that is, at the end of the last time step in the stress 

period). These constraints can be used to control excessive lowering (eq. 29) or mounding 

(eq. 30) of the potentiometric surface.

The second type of head constraint is drawdown of head at a specific location and 

stress period. Drawdowns are defined by and are equal to the difference between an 

initial head at location at the end of stress period , , and the head calculated 

at the location at the end of stress period after implementation of the optimal 

management strategy, :

 . (31)

The drawdown constraints are written as

(32)

and

 , (33)

where and are specified lower and upper bounds on drawdown at location 

 at the end of stress period (fig. 2B). The use of drawdown constraints in a GWM 

formulation requires that a reference withdrawal or injection rate be specified for each of 

the candidate well sites. These reference rates are used to calculate the initial head, , 

at each drawdown constraint location.

The third type of head constraint is a lower bound on the difference in head between 

two model locations and at the end of stress period :

, (34)

where and are heads calculated for locations and at the end 

of stress period and is a specified lower bound on the difference in head 

between locations 1 and 2 at the end of stress period (fig. 2C). The two specified locations 

may be vertically or horizontally separated and need not be adjacent cells. GWM requires 

that the head calculated at the second location be lower than the head at the first location by 

an amount that is at least .

hijkt ,, , hijkt ,, ,

l ≥

hijkt ,, , hijkt ,, ,

u ≤

hijkt ,, ,

l hijkt ,, ,u ijk , ,

t

ddijkt ,, ,

ijk , , t hijkt ,, , ( )0

t

hijkt ,, ,

ddijkt ,, , hijkt ,, , ( )0 hijkt ,, , = –

ddijkt ,, , dd ijkt ,, ,

l ≥

ddijkt ,, , dd ijkt ,, ,

u ≤

dd ijkt ,, ,

l dd ijkt ,, ,u

ijk , , t

hijkt ,, , ( )0

( ) ijk , , 1 ( ) ijk , , 2 t

hijkt ,, , ( )1 hijkt ,, , ( )2 – hdifijkt ,, , ( )1 2, ≥

hijkt ,, , ( )1 hijkt ,, , ( )2 ( ) ijk , , 1 ( ) ijk , , 2

t hdifijkt ,, , ( )1 2,

t

hdifijkt ,, , ( )1 2,
16 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Figure 2. Types of hydraulic-head constraints supported by GWM: A, Lower ( ) or upper ( ) 

bound on head ( ) at location at the end of stress period ; B, Lower ( ) or upper 

( ) bound on drawdown ( ) at location at the end of stress period ; is the 

initial head at ; C, Specified, , and model-calculated, ( ), difference 

in head between two model locations and at the end of stress period ; and D, Specified 

, and model-calculated , gradient in head between two model 

locations and at the end of stress period ; is the specified distance between the two 

locations.

hijkt ,, ,

l hijkt ,, ,u

hijkt ,, , ijk , , t dd ijkt ,, ,

l

dd ijkt ,, ,

u ddijkt ,, , ijk , , t hijkt ,, , ( )0

ijk , , hdifijkt ,, , ( )1 2, hijkt ,, , ( )1 hijkt ,, , ( )2 –

( ) ijk , , 1 ( ) ijk , , 2 t

gradijkt ,, , ( )1 2, hijkt ,, , ( )1

 hijkt ,, , ( )2 [ ] – ⁄ Δx

( ) ijk , , 1 ( ) ijk , , 2 t x Δ
Formulation of Ground-Water Management Problems with GWM 17

The fourth type of head constraint is a lower bound on the gradient in head between 

two model locations and at the end of stress period . The head gradient 

between two model locations ( ) can be written as

, (35)

where and are heads calculated for locations and , and is 

the distance between the two locations (fig. 2D). The constraint that is implemented in 

GWM is

. (36)

The user must specify two model locations (with head at the second location lower than 

head at the first location), a distance between the two model locations (which could be the 

distance between the midpoint of each of the two model cells), and a specified lower bound 

on the gradient ( ). The two model locations need not be adjacent.

Streamflow Constraints

Two types of streamflow constraints can be specified in a GWM management model, 

both of which are based on streamflows calculated by the Streamflow-Routing (STR1) 

Package developed for MODFLOW by Prudic (1989). The first type of constraint is an 

absolute lower and upper bound placed on streamflow at a specific stream location and 

stress period:

(37)

and

, (38)

where and are the specified lower and upper bounds on streamflow at stream 

location at the end of stress period . 

The second type of streamflow constraint is a streamflow depletion at a specific stream 

location and stress period. Streamflow depletion is defined by , and is equal to the 

difference between an initial streamflow at stream location at the end of stress period , 

, and the streamflow calculated at the location at the end of stress period after 

implementation of the optimal management strategy, :

. (39)

The streamflow-depletion constraints are written as

(40)

and

, (41)

where and are specified lower and upper bounds on streamflow depletion at 

location at the end of stress period . The use of streamflow-depletion constraints in a 

GWM formulation requires that a reference withdrawal or injection rate be specified for 

each of the candidate well sites. These reference rates are used to calculate the initial 

streamflow, , at each streamflow-depletion constraint location.

( ) ijk , , 1 ( ) ijk , , 2 t

gradijkt ,, , ( )1 2,

hijkt ,, , ( )1 hijkt ,, , ( )2 –

Δx ----------------------------------------------

hijkt ,, , ( )1 hijkt ,, , ( )2 ( ) ijk , , 1 ( ) ijk , , 2 Δx

hijkt ,, , ( )1 hijkt ,, , ( )2 –

Δx ---------------------------------------------- gradijkt ,, , ( )1 2, ≥

gradijkt ,, , ( )1 2,

Qsfr

Qsfr Qsf r

l ≥

Qsfr Qsf r

u ≤

Qsf r

l Qsf ru

r t

Qsdr

r t

Qsfr ( )0 t

Qsfr

Qsdr Qsfr ( )0 Qsfr = –

Qsdr Qsd r

l ≥

Qsdr Qsd r

u ≤

Qsd r

l Qsd ru

r t

Qsfr ( )0
18 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Complete Statement of Ground-Water Management 

Formulation Solved by GWM

Combining all objective and constraint functions described above, the complete statement of the ground-water management problem that can be formulated and solved by use 

of the GWM Process is

Maximize or minimize

, (42)

subject to

, (43)

, (44)

, (45)

, (46)

, (47)

, (48)

, (49)

, (50)

, (51)

, (52)

, (53)

, (54)

, (55)

, (56)

, (57)

, (58)

and

. (59)

βnQwnTQwn

 

n 1 =

N

∑ γ mExmTExm 

m 1 =

M

∑ κl

Il

l 1 =

L

+ + ∑

0 Qwn Qwn

u ≤ ≤

0 Exm Exm

u ≤ ≤

Qwn Qwn

l

I ≥ l

Exm Exm

l I ≥ l

apGVp b ≤ p

p 1 =

P

∑

apGVp b ≥ p

p 1 =

P

∑

apGVp 

p 1 =

P

∑ b = p

hijkt ,, , h ijkt ,, ,

l ≥

hijkt ,, , hijkt ,, ,

u ≤

ddijkt ,, , dd ijkt ,, ,

l ≥

ddijkt ,, , dd ijkt ,, ,

u ≤

hijkt ,, , ( )1 hijkt ,, , ( )2 – hdifijkt ,, , ( )1 2, ≥

hijkt ,, , ( )1 hijkt ,, , ( )2 –

Δx ---------------------------------------------- gradijkt ,, , ( )1 2, ≥

Qsfr Qsf r

l ≥

Qsfr Qsf r

u ≤

Qsdr Qsd r

l ≥

Qsdr Qsd r

u ≤
Solution of Ground-Water Management Problems with GWM 19

If binary variables are not used, then each binary variable in equations 42, 45, and 

46 is set to 0. If binary variables are present, then GWM automatically adds additional constraints, not shown here, that ensure that binary variables take only a value of 0 or 1, and 

that flow-rate and external variables can only be nonzero when their associated binary 

variable has a value of 1.

Solution of Ground-Water Management Problems with GWM

The complete ground-water management formulation defined by equations 42–59 

includes linear-programming, nonlinear-programming, and binary-programming aspects. 

A linear program is an optimization formulation in which the objective function and all 

constraints are linear functions of the decision variables. A nonlinear program is one in 

which the objective function and (or) one or more of the constraints are nonlinear functions 

of the decision variables. A binary program is one that includes binary decision variables, 

which in this case are the variables. An optimization formulation that combines a linear 

program with binary variables is called a mixed-binary linear program. A specific management formulation that is developed from the complete formulation will include one or more 

of these program types. GWM provides capabilities to solve linear, nonlinear, and mixedbinary linear formulations. Mixed-binary nonlinear formulations also can be solved if the 

nonlinearities of the system are mild. The presence of binary variables in a nonlinear formulation, however, can make the solution process difficult, and may lead to suboptimal 

solutions or spurious results.

This section describes the three approaches used by the RMS (Response Matrix Solution) Package of GWM to solve linear, nonlinear, and mixed-binary linear (and mildly nonlinear) formulations. Users of GWM are encouraged to read the section on solution of linear 

formulations before proceeding to the other solution techniques. Additional detailed information on the theoretical background of the solution techniques can be found in Gorelick 

and others (1993) and Ahlfeld and Mulligan (2000).

Linear Formulations

If the ground-water management problem to be solved is linear, then a highly efficient, 

reliable, and widely used method for solving linear optimization problems called the simplex method can be used to solve the management problem. Use of the algorithm requires 

that the ground-water management formulation be stated in the form of a linear program. 

With respect to the complete ground-water management formulation described by equations 42–59, this means, first, that the state variables of the ground-water system (heads and 

streamflows) must respond linearly to changes in the stress rates imposed at the flow-rate 

decision variables (that is, at each ), and, second, that no binary variables be present in 

the formulation. For the first condition to be met, the ground-water flow system must be 

simulated with constant transmissivity in every layer (LAYCON equal 0 or 2), all boundary 

conditions must be linear, and GWF Process packages that contain other nonlinearities 

must not be present. Head-dependent boundary conditions such as those specified along 

simulated streams and illustrated by the graphical relation in figure 1 can induce nonlinear 

responses. Therefore, although there may be instances where seepage rates between a 

ground-water system and adjoining streams are linear (or mildly nonlinear) functions of 

ground-water heads and imposed flow rates, the constraints on streamflow and streamflow 

depletion (eqs. 56–59) are dropped from the management formulation for strictly linear 

programs; these constraints are described below for nonlinear formulations.

Il

Il

Qwn
20 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Eliminating the binary variables from the formulation results in constraint equations 

45 and 46 being dropped and a modified objective function

Maximize or minimize

, (60)

subject to constraint equations 43, 44, and 47–55.

The final step in transforming the ground-water management problem into a form that 

can be solved by use of the simplex method is to establish functional relations between the 

stresses imposed at the managed wells (the flow-rate decision variables) and the resulting 

changes in heads at the constraint locations. These functional relations, which were 

described in a general form by equation 2, are necessary to rewrite the head constraints 

(eqs. 50–55) in terms of the flow-rate decision variables. This transformation requires use 

of a first-order Taylor series expansion to define head at each constraint location and stress 

period as a function of the new withdrawal and injection stresses:

, (61)

where

is head at constraint location and stress period for a new vector 

(that is, a new set) of withdrawal and injection flow rates having 

individual elements ; 

is head at constraint location and stress period for an original 

vector (that is, a base-condition set) of withdrawal and injection flow 

rates having individual elements ;

is the change in head at location and stress period for a change

in withdrawal or injection flow rate for the nth flow-rate decision 

variable, evaluated at the original vector of flow rates ; and

is the total number of flow-rate decision variables.

Equation 61 states that the head at each constraint location is equal to the head at the 

constraint location for a base condition of withdrawal and injection rates plus the sum 

of the head changes that result from the changes in withdrawal or injection rate at each of 

the well sites. The linear summation defined by the right-hand side of the equation 

reflects the assumed linearity of the ground-water flow system and consequent linear 

response of heads to changes in withdrawal or injection rates. The partial derivatives in 

equation 61, ( ), which are called the response coefficients, provide information on the response of ground-water heads to stresses at each withdrawal or injection site. 

For linear systems, each partial derivative is a constant whose value does not change with 

changes in the distribution of withdrawal and injection rates throughout the ground-water 

flow system.

Equation 61 is substituted for each head term in the head constraints described by 

equations 50 and 51; similar substitutions are made for other head constraints described by 

equations 52–55. For a more detailed description of the mathematical steps involved in 

these substitutions, the reader is directed to Ahlfeld and Mulligan (2000, p. 64–66).

βnQwnTQwn 

n 1 =

N

∑ γ mExmTExm

m 1 =

M

+ ∑

hijkt ,, , ( ) Qw h

0

ijkt ,, , Qw0 ( ) hijkt ,, , ∂

Qwn ∂

--------------- Qw0 ( ) Qwn Qwn

0 ( ) –

n 1 =

N

= + ∑

hijkt ,, , ( ) Qw ijk , , t

Qw

Qwn

h

0

ijkt ,, , Qw0 ( ) ijk , , t

Qw0 Qwn

0

hijkt ,, , ∂

Qwn ∂

--------------- Qw0 ( ) ijk , , t

Qw0

N

Qw0

n

hijkt ,, , ∂ Qwn ⁄ ∂
Solution of Ground-Water Management Problems with GWM 21

Calculation of Response Coefficients

The partial derivatives that define the response coefficients are not calculated directly; 

instead, they are approximated by a first-order, finite-difference perturbation method. The 

derivative of head with respect to each flow-rate decision variable is approximated by the 

forward-difference equation

, (62)

where is the perturbation value for the nth flow-rate decision variable and 

is the head at constraint location and stress period computed by using a vector of 

withdrawal and injection stress rates that differs from the original vector of stress 

rates only in the nth element, which is changed by an amount .

To calculate each response coefficient defined by equation 62, the GWF Process of 

MODFLOW is run a total of N+1 times. In the first run, which is called the base-condition 

run, the head is calculated at each constraint location and stress period for the set of basecondition withdrawal and injection rates [that is, each ]. In each of the remaining N runs, the head for each constraint location is calculated on the basis of the change 

(perturbation) in the withdrawal or injection rate for the nth flow-rate decision variable. For 

each of these runs, the withdrawal or injection rate at each of the remaining N-1 well sites 

is kept at the base-condition value. The N runs and consequent computations of the 

response coefficients result in a matrix of response coefficients that is used to solve the 

optimization problem.

The perturbation values are calculated from the equation

, (63)

where is a user-specified perturbation variable and is the specified upper bound on 

withdrawal or injection rate at well site . The perturbation variable can be a positive or 

negative value: a positive value of implies an increase in flow rate for (referred to 

as a forward-difference calculation), whereas a negative value of implies a decrease in 

flow rate for (referred to as a backward-difference calculation).

As described above, calculation of the response coefficients requires conducting one 

flow-process run to obtain base values of the state variables and one additional flowprocess run for each perturbed flow-rate decision variable. Because these runs are done 

automatically by GWM, it is important that the GWF Process input files produce runs that 

are stable for a range of values of the flow-rate decision variables. If the flow-process 

simulation is subject to convergence failure, excessive dewatering, or other instabilities, 

GWM may not be able to successfully compute response coefficients.

GWM has mechanisms for adjusting the perturbation value to attempt to overcome 

instabilities in the GWF Process. For each flow-process run, GWM checks for successful 

completion of the GWF Process. If failure occurs on a perturbation run, then GWM changes 

the perturbation value by a user-specified factor and re-executes the GWF Process. For a 

single flow-rate decision variable, re-adjustment of the perturbation value can be repeated 

until a successful run is achieved or a user-specified maximum is reached. See “Description 

of Selected Conventions, Options, and Variables in GWM” for more details on this 

feature.

hijkt ,, , ∂

Qwn ∂

--------------- Δhijkt ,, ,

ΔQwn

---------------- hijkt ,, , QwΔn ( ) h

0

ijkt ,, , Qw0 – ( )

QwΔn

≈ = --------------------------------------------------------------------

QwΔn hijkt ,, , QwΔn ( )

ijk , , t

QwΔn

Qw0 QwΔn

h

0

ijkt ,, , Qw0 ( )

QwΔn

QwΔn δ

0

Qwn

u =

δ

0 Qwn

u

n δ

0

δ

0 QwΔn

δ

0

QwΔn
22 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Accuracy and Precision of the Response Coefficients

The response coefficients generated by the MODFLOW GWF Process are a critical 

link between the physics of the ground-water flow system and the results of a ground-water 

management model represented in GWM. As a consequence, the accuracy and precision of 

the response coefficients play an important role in the solution of a ground-water management problem. The accuracy of the response coefficients—that is, their ability to reflect the 

actual response of the aquifer—depends on at least two factors.

First, the accuracy of the heads calculated by MODFLOW for a particular groundwater flow system is a reflection of the accuracy with which the geologic framework, 

hydraulic properties, boundary conditions, and other flow-system characteristics are 

known and represented in the model. Clearly, the accuracy of the response coefficients 

increases as the level of understanding of the flow system and representation of the flow 

system with a numerical model improve.

The second issue that affects response-coefficient accuracy is approximation of the 

partial-differential equation of ground-water flow by a set of finite-difference equations. As 

described previously in this report (“Numerical Modeling of Ground-Water Flow with 

MODFLOW”), solution of the finite-difference equations that are used by MODFLOW to 

describe the head at each cell within the model domain gives only approximate values of 

the “true” head distribution in the aquifer as described by the partial-differential equation 

of ground-water flow. The error that is introduced by this approximation is proportional to 

the size of the grid cells that compose the discretized model domain and the size of the time 

steps used for transient conditions. In general, the approximation error increases as the grid 

becomes coarser and (or) the time steps become longer. Therefore, approximation error can 

be decreased by decreasing the size of the model grid cells and (or) time steps, although 

such actions typically increase computer run time.

The precision of the response coefficients is an indication of their ability to reflect the 

actual response of the calculated system state to changes in stress. A measure of responsecoefficient precision is the number of significant digits in the value of the response coefficient. Maintaining adequate precision in computed response coefficients is essential for 

successful solution of the optimization problem by the RMS Package. Even if the flowprocess run is not particularly accurate, it is important that the precision of each response 

coefficient be maintained. Precision of the response coefficients is affected in part by the 

size of the stress-rate perturbation values (that is, the size of each ) used in equation 

62. Although the response coefficients for the linear systems considered here are constants 

whose values are independent of the size of the perturbation values, the precision of the 

response coefficients, and therefore the precision of the management solution, depends on 

the number of significant digits carried for each real number in the computation of each 

response coefficient. The relation between perturbation size and response-coefficient 

precision is related to roundoff error when the difference in heads is taken in the numerator 

of equation 62. If the two computed heads are very close in value, then significant precision 

can be lost. The issue is even more important when considering the use of perturbation for 

calculating the response coefficients for head-difference or gradient constraints (eqs. 54 

and 55).

Most of the methods available to solve the flow-process finite-difference equations are 

iterative. Values of head are iteratively generated until the maximum calculated change in 

head at any model cell is less than a specified convergence criterion between iterations. The 

precision of the resulting heads can be estimated to be of the same magnitude as the convergence criterion. As a result, the precision of the response coefficients depends upon the 

convergence criterion used by the flow process.

QwΔn
Solution of Ground-Water Management Problems with GWM 23

There are two approaches for improving the precision of the response coefficients calculated for linear systems. The first approach is to lower the head convergence criterion in 

the GWF Process iterative solver, although this is likely to increase the computer time 

required to solve each GWF Process run. The second approach is to use relatively large perturbation values to ensure large head differences in the numerator of equation 62 (Reifler 

and Ahlfeld, 1996). Perturbation values equal to 500 percent of the expected solution stress 

rate are recommended as an initial guess; the perturbation values can then be varied in a 

series of GWM runs to evaluate how the management solutions are affected, if at all, by the 

size of the stress-rate perturbations.

Simplex Algorithm

Generation of the response-coefficient matrix by the RMS Package completes the 

transformation of the ground-water management problem into a form that can be solved 

using linear programming techniques. The revised linear formulation of the ground-water 

management problem consists of the objective function defined by equation 60, subject to 

constraints defined by equations 43, 44, 47–49, and revised constraints 50–55, in which 

heads have been rewritten in terms of the flow-rate decision variables (eq. 61). The revised 

linear formulation can be expressed in an equivalent vector form as

Minimize (64)

subject to (65)

, (66)

where Z is the value of the objective function; is a transposed column vector of objectivefunction coefficients associated with the decision variables; is a column vector of decision variables with upper bounds ; is a matrix of coefficients that includes the response 

matrix for head; and is a column vector of right-hand-side coefficients associated with 

the constraints. Constraints that have the form of inequalities are transformed to equalities 

by the addition of slack and surplus variables. The RMS Package solves the minimization 

problem described by equations 64–66; for maximization problems, a simple transformation is done internally by multiplying the objective function ( ) by -1.

The RMS Package solves the ground-water management problem by using the simplex 

algorithm, which iteratively solves for the optimal solution, . The mathematical details 

of the simplex algorithm as implemented in GWM are described in Appendix 2. Although 

the simplex method is reliable and stable, numerical conditions can arise in which the algorithm will cycle among points and will not converge. Such conditions are rare for wellposed linear problems, but to prevent indefinite cycling, the user must specify the maximum number of iterations that is allowed by the algorithm (input variable LPITMAX). A 

recommended value for LPITMAX is ten times the number of constraints, although experience suggests that the number of iterations to convergence is typically less than two times 

the number of functional constraints (Ahlfeld and Mulligan, 2000); functional constraints 

include all constraints that are not simple bounds on the decision variables.

There are four possible outcomes to the solution of each linear formulation. The first 

outcome is a single optimal solution for the elements of the decision-variable vector that 

minimizes (or maximizes) the objective function and meets all of the constraints. Specifically, this solution consists of a set of optimal withdrawal and injection rates for each flowrate decision variable and a set of optimal flow rates for each external variable . 

The second possible outcome is one in which there is no set of decision variables that 

Z c

t = x

Ax = b

0 ≤ ≤ x u

c

t

x

u A

b

c

t

x

x

*

x

Qwn Exm
24 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

simultaneously satisfies all constraints; in this case, the formulation is said to be infeasible. 

The third outcome is one in which the formulation is determined to be unbounded and the 

optimal objective function is either positive or negative infinity. The last possible outcome 

is one in which there are multiple optimal solutions to the formulation, in which the values 

of the decision variables will be different but the objective-function values will be identical. 

Ahlfeld and Mulligan (2000) describe some of the ways in which infeasible, unbounded, 

and multiple-optima problems can arise in ground-water management problems, as well as 

ways to prevent these outcomes.

One of the benefits of using the simplex method to solve a linear program is that specific sensitivity-analysis information can be determined easily by using the information in 

the final iteration of the algorithm. Such information includes the sensitivity of the optimal 

solution to the objective-function coefficients and the right-hand-side values of the constraints. Sensitivity analysis is implemented in GWM for linear problems and reported 

along with the optimal solution. See the "Output Files" section of this report for additional 

information about management-model results.

Nonlinear Formulations

Nonlinearities can arise in the ground-water management formulation (eqs. 42–59) 

indirectly as a result of two common characteristics of ground-water flow models. The first 

is the presence of layers in which transmissivity is a function of head (LAYCON 

equals 1 or 3). As noted previously (p. 4), water-table conditions cause a nonlinear relation 

between the position of the water table and withdrawal or injection stresses. The second 

characteristic is the presence of head-dependent boundary conditions such as streams, 

drains, evapotranspiration, and so forth. These boundary conditions can create nonlinear 

relations between ground-water heads and flow rates to or from the simulated boundary; 

an example of a nonlinear relation between ground-water heads and a simulated headdependent boundary condition is illustrated in figure 1 for seepage between an aquifer and 

stream. Other GWF Process packages may also induce nonlinear responses.

The RMS Package provides a solution approach designed to address many of the nonlinear features that can arise in ground-water management problems. For purposes of 

describing this approach, a subset of the general management formulation will be considered. This formulation excludes binary variables and consists of the objective function 

described by equation 60 subject to constraints 43, 44, 47–55, and the streamflow and 

streamflow-depletion constraints described by equations 56–59.

The approach for solving nonlinear formulations is based on the sequential linear programming (SLP) algorithm and is referred to as the SLP approach. It is based on repeated 

linearization of the nonlinear features in the management problem, and is implemented by 

recalculating the response matrix for each sequential linear program. The first-order Taylor 

series expansion for head given by equation 61 is assumed to be accurate for each sequential linear program, but in contrast to the linear case, the vector of base flow rates changes 

at each iteration. The head at any location is estimated by

, (67)

where the superscript represents an iteration level, so that is the head obtained 

when the vector (set) of withdrawal and injection rates is applied. For values of 

that are close to , the error in this approximation is small, and the values of predicted by this equation are relatively accurate.

hijkt ,, , ( ) Qw hijkt ,, ,

v Qwv ( ) hijkt ,, ,

v ∂

Qwn

v ∂

--------------- Qwv ( ) Qwn Qwn

v ( ) –

n 1 =

N

= + ∑

v hijkt ,, ,

v

Qwv Qw

Qwv hijkt ,, ,
Solution of Ground-Water Management Problems with GWM 25

A first-order Taylor series expansion also is used to establish functional relations 

between streamflow and the managed withdrawal or injection stresses. These functional relations, which were described in a general form by equation 4, are necessary to rewrite the 

streamflow constraints (eqs. 56 and 57) in terms of the flow-rate decision variables. The firstorder Taylor series approximation of streamflow for each sequential linear program is

, (68)

where

is streamflow at the rth stream reach and stress period for a new vector 

(that is, a new set) of withdrawal and injection flow rates, , having 

individual elements ; 

is streamflow at the rth stream reach and stress period for the vector of 

withdrawal and injection flow rates at iteration level , , having 

individual elements ;

is the change in streamflow at the rth stream reach and stress period for a

change in withdrawal or injection flow rate for the nth flow-rate 

decision variable, evaluated for the set of withdrawal and injection rates

; and

is the total number of flow-rate decision variables.

The partial derivatives in equation 68, ( ), are the response coefficients for the 

response of streamflow to stresses at each withdrawal or injection site. Equation 68 is substituted for each streamflow term in constraints 56 and 57; similar substitutions are made for the 

streamflow-depletion constraints described by equations 58 and 59.

At each iteration of the SLP algorithm, a linear program is constructed on the basis of the 

first-order Taylor series approximation and the associated response coefficients, and the formulation is solved by using the simplex method described previously. Because the head and 

streamflow responses may be nonlinear, the response coefficients in equations 67 and 68, 

 and , may no longer be constant. Therefore, the response coefficients for heads and streamflows must be recalculated at each iteration . This calculation 

uses a new vector of base-condition flow rates that is derived from the optimal flow rates 

obtained from the linear-program solution from the prior iteration.

The sequential process is continued until two convergence criteria are met. The first 

requires that the change in flow-rate variable values from the prior iteration to the current iteration be less than a fraction of the magnitude of the flow-rate variables at the current iteration:

 . (69A)

The infinity norm is used in this expression and is specified by the user as input variable 

SLPVCRIT. Note that 1 is added to the norm on the right-hand side of equation 69A. While it 

is likely that in most cases the norm of decision variables will be much larger than 1, there 

may be cases where the norm of the decision variables is 0 or nearly 0. In these cases, the addition of 1 ensures that a reasonable convergence criterion is used. The input variable can be 

considered a measure of the number of correct figures desired in the solution. For example, a 

value of 10-5 indicates that the solution will be correct to five significant digits or to the fifth 

decimal place, if the norm of the decision variables is less than 1 (Gill and others, 1981).

Qsfr( ) Qw Qsfr

v Qwv ( ) Qsfr

v ∂

Qwn

v ∂

------------ Qwv ( ) Qwn Qwn

v ( ) –

n 1 =

N

= + ∑

Qsfr( ) Qw

Qw

Qwn

Qsf r

v Qwv ( )

v Qwv

Qwn

v

Qsf r

v ∂

Qwn

v ∂

------------- Qwv ( )

Qwv

N

Qsf r

v ∂ Qwnv ⁄ ∂

hijkt ,, ,

v ∂ Qwnv ⁄ ∂ Qsfrv ∂ Qwnv ⁄ ∂

v

 Qwv 1 + Qwv – ε1 1 Qwv 1 + ≤ ( ) + 

ε1

ε1
26 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

The second criterion requires that the change in objective function value, Z, be less 

than a specified fraction of the magnitude of the objective function value

. (69B)

The fraction, , is specified by the user as input variable SLPZCRIT.

The SLP algorithm can be summarized as:

Step 0: Set and set ;

Step 1: Compute response coefficients from base vector of rates ;

Step 2: Assemble and solve the linear program; solution vector is assigned to ;

Step 3: If convergence tests (eqs. 69A and 69B) are met, stop; else,

Step 4: Set and go to Step 1.

Response coefficients are calculated for each iteration by use of a forwarddifference equation similar to that in equation 62. For heads, this equation is

, (70)

where is the perturbation value for the nth flow-rate decision variable for iteration , 

and is the head at constraint location and stress period computed by 

using a vector of withdrawal and injection stress rates that differs from the previous 

vector of stress rates only in the nth element, which is changed by an amount . 

An equation similar to 70 also is written for the streamflow-response coefficients 

.

In the first iteration ( ), response coefficients are calculated on the basis of a 

set of base-condition withdrawal and injection rates ( ) and initial perturbation values 

 defined by the user in the same way as they are defined for linear formulations. 

Because of the nonlinearity of the problem and the need to solve the problem iteratively, 

the specified base-condition withdrawal and injection rates should be as close as possible 

to the expected optimal solution.

The same issues that affect the accuracy and precision of the response coefficients for 

linear formulations also apply for nonlinear formulations. However, because of the nonlinearity of the systems discussed here, the user must be careful to choose initial perturbation 

values small enough to ensure that the response coefficients for heads and streamflows calculated by the forward-difference equations are accurate. If the perturbation value used is 

too large, then the derivative approximation may be poor, which can lead to convergence 

problems in the sequential linearization process. With small perturbation values, heads 

determined by the GWF Process must be calculated with high precision to ensure adequate 

precision in the response matrix.

 Zv 1 + Z

v

 ε2 1 Zv 1 + – ≤ ( ) + 

ε2

v 0 = Qw0

Qwv

Qwv+1

v = v 1 +

v

hijkt ,, ,

v Δ

Qwn

v Δ

---------------- hijkt ,, ,

v QwΔnv ( ) hijkt ,, ,v Qwv – ( )

QwΔn

v = ------------------------------------------------------------------

QwΔn

v v

hijkt ,, , QwΔn

v ( ) ijk , , t

QwΔn

v

Qwv QwΔn

v

Qsfr

v ∂ Qwnv ⁄ ∂

v 0 =

Qw0

QwΔn ( )initial
Solution of Ground-Water Management Problems with GWM 27

Because of the importance of the perturbation value to the success of the SLP algorithm, several features are provided to adjust the perturbation value automatically during 

the course of the algorithm. At each iteration, the perturbation values are calculated by 

using

, (71)

where is the perturbation variable for iteration and is the specified upper bound 

on the withdrawal or injection rate at the nth well site. The initial value of the perturbation 

variable is set with , as described for linear problems. Additional input variables are 

available to cause the perturbation variable to change in subsequent iterations. The following formula (Minihane, 2002) is used to compute the perturbation variable at each iteration:

, (72)

where is the iteration level, is an initial perturbation variable for the nth well site, 

 is the minimum perturbation variable for the nth well site, and is a scaling 

factor that determines the rate of decrease in the perturbation variable. Note that the superscript on indicates exponentiation, so that when is set to a value larger than 1, the 

perturbation variable decreases as the sequential process proceeds and approaches the 

minimum value .

During each iteration of the SLP algorithm, a response matrix is calculated. As 

described in the section “Calculation of Response Coefficients,” it is possible that a flowprocess run may fail. If a run using the perturbed stress rates does fail, then the perturbation 

value is changed automatically during an SLP iteration. This has no effect on the perturbation value at the next SLP iteration. That is, at the beginning of each iteration, is calculated according to equations 71 and 72 regardless of the changes that may have been made 

to the perturbation value at the prior iteration.

At each iteration of the SLP algorithm, the response coefficients are calculated from a 

different set of base conditions. If the first base run was successful but a later base run fails, 

then it is possible to adjust the base vector of flow-rate variables rather than terminate the 

algorithm. This is done by use of a relaxation parameter , which moves the base solution 

closer to the prior (successful) base solution. For each of the flow-rate decision variables, 

a temporary base value is calculated according to

, (73)

where is the temporary value of the flow-rate decision variable. The temporary base 

solution, which consists of values for each flow-rate decision variable calculated according 

to equation 73, is used for a new attempt at a successful base run. If the flow-process run 

fails again, a new temporary base solution is determined with the prior temporary base solution substituted for in equation 73. A sequence of applications of equation 73 results 

in a sequence of temporary base solutions that move closer to the prior base solution .

QwΔn

v δ

v

Qwn

u =

δ

v v Qwn

u

δ

0

δ

v

 δinitial δminimum –

δscale ( )v ------------------------------------ δ = + minimum

v δinitial

δminimum δscale

δscale δscale

δminimum

δ

v

α

Q

ˆ wn

v

 1 ( ) – α Q

ˆ wn

v

 αQ

ˆ wn

v 1 – = +

Q

ˆ wn

Q

ˆ wn

v

Q

ˆwn

v 1–
28 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Mixed-Binary Linear Formulations

In GWM, mixed-binary linear formulations are those that include the binary variables 

 in a linear formulation. When binary variables are specified, the RMS Package uses a 

branch and bound solution method to solve the formulation. The branch and bound method 

solves a series of linear programs. In each linear program, the binary variables are forced 

to an assumed binary value (0 or 1) or are allowed to vary between zero and one (that is, 

the binary variables are relaxed). Comparisons among the solutions of each of these linear 

programs leads to a solution without the necessity of checking all possible combinations of 

binary-variable values (Nemhauser and Wolsey, 1988; Hillier and Lieberman, 2001). The 

simplex algorithm is used to solve each of the linear programs. Implementation of the 

branch and bound algorithm is described in Appendix 2.

The branch and bound algorithm is based on repeated solution of linear programs. The 

user must specify the maximum number of iterations, or linear-program solutions, allowed 

(input variable BBITMAX). This parameter defines the computer storage space allocated 

for storing information related to each linear-program solution. Because a linear program 

is solved at each iteration of the branch and bound algorithm, the user must also specify a 

value of LPITMAX for a mixed-binary linear formulation. Output describing the progress 

of the branch and bound algorithm and the results of each linear program solved can be 

obtained by using the BBITPRT input variable.

Input Instructions and Output Files

This section describes the input instructions and output files for GWM. Before solving 

an optimization problem with GWM, the user must have already developed a ground-water 

flow model of the study area based on the MODFLOW GWF Process. As described in the 

Introduction of this report, the GWM Process can only be used with the Global and GWF 

Processes; it cannot be used with the Ground-Water Transport, Observation, Sensitivity, or 

Parameter-Estimation Processes. If, however, the GWF Process input files contain parameters defined for the Parameter-Estimation Process and the input file for the Sensitivity 

Process is included in the name file, then the specified parameter values will be used in the 

GWF Process (see “Sample Problem 2: SEAWATER”). The Global Process controls overall program operation and sets up data structures that can be used by all MODFLOW processes. The only Global Process file that needs to be modified for a GWM run is the name 

file; the required modifications to the file are described in detail below.

Name File

The name file contains the names of most input and output files used by MODFLOW, 

and determines which MODFLOW program options are activated (Harbaugh and others, 

2000). The name file is read on unit 99. The name file contains one record of information 

for each input and output file. (A record is a line in a file.) Each record consists of three 

variables, which are read in free format; the length of each record must be 199 characters 

Il
Input Instructions and Output Files 29

or less. Comment records can be used in the name file and are indicated by the # character 

in column one; comment records can be placed anywhere in the name file. Any text characters can follow the # character. Comment records have no effect on the run; their purpose 

is to allow users to provide documentation about a particular run.

Each record has the following format:

Ftype Nunit Fname

Explanation of the variables:

Ftype—is the file type. Ftype may be entered in all uppercase, all lowercase, or any combination thereof.

Nunit—is the FORTRAN unit to be used when reading from or writing to the file. Any legal 

unit number on the computer being used can be specified except units 96–99. Also, the unit 

number for the file must be unique; that is, it cannot be equal to any of the unit numbers 

used for other files specified in the name file.

Fname—is the name of the file, which is a character value. Pathnames may be specified as 

part of Fname.

GWM uses both the GLOBAL and LIST output files (see further discussion in the 

“Output Files” subsection below); therefore, both of these file types must be specified in 

the name file. The GLOBAL record must be the first non-comment record in the name file, 

and the LIST record must be the second non-comment record. A record also must be added 

to specify that the GWM Process is active. The record must specify Ftype GWM (bold text 

indicates a MODFLOW keyword). The file identified in this record contains information 

needed for the GWM run.

Example input records for the GLOBAL, LIST, and GWM file types are:

GLOBAL 1 global.gwm

LIST 6 list.gwm

GWM 55 input.gwm

GWM Process Files

Input files for the GWM Process consist of the GWM file and several supporting files. 

The GWM file is used to activate the GWM Process and to identify the files that will be 

opened for the GWM run. Four types of information about the management problem are 

specified in the input files: the decision variables, objective function, constraints of the 

management problem, and the solution and output-control parameters. With the exception 

of the GWM file, each of the input files is read from FORTRAN unit 99; each file is opened 

on unit 99 just prior to reading input from the file and closed immediately after reading the 

input. 

For the most part, the general structure of the input formats for the GWM Process files 

are consistent with other MODFLOW processes; users of GWM should review the input 

instructions for MODFLOW given in Harbaugh and others (2000). Input for each GWM 
30 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Process file is grouped by numbered items, and each item consists of input variables. The 

first item in each of the input files is Item 0 (#Text), which can be used for comment lines 

but is optional. Some items consist of several variables, and the item can be repeated multiple times. The input data for each item must start on a new record. Each record is limited 

to a length of 199 characters. An input variable may include a single value or multiple values. Variables are defined after all the items are listed.

Each input variable has a data type, which can be Real, Integer, or Character. Integers 

are whole numbers and must not include a decimal point or exponent. Real numbers can 

include a decimal point and an exponent; if no decimal point is included in the entered 

value, then the decimal point is assumed to be at the right side of the value. Any printable 

character is allowed for character variables. Unlike the GWF Process, variables used by 

GWM that start with the letters I–N are not necessarily integers and those that start with the 

letters A–H and O–Z are not necessarily real numbers. Data types are specified for each 

input variable.

Free formatting is used for GWM input. With free format, values are not required to 

occupy a fixed number of columns in a record. Each value can occupy one or more columns 

as required to represent it; however, the values must still be included in the prescribed 

order. One or more spaces, or a single comma optionally combined with spaces, must separate adjacent values. Also, a numeric value of zero must be explicitly represented with 0 

and not by one or more spaces when free format is used.

Units of values used in the GWM Process should be consistent with the units used in 

the other MODFLOW data-input files.

Description of Selected Conventions, Options, and Variables in GWM

As an aid in the preparation of GWM input files, some of the conventions, options, and 

variables of GWM are described below. Program variables are shown in plain upper-case 

text and file names are shown in bold upper-case text.

Sign conventions for flow-rate and external decision variables and their 

coefficients: As described previously, all flow-rate and external decision variables in 

GWM are treated as positive values, whether they represent a withdrawal, injection, export, 

or import of water. A flow-rate decision variable is defined with input variable FTYPE in 

the Decision Variables (DECVAR) file: If FTYPE is W, the variable is used for withdrawal; if FTYPE is I, the variable is used for injection. Because the user specifies whether 

a variable is for withdrawal or injection, GWM will know how to treat the variable 

internally. Therefore, when specifying the minimum or maximum flow rate at a site (input 

variables FVMIN and FVMAX), as well as the reference flow rate for the variable (input 

variable FVREF), values greater than 0 should be used for both withdrawal and injection 

variables. Likewise, GWM will know whether an external decision variable is an import or 

export of water by the definition of input variable ETYPE in the DECVAR file: If ETYPE 

is IM, the variable is a source (import) of water; if ETYPE is EX, the variable is a sink 

(export) of water. Therefore, the user should specify values greater than 0 for the minimum 

(EVMIN) and maximum (EVMAX) flow rates for either type of external decision variable. 

When specifying the objective-function coefficients for the decision variables with 

input variables FVOBJC, EVOBJC, and BVOBJC in the Objective Function (OBJFNC) 

file, the user must specify the sign of each coefficient to achieve the desired objective-function form. For example, if the objective is to maximize economic benefit of the water withdrawn from or imported to a basin, and if there is a benefit of $1 per unit of water withdrawn 

at a particular flow-rate site but a cost of $1 per unit of water imported from an external 
Input Instructions and Output Files 31

variable, then the user would specify a positive coefficient for the withdrawal site and 

a negative coefficient for the external variable. The same reasoning applies to linearsummation constraints specified in the SUMCON file. For example, if the formulation 

requires that the difference in volumetric flow rates between an import and export 

be greater than or equal to a specified value , then the user would define a constraint with 

positive coefficient on the import variable and a negative coefficient on the export variable, 

that is, .

Lower bounds on flow-rate and external decision variables: Two approaches are 

used to define lower bounds on flow-rate and external decision variables, depending on 

whether or not the decision variables have been associated with binary variables. If a flowrate or external decision variable has been associated with a binary variable (that is, eqs. 12 

and 13 apply), then nonzero lower bounds can be defined for the variable in the DecisionVariables Constraints (VARCON) file by using input variables FVMIN for flow-rate decision variables and EVMIN for external decision variables.

If a flow-rate or external decision variable has not been associated with a binary variable (that is, eqs. 10 and 11 apply), then a lower bound of 0 must be specified for the variable by use of FVMIN and EVMIN in the VARCON file. The user has the option, however, 

of specifying a nonzero lower bound for a flow-rate or external decision variable not 

associated with a binary variable by use of linear-summation constraints. In this case, the 

user would specify FVMIN or EVMIN equal to 0 and then define a linear-summation 

constraint in the SUMCON file. For example, a minimum withdrawal rate of 1,000 ft3/d 

for hypothetical flow-rate decision variable Q1 could be specified in the SUMCON file by 

use of equation 15; specifically, .

Background stresses and reference flow rates: As described in the definition of 

flow-rate decision variables (p. 8), the user can specify withdrawal or injection flow rates 

at unmanaged wells in a GWM run. These flow rates are referred to as background stresses, 

and are specified by use of either the WEL or MNW Packages of the GWF Process of 

MODFLOW. Flow rates specified in these packages will not be modified during the GWM 

run.

Reference flow rates must be specified for each flow-rate decision variable when either 

drawdown constraints (eqs. 32 and 33) or streamflow-depletion constraints (eqs. 40 and 41) 

are used in a GWM run. Reference flow rates are specified by use of input variable FVREF 

in the Decision-Variables Constraints (VARCON) file. Note, however, that a user could 

implicitly define a reference flow rate as a background stress by use of an unmanaged well 

in either the WEL or MNW Packages. In this case, the user would set FVREF equal to 0 

for each decision variable for which a reference flow rate has been defined in the WEL or 

MNW Packages.

Variables related to calculation of response coefficients: The user must specify a 

number of input variables that are used by GWM to calculate the response coefficients; 

most of these variables are specified in the Solution and Output-Control Parameters 

(SOLN) file. The first of these variables are the perturbation variables ( for linear formulations and , , and for nonlinear formulations). For linear formulations, the 

perturbation variable in equation 63 must be specified with input variable DELTA. The 

perturbation variable DELTA can be a positive or negative value: a positive value implies 

an increase in flow rate (referred to as a forward-difference calculation), whereas a negative 

value implies a decrease in flow rate (referred to as a backward-difference 

calculation). Moreover, DELTA is not limited to values between +1.0 and -1.0. The user 

may need to experiment with different values of DELTA in a series of GWM runs to 

Im1 Ex1

b

Im1 Ex1 – ≥ b

Q1 1000 ≥

δ

0

δinitial δminimum δscale

δ

0
32 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

determine the value of DELTA that is most appropriate for a particular problem. Guidelines 

for the selection of DELTA are discussed in the "Accuracy and Precision of the Response 

Coefficients" section of this report.

If the user selects the SLP solution type in the SOLN file (that is, the optimization formulation is solved using sequential linear programming), then the three perturbation variables , , and in equation 72 must be specified with input variables DINIT, 

DMIN, and DSC, respectively. The selection of DINIT ( ) follows the same reasoning 

as the selection of DELTA in the linear case (see "Accuracy and Precision of the Response 

Coefficients" section of this report). DMIN ( ) should be selected such that the final 

perturbation value is large enough to produce several significant digits in the response 

matrix, but small enough to secure an accurate approximation for the nonlinear response. 

An initial estimate for DMIN of 0.5 percent of the expected optimal solution value is suggested. The scaling factor DSC ( ) is normally selected to be larger than 1.0 so that the 

perturbation parameter decreases with increasing iterations. A value of 5.0 has been shown 

to give good results (Minihane, 2002). As with the selection of DELTA, the user may need 

to experiment with different values of DINIT, DMIN, and DSC to determine the set of variables that works best for a particular GWM problem.

All solution types require specification of an upper bound on the withdrawal or injection rate at each flow-rate decision variable (variable ) for calculation of the perturbation values for linear formulations (eq. 63) or for nonlinear formulations 

(eq. 71). Upper bounds on flow-rate decision variables are specified with input variable 

FVMAX in the Decision-Variables Constraints (VARCON) file.

All solution types also require specification of the base-condition withdrawal or injection flow rates for the flow-rate decision variables that are used for calculation of the 

response coefficients (that is, ). For convenience, two options are provided for specifying the base-condition flow rates. The particular option that is selected by the user is 

defined by use of input variable IBASE in the SOLN file: if IBASE is set to 0, then input 

variable FVREF in the VARCON file will be used to specify the base-condition flow rate 

for each decision variable; if IBASE is set to 1, then input variable FVBASE in the SOLN

file will be used to specify the base-condition flow rate for each decision variable.

Testing of response-coefficient precision: Users of GWM should be aware of the 

precision of the response matrix. As an aid in determining the precision of the response 

matrix that is generated, GWM output includes the average number of significant digits in 

the computed response matrix. The number of significant digits in a response coefficient is 

estimated by dividing the change in system state caused by the perturbation by the head 

convergence criterion used by the flow process. The number of digits to the left of the decimal point in the resulting ratio is used as a measure of the number of significant digits. This 

ratio is computed by using the GWF Process variable HCLOSE to represent the headconvergence criterion. For this ratio to be meaningful, the value of HCLOSE must be properly specified. When using either the Preconditioned Conjugate-Gradient (PCG), Direct 

Solver (DE4), or other flow-process solution packages that have multiple convergence 

criteria, the user should ensure that HCLOSE is assigned a meaningful value.

Because heads will tend to respond very little to pumping from flow-rate variables that 

are far away, the corresponding response coefficients will be small and may have poor precision. As long as other flow-rate variables produce significant response with adequate precision, however, the response matrix should be sufficiently precise to achieve a solution. 

The RMS Package tests each column of the response matrix as it is computed. At least one 

response coefficient in the column must have sufficient precision to predict a meaningful 

response. In other words, each flow-rate variable must have a numerically significant effect 

δinitial δminimum δscale

δinitial

δminimum

δscale

Qwn

u

QwΔn QwΔn

v

Qw0
Input Instructions and Output Files 33

on at least one constraint. The test is conducted by requiring that the number of significant 

digits in the ratio of system state to HCLOSE be greater than or equal to the user-defined 

variable NSIGDIG for at least one entry in the column. If a column fails this test, then 

GWM can automatically recalculate the response-matrix column with a new perturbation 

value.

Correction of a failed GWF Process run: The GWF process may fail during 

a base-condition run or during one of the perturbation runs. This can occur for several 

reasons, such as the GWF Process solution algorithm fails to converge, a cell included in a 

constraint dewaters, or GWM detects that the simulated perturbed state value is inadequate 

or imprecise. If the flow process fails during a base-condition run for solution types NS, 

MPS, or LP, then the program will terminate; this will occur if the base values of the flowrate decision variables supplied by the user do not produce a successful run. In this case, 

the user must modify the base-condition stresses and (or) GWF solution parameters to 

improve model convergence and (or) stability. For solution type SLP, however, the 

response coefficients are calculated from a different set of base-condition withdrawal or 

injection flow rates at each iteration of the SLP algorithm. As described in the section on 

solution of nonlinear formulations, if the first base-condition run was successful but a later 

base run fails, it is possible to adjust the base vector of flow-rate decision variables (rather 

than terminate the algorithm) by use of a relaxation parameter (eq. 73), which is specified 

with input variable AFACT in the SOLN file. AFACT is restricted to a value of between 0 

and 1; a value of 0.5 is suggested.

The flow process also may fail during a perturbation run. This might occur because of 

a poor choice for the perturbation value calculated in equations 63 or 71. In this case, GWM 

can automatically adjust the perturbation value and re-execute the flow process. Two userspecified input variables control adjustment of the perturbation value if the flow process 

fails during a perturbation run. The first parameter, PGFACT, is the factor by which the 

perturbation value will be adjusted; PGFACT is restricted to values between 0 and 1. 

Depending on the type of failure, GWM will either increase or decrease the perturbation 

value. A suggested value of PGFACT is 0.5. Use of this value will result in a doubling or 

halving of the perturbation value at each subsequent perturbation attempt. The perturbation 

value for each flow-rate decision variable is adjusted independently, so that any adjustments made for one decision variable do not affect the perturbation value for other variables. The second parameter, NPGNMX, is the maximum number of perturbationadjustment attempts that will be made for a given flow-rate decision variable. If this value 

is exceeded for any of the decision variables, then the algorithm will terminate.

For nonlinear formulations in which the SLP algorithm is used, the simplex algorithm 

may determine that the linear program is infeasible or unbounded at any of the SLP iterations. This does not necessarily mean, however, that the original problem is either infeasible or unbounded. Rather, it may be that the particular set of response coefficients computed at that iteration produces an infeasible or unbounded constraint set. When a linear 

program fails during an SLP iteration, the RMS Package repeats the iteration, calculating 

the response coefficients with a different perturbation value. This is repeated up to 

NINFMX times, where NINFMX is user-specified. If no solution to each of the NINFMX 

linear programs can be found, then the RMS Package concludes that the original problem 

is infeasible or unbounded and halts the algorithm.

α
34 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

GWM File

The following items are read for each GWM run. The items must be listed in the order 

shown below. Values in bold are keywords that can be specified in either uppercase or lowercase letters. Input instructions for the files identified after each keyword are described in 

following subsections. Each GWM run requires specification of a DECVAR, OBJFNC, 

VARCON, and SOLN file type; the SUMCON, HEDCON, and STRMCON file types 

are optional. GWM automatically selects file units, so the user does not specify file units 

for these files.

Input items:

0. #Text

1. Each management problem must include a file that provides information about the 

decision variables and is read from the file whose name is specified by Fname:

DECVAR Fname

2. Each management problem must include a file that provides information about the 

objective function and is read from the file whose name is specified by Fname:

OBJFNC Fname

3. Each management problem must include a file that provides information on the lower 

and upper bounds specified for the flow-rate and external decision variables and is 

read from the file whose name is specified by Fname:

VARCON Fname

4. Each management problem may also include up to three additional files that provide 

information about the other types of constraints of the management problem that are 

allowed in GWM (summation constraints, head constraints, and streamflow constraints). These files are specified by the following records, which can be listed in any 

order:

SUMCON Fname

HEDCON Fname

STRMCON Fname

5. Each management problem must include a file that provides information about the 

solution and output-control parameters and is read from the file whose name is specified by Fname:

SOLN Fname

Explanation of variables:

Text—is a character variable (199 characters) that starts in column 2. Any characters can 

be included in Text. The “#” character must be in column 1. Lines beginning with “#” are 

restricted to these first lines of the file. Text is printed when the file is read. Item 0 can be 

repeated multiple times.

Fname—is a character variable that specifies the name of an existing file. Pathnames may 

be specified as part of Fname. 
Input Instructions and Output Files 35

Decision Variables (DECVAR) File

This file is used to define the decision variables of the management model. There are 

three types of decision variables. The primary decision variables are the flow rates (either 

withdrawal or injection) at each managed well site. A single well site may have more than 

one flow-rate decision variable associated with it; moreover, a single flow-rate decision 

variable can extend over one or more model cells and can be active during one or more 

stress periods. The second types of decision variables are external flow rates representing 

imported or exported water to the model domain. External variables do not have a direct 

effect on the system state variables and are not assigned to a specific location in the model. 

The third types of decision variables are binary variables used to define the status of each 

flow site or external variable as active (for example, the site is constructed) or inactive (the 

site is not constructed). Binary variables have a value of 0 (inactive site) or 1 (active site). 

One or more flow-rate and external decision variables are associated with each binary variable.

Input items:

0. #Text

1. IPRN

2. NFVAR NEVAR NBVAR

3a. The following records are read for each of NFVAR decision variables:

FVNAME NC LAY ROW COL FTYPE FSTAT WSP

3b. If NC > 1 in record 3a, then the following record is read NC times and the values of 

LAY, ROW, and COL read in record 3a are ignored:

RATIO LAY ROW COL

4. The following record is read for each of NEVAR external decision variables:

EVNAME ETYPE ESP

5. The following record is read for each of NBVAR binary decision variables:

BVNAME NDV BVLIST

Explanation of the variables:

Text—is a character variable (199 characters) that starts in column 2. Any characters can 

be included in Text. The “#” character must be in column 1. Lines beginning with “#” are 

restricted to these first lines of the file. Text is printed when the file is read. Item 0 can be 

repeated multiple times.

IPRN—is an integer variable that describes the type of input echo that is written to the 

GLOBAL file. IPRN must be specified as either 0 or 1. When IPRN equals 0, a minimum 

amount of information about the decision variables is written to the GLOBAL output file; 

when IPRN equals 1, detailed information about the decision variables is written to the 

GLOBAL output file.

NFVAR—is an integer variable equal to the number of flow-rate decision variables. 

NFVAR must be greater than 0. Only one flow-rate decision variable can be defined for a 

particular well site and set of stress periods, with the exception that both a withdrawal variable (FTYPE=W) and an injection variable (FTYPE=I) can be defined for the site.
36 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

NEVAR—is an integer variable equal to the number of external decision variables. NEVAR 

must be greater than or equal to 0. 

NBVAR—is an integer variable equal to the number of binary variables. NBVAR must be 

greater than or equal to 0. If NBVAR is 0, binary variables are not included in the management formulation.

FVNAME—is a character variable (maximum length of 10 characters) that is a unique 

name designated for the flow-rate decision variable. No spaces are allowed in the name. 

The end of the name is designated by a blank space.

NC—is an integer variable equal to the number of model cells over which the flow rate for 

decision-variable FVNAME is distributed. NC must be greater than or equal to 1. If NC 

equals 1, then all of the water withdrawn or injected at decision variable FVNAME is 

applied at the single model cell LAY, ROW, COL. If NC is greater than 1, then the flow rate 

calculated for decision variable FVNAME is distributed over the NC cells specified in 

record 3b.

LAY, ROW, COL—are integer variables equal to the layer, row, and column number of the 

model cell to which flow for decision-variable FVNAME will be assigned.

FTYPE—is a character variable that indicates whether the decision variable is a withdrawal 

or injection site. If FTYPE is W, the site is used for withdrawal; if FTYPE is I, the site is 

used for injection. If either withdrawal or injection is allowed at the site, two decision variables must be defined for the site, one for withdrawal (that is, with FTYPE = W) and one 

for injection (FTYPE = I).

FSTAT—is a character variable that indicates whether the decision variable will be considered in the management problem. If FSTAT is Y, the decision variable is available; if FSTAT 

is N, the decision variable is unavailable. If the decision variable is unavailable, then no 

withdrawal or injection will be calculated at the decision-variable location. For linear-optimization problems, FSTAT can be used to remove a well from the candidate set of decision 

variables without having to recalculate the response matrix (in this case, IRM = 0; see 

instructions for the Solution and Output-Control Parameters File below).

WSP—is a character string (up to 120 characters long) that indicates the stress periods 

associated with decision variable FVNAME. A single flow rate will be determined by 

GWM for all the stress periods included in WSP. The string must not contain any blank 

spaces. Multiple stress periods are listed using colons (:) or hyphens (-). For example,

1 indicates that stress period 1 is the only stress period associated with the decision 

variable;

1:3 indicates that the flow rate is the same for stress periods 1 and 3; and

1–12 indicates that the flow rate is the same for stress periods 1 through 12.

RATIO—is a real variable. RATIO is the fraction of the total flow rate for decision variable 

FVNAME that is distributed to cell LAY, ROW, COL. The sum of RATIO must equal 1.0 

for all of the NC cells specified for FVNAME; if the sum does not equal 1.0, GWM calculates the fraction for each cell by dividing the RATIO value specified for each cell by the 

sum of the RATIO values specified for all cells within FVNAME.
Input Instructions and Output Files 37

EVNAME—is a character variable (maximum length of 10 characters) that is a unique 

name designated for the external decision variable. No spaces are allowed in the name. The 

end of the name is designated by a blank space.

ETYPE—is a character variable that indicates whether the external variable is a source 

(import) or sink (export) of water. If ETYPE is IM, the variable is a source (import) of 

water; if ETYPE is EX, the variable is a sink (export) of water. Both types of external 

variables can be used in a management problem.

ESP—is a character string (up to 120 characters long) that indicates the stress periods associated with external variable EVNAME. A single flow rate will be determined by GWM for 

all the stress periods included in ESP. The string must not contain any blank spaces. Multiple stress periods are listed using colons (:) or hyphens (-). For example,

1 indicates that stress period 1 is the only stress period associated with the decision 

variable;

1:3 indicates that the flow rate for the external variable is the same for stress periods 

1 and 3; and

1-12 indicates that the flow rate for the external variable is the same for stress periods 

1 through 12.

BVNAME—is a character variable (maximum length of 10 characters) that is a unique 

name designated for the binary decision variable. The use of BVNAME, NDV, and 

BVLIST allows the user to associate one or more FVNAME or EVNAME decision variables with a single binary-variable identifier. For example, the user may want to define 12 

decision variables that are the monthly withdrawal rates at a single well site. If any one of 

the 12 decision variables is selected in the optimal solution, then an installation cost associated with the binary variable for the well site must be incurred. Also see example 3 in 

Appendix 1.

NDV—is an integer variable equal to the number of flow-rate or external decision variables 

associated with BVNAME.

BVLIST—is a list of the flow-rate and external decision variables associated with binary 

variable BVNAME. The list is drawn from the character names of these variables, 

FVNAME and EVNAME, defined in records 3a and 4. Each character variable in the list 

must be separated by a space, and there must be a total of NDV variables listed. The list can 

include any combination of decision variables, irrespective of well-site locations or stress 

period.

Objective Function (OBJFNC) File

This file is used to define the type of objective function that is to be solved and the 

coefficients for each decision variable in the objective function. Note that it is not necessary 

to include all decision variables in the objective function. In other words, some management formulations will have objective functions that do not include all of the decision 

variables defined for the problem.
38 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Input items:

0. #Text

1. IPRN

2. OBJTYP FNTYP

3. NFVOBJ NEVOBJ NBVOBJ

4. The following record is repeated for each of NFVOBJ flow-rate decision variables:

FVNAME FVOBJC

5. The following record is repeated for each of NEVOBJ external decision variables:

EVNAME EVOBJC

6. The following record is repeated for each of NBVOBJ binary decision variables:

BVNAME BVOBJC

Explanation of the variables:

Text—is a character variable (199 characters) that starts in column 2. Any characters can 

be included in Text. The “#” character must be in column 1. Lines beginning with “#” are 

restricted to these first lines of the file. Text is printed when the file is read. Item 0 can be 

repeated multiple times.

IPRN—is an integer variable that describes the type of input echo that is written to the 

GLOBAL file. IPRN must be specified as either 0 or 1. When IPRN equals 0, a minimum 

amount of information about the objective function is written to the GLOBAL output file; 

when IPRN equals 1, detailed information about the objective function is written to the 

GLOBAL output file.

OBJTYP—is a character variable used to define whether the objective is to maximize or 

minimize the objective function. OBJTYP must be defined as either MIN (for minimize) or 

MAX (for maximize).

FNTYP—is a character variable used to define the type of objective function. Currently, 

only one type of function is available, which is WSDV for weighted sum of decision variables. The user must specify WSDV on the input record.

NFVOBJ—is an integer variable equal to the number of flow-rate decision variables in the 

objective function and must have a value less than or equal to NFVAR specified in the decision-variables file.

NEVOBJ—is an integer variable equal to the number of external decision variables in the 

objective function and must have a value less than or equal to NEVAR specified in the decision-variables file.

NBVOBJ—is an integer variable equal to the number of binary decision variables in the 

objective function and must have a value less than or equal to NBVAR specified in the decision-variables file.

FVNAME—is a character variable (maximum length of 10 characters) that is one of the 

flow-rate decision-variable names. Each of the FVNAME variables listed must be defined 

in the DECVAR file. A flow-rate decision-variable name can only be listed once in the 

OBJFNC file.

FVOBJC—is a real variable that is a coefficient on each flow-rate decision variable 

FVNAME. For example, FVOBJC could represent the cost per unit volume of water withdrawn or injected at the management site.
Input Instructions and Output Files 39

EVNAME—is a character variable (maximum length of 10 characters) that is one of the 

external decision-variable names. Each of the EVNAME variables listed must be defined 

in the DECVAR file. An external decision-variable name can only be listed once in the 

OBJFNC file.

EVOBJC—is a real variable that is a coefficient on each external decision variable 

EVNAME. For example, EVOBJC could represent the cost per unit volume of water 

associated with the external variable.

BVNAME—is a character variable (maximum length of 10 characters) that is one of the 

binary-variable names. Each of the BVNAME variables listed must be defined in the 

DECVAR file. A binary-variable name can only be listed once in the OBJFNC file.

BVOBJC—is a real variable that is a coefficient on each binary variable BVNAME. For 

example, BVOBJC could represent the cost for installation of the management site. In 

typical usage of binary variables, the coefficients will be positive when OBJTYP is MIN 

and negative when OBJTYP is MAX. This will ensure that the binary variables are only 

active when their associated flow-rate and external decision variables are active.

Constraint Files

Four general types of constraints can be specified in a GWM management problem: 

constraints on the lower and upper bounds on the decision variables themselves by use of 

the VARCON file; linear-summation constraints by use of the SUMCON file; hydraulichead constraints by use of the HEDCON file; and streamflow constraints by use of the 

STRMCON file. A VARCON constraint file must be specified for each management 

problem, but the last three constraint type files are optional. Each of these four types 

of constraints is described below.

Decision-variable constraints (VARCON) file

The decision-variable constraints file is used to define lower and upper bounds for the 

flow-rate and external decision variables, and the reference flow rates to be used in the first 

ground-water flow run conducted by GWM. Records must be specified for all NFVAR and 

NEVAR decision variables defined in the DECVAR file. 

Input items:

0. #Text

1. IPRN

2. The following record is read for each of NFVAR decision variables:

FVNAME FVMIN FVMAX FVREF

3. The following record is read for each of NEVAR decision variables:

EVNAME EVMIN EVMAX

Explanation of the variables:

Text—is a character variable (199 characters) that starts in column 2. Any characters can 

be included in Text. The “#” character must be in column 1. Lines beginning with “#” are 

restricted to these first lines of the file. Text is printed when the file is read. Item 0 can be 

repeated multiple times.
40 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

IPRN—is an integer variable that describes the type of input echo that is written to the 

GLOBAL file. IPRN must be specified as either 0 or 1. When IPRN equals 0, a minimum 

amount of information about the decision-variable constraints is written to the GLOBAL 

output file; when IPRN equals 1, detailed information about the decision-variable 

constraints is written to the GLOBAL output file.

FVNAME—is a character variable (maximum length of 10 characters) that is one of the 

flow-rate decision-variable names. Each of the FVNAME variables listed must be defined 

in the DECVAR file. A flow-rate decision-variable name can only be listed once in the 

VARCON file.

FVMIN, FVMAX—are real variables that are equal to the minimum (FVMIN) and maximum (FVMAX) flow rate allowed for the decision variable. Values greater than or equal to 

0 must be specified for FVMIN and FVMAX; GWM will know whether the pumping rates 

are withdrawal or injection rates from the specification of FTYPE given in the Decision 

Variable File. FVMIN must be less than or equal to FVMAX. Note that a nonzero value of 

FVMIN implies that the decision variable has been associated with a binary variable in the 

DECVAR file. If the decision variable is not associated with a binary variable, then the nonzero value of FVMIN is ignored by GWM. The user can specify a nonzero lower 

bound for a flow-rate decision variable not associated with a binary variable by use of a 

linear-summation constraint (see description of SUMCON File).

FVREF—is a real variable equal to the flow rate for the decision variable that is used by 

GWM to calculate the reference values of the state variables. These include heads 

at drawdown-constraint locations if drawdown constraints are used (see input instructions 

for HEDCON file) and reference streamflows at streamflow-constraint locations if 

streamflow-depletion constraints are used (see input instructions for STRMCON file). 

FVREF also may be used to calculate base conditions for the calculation of the response 

matrix (see discussion of variable IBASE in the SOLN file). If no value is entered for 

FVREF, it is assigned a value of 0.

EVNAME—is a character variable (maximum length of 10 characters) that is one of the 

external decision-variable names. Each of the EVNAME variables listed must be defined 

in the DECVAR file. An external decision-variable name can only be listed once in the 

VARCON file.

EVMIN, EVMAX—are real variables that are equal to the minimum (EVMIN) and maximum (EVMAX) flow rate allowed for the external decision variable. Values greater than or 

equal to 0 must be specified for EVMIN and EVMAX; GWM will know whether the flow 

rates for the external variable are imported or exported flow rates from the specification of 

ETYPE given in the Decision Variable File. EVMIN must be less than or equal to EVMAX. 

Note that a nonzero value of EVMIN implies that the decision variable has been associated 

with a binary variable in the DECVAR file. If the decision variable is not associated with a 

binary variable, then the nonzero value of EVMIN is ignored by GWM. The user can 

specify a nonzero lower bound for an external decision variable not associated with a binary 

variable by use of a summation constraint (see description of SUMCON File).

Linear-summation constraints (SUMCON) file

The linear-summation constraints file is used to define linear relations among decision 

variables. Examples of the use of these constraints are given by equations 17–28 in the 

“Formulation of Ground-Water Management Problems” section of the report, and in the 

sample problems. 
Input Instructions and Output Files 41

Input items:

0. #Text

1. IPRN

2. SMCNUM

3a. Records 3a and 3b are read for each of the SMCNUM constraints:

SMCNAME NTERMS TYPE RHS

3b. The following record is repeated NTERMS times for each of the NTERMS specified 

in record 3a:

GVNAME GVCOEFF

Explanation of the variables:

Text—is a character variable (199 characters) that starts in column 2. Any characters can 

be included in Text. The “#” character must be in column 1. Lines beginning with “#” are 

restricted to these first lines of the file. Text is printed when the file is read. Item 0 can be 

repeated multiple times.

IPRN—is an integer variable that describes the type of input echo that is written to the 

GLOBAL file. IPRN must be specified as either 0 or 1. When IPRN equals 0, a minimum 

amount of information about the summation constraints is written to the GLOBAL output 

file; when IPRN equals 1, detailed information about the summation constraints is written 

to the GLOBAL output file.

SMCNUM—is an integer variable equal to the number of summation constraints defined 

in the file.

SMCNAME—is a character variable (maximum length of 10 characters) that is a unique 

name designated for the constraint. No spaces are allowed in the name. The end of the name 

is designated by a blank space.

NTERMS—is an integer variable equal to the number of terms on the left-hand side of the 

constraint. All of the terms are added together to form the left-hand side of the constraint.

TYPE—is a character variable used to specify the type of constraint. Three options are 

allowed:

LE: the left-hand side of the equation is less than or equal to the right-hand side of the 

constraint;

GE: the left-hand side of the equation is greater than or equal to the right-hand side of 

the constraint;

EQ: the left-hand and right-hand sides of the constraint are equal.

RHS—is a real variable equal to the value of the right-hand side of the constraint.

GVNAME—is a character variable (maximum length of 10 characters) that is one of the 

decision-variable names defined in the DECVAR file for either an FVNAME, EVNAME, 

or BVNAME variable. Any combination of flow-rate, external, and binary decision variables may be present in a constraint. The user must ensure that the variables used are 

logically consistent.

GVCOEFF—is a real variable equal to the value of the coefficient in front of variable 

GVNAME. The user must ensure that a consistent set of units is used for all GVCOEFF and 

RHS terms.
42 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Head constraints (HEDCON) file

The head constraints file is used to define head constraints at model cells. These 

include upper and lower bounds on heads, drawdowns, head differences between two cells, 

and gradients between two cells.

Input items:

0. #Text

1. IPRN

2. NHB NDD NDF NGD

3. The following record is read for each of the NHB constraints:

HBNAME LAYH ROWH COLH TYPH BND NSP

4. The following record is read for each of the NDD constraints:

DDNAME LAYD ROWD COLD TYPD BND NSP

5. The following record is read for each of the NDF constraints:

HDIFNAME LAY1 ROW1 COL1 LAY2 ROW2 COL2 HD NSP

6. The following record is read for each of the NGD constraints:

GRADNAME LAY1 ROW1 COL1 LAY2 ROW2 COL2 LEN GRAD NSP

Explanation of the variables:

Text—is a character variable (199 characters) that starts in column 2. Any characters can 

be included in Text. The “#” character must be in column 1. Lines beginning with “#” are 

restricted to these first lines of the file. Text is printed when the file is read. Item 0 can be 

repeated multiple times.

IPRN—is an integer variable that describes the type of input echo that is written to the 

GLOBAL file. IPRN must be specified as either 0 or 1. When IPRN equals 0, a minimum 

amount of information about the head constraints is written to the GLOBAL output file; 

when IPRN equals 1, detailed information about the head constraints is written to the 

GLOBAL output file.

NHB—is an integer variable equal to the number of head-bound constraints that need to be 

satisfied in the management model.

NDD—is an integer variable equal to the number of drawdown constraints that need to be 

satisfied in the management model.

NDF—is an integer variable equal to the number of head difference constraints that need 

to be satisfied in the management model.

NGD—is an integer variable equal to the number of gradient constraints that need to be 

satisfied in the management model.

Head-bound constraints:

HBNAME—is a character variable (maximum length of 10 characters) that is a unique 

name designated for the head-bound constraint. No spaces are allowed in the name. The end 

of the name is designated by a blank space.

LAYH, ROWH, COLH—are integer variables equal to the layer, row, and column number 

of the model cell in which the head-bound constraint is located.
Input Instructions and Output Files 43

TYPH—is a character variable used to specify the type of head bound. Two options are 

allowed:

LE: head calculated by the model must be less than or equal to the value specified by 

BND;

GE: head calculated by the model must be greater than or equal to the value specified 

by BND.

BND—is a real variable equal to the specified upper ( ) or lower ( ) bound on 

head at the model cell at the end of the stress period (fig. 2A).

NSP—is an integer variable that indicates the stress period during which the constraint is 

imposed. If the constraint is imposed over multiple stress periods, then a separate record 

must be provided for each stress period.

Drawdown constraints:

DDNAME—is a character variable (maximum length of 10 characters) that is a unique 

name designated for the drawdown constraint. No spaces are allowed in the name. The end 

of the name is designated by a blank space.

LAYD, ROWD, COLD—are integer variables equal to the layer, row, and column number 

of the model cell in which the drawdown constraint is located.

TYPD—is a character variable used to specify the type of head-drawdown bound. Two 

options are allowed:

LE: drawdown calculated by the model at the model cell must be less than or equal to 

the value specified by BND;

GE: drawdown calculated by the model must be greater than or equal to the value 

specified by BND.

BND—is a real variable equal to the specified upper ( ) or lower ( ) bound on 

drawdown at the model cell at the end of the stress period (fig. 2B).

NSP—as defined for Record 3.

Head-difference constraints:

HDIFNAME—is a character variable (maximum length of 10 characters) that is a unique 

name designated for the head-difference constraint. No spaces are allowed in the name. The 

end of the name is designated by a blank space.

LAY1, ROW1, COL1—are integer variables equal to the layer, row, and column number of 

the model cell corresponding to the first location, , in the head-difference constraint 

(fig. 2C). 

LAY2, ROW2, COL2—are integer variables equal to the layer, row, and column 

number of the model cell corresponding to the second location, , in the headdifference constraint (fig. 2C). GWM requires that the head at the second head-difference 

location be lower than the head at the first location by an amount of at least HD.

HD—is a real variable equal to the specified difference in heads, , between the 

first and second model cells at the end of the stress period (fig. 2C).

NSP—as defined for Record 3.

hijkt ,, ,

u hijkt ,, ,l

dd ijkt ,, ,

u dd ijkt ,, ,l

( ) ijk , , 1

( ) ijk , , 2

hdifijkt ,, , ( )1 2,
44 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Gradient constraints:

GRADNAME—is a character variable (maximum length of 10 characters) that is a unique 

name designated for the gradient constraint. No spaces are allowed in the name. The end of 

the name is designated by a blank space.

LAY1, ROW1, COL1—are integer variables equal to the layer, row, and column 

number of the model cell corresponding to the first location, , in the gradient 

constraint (fig. 2D).

LAY2, ROW2, COL2—are integer variables equal to the layer, row, and column number 

of the model cell corresponding to the second location, , in the gradient constraint 

(fig. 2D). GWM requires that the head at the second head-difference location be lower than 

the head at the first location.

LEN—is a real variable equal to the distance between the first and second model cells (

shown in fig. 2D).

GRAD—is a real variable equal to the specified gradient, , between the first 

and second model cells at the end of the stress period (fig. 2D).

NSP—as defined for Record 3.

Streamflow constraints (STRMCON) file

The streamflow constraints file is used to define streamflow constraints when the 

STR1 Package of MODFLOW (Prudic, 1989) is used. Two types of streamflow constraints 

are allowed—constraints on the upper and lower bounds on streamflow and constraints on 

the upper and lower bounds on streamflow depletion. 

Input items:

0. #Text

1. IPRN

2. NSF NSD

3. The following record is read for each of the NSF constraints:

SFNAME SEG REACH TYPSF BND NSP

4. The following record is read for each of the NSD constraints:

SDNAME SEG REACH TYPSD BND NSP

Explanation of the variables:

Text—is a character variable (199 characters) that starts in column 2. Any characters can 

be included in Text. The “#” character must be in column 1. Lines beginning with “#” are 

restricted to these first lines of the file. Text is printed when the file is read. Item 0 can be 

repeated multiple times.

IPRN—is an integer variable that describes the type of input echo that is written to the 

GLOBAL file. IPRN must be specified as either 0 or 1. When IPRN equals 0, a minimum 

amount of information about the streamflow constraints is written to the GLOBAL output 

file; when IPRN equals 1, detailed information about the streamflow constraints is written 

to the GLOBAL output file.

( ) ijk , , 1

( ) ijk , , 2

Δx

gradijkt ,, , ( )1 2,
Input Instructions and Output Files 45

NSF—is an integer variable equal to the number of streamflow constraints that need to be 

satisfied in the management model.

NSD—is an integer variable equal to the number of streamflow-depletion constraints that 

need to be satisfied in the management model.

Streamflow constraints:

SFNAME—is a character variable (maximum length of 10 characters) that is a unique name 

designated for the streamflow constraint. No spaces are allowed in the name. The end of the 

name is designated by a blank space.

SEG, REACH—are integer variables equal to the segment and reach number, as specified 

in the STR1 Package, of the model cell in which the streamflow constraint is located 

(Prudic, 1989).

TYPSF—is a character variable used to specify the type of streamflow constraint. Two 

options are allowed:

LE: streamflow at the stream site must be less than or equal to the value specified by 

BND;

GE: streamflow at the stream site must be greater than or equal to the value specified 

by BND.

BND—is a real variable equal to the specified amount of streamflow allowed at the stream 

site at the end of the stress period.

NSP—is an integer variable that indicates the stress period during which the constraint is 

imposed. To impose the constraint for multiple stress periods, define additional constraints.

Streamflow-depletion constraints:

SDNAME—is a character variable (maximum length of 10 characters) that is a unique 

name designated for the streamflow-depletion constraint. No spaces are allowed in the 

name. The end of the name is designated by a blank space.

SEG, REACH—are integer variables equal to the segment and reach number of the model 

cell in which the streamflow-depletion constraint is located (Prudic, 1989).

TYPSD—is a character variable used to specify the type of streamflow-depletion 

constraint. Two options are allowed:

LE: streamflow depletion at the stream site must be less than or equal to the value 

specified by BND;

GE: streamflow depletion at the stream site must be greater than or equal to the value 

specified by BND.

BND—is a real variable equal to the specified amount of streamflow depletion allowed at 

the stream site at the end of the stress period.

NSP—as defined for Record 3.
46 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Solution and Output-Control Parameters (SOLN) File

The solution and output-control parameters file is used to define several variables that 

control the solution algorithm for the optimization problem and the type and amount of output that is printed to the output files. 

Input items:

0. #Text

1. SOLNTYP

If SOLNTYP is NS then:

2a. DELTA

2b. NSIGDIG NPGNMX PGFACT

2c. RMNAME

Skip to record 6a.

If SOLNTYP is MPS then:

3a. DELTA

3b. NSIGDIG NPGNMX PGFACT

3c. MPSNAME

Skip to record 6a.

If SOLNTYP is LP then:

4a. IRM

4b. LPITMAX BBITMAX

4c. DELTA

4d. NSIGDIG NPGNMX PGFACT

4e. BBITPRT RANGE

The following record is read if IRM equals 0 or 1:

4f. RMNAME

Skip to record 6a.

If SOLNTYP is SLP then:

5a. SLPITMAX LPITMAX BBITMAX

5b. SLPVCRIT SLPZCRIT DINIT DMIN DSC

5c. NSIGDIG NPGNMX PGFACT AFACT NINFMX

5d. SLPITPRT BBITPRT RANGE

6a. IBASE

If IBASE equals 1, the following record is read for each of NFVAR decision variables:

6b. FVNAME FVBASE
Input Instructions and Output Files 47

Explanation of the variables:

Text—is a character variable (199 characters) that starts in column 2. Any characters can 

be included in Text. The “#” character must be in column 1. Lines beginning with “#” are 

restricted to these first lines of the file. Text is printed when the file is read. Item 0 can be 

repeated multiple times.

SOLNTYP—is a character variable. SOLNTYP must be specified as either NS, MPS, LP, 

or SLP:

NS: no solution to the management formulation will be found. GWM will calculate the 

response matrix and write it to the response-matrix file specified by RMNAME in record 

2c; GWM will then stop.

MPS: no solution to the management formulation will be found. GWM will write the 

management formulation in MPS (Mathematical Programming System) format to file 

specified by MPSNAME in record 3c; GWM will then stop.

LP (with or without binary variables): the optimization formulation is solved by using 

linear programming, or, for problems with binary variables, linear programming and the 

branch and bound method. A SOLNTYP equal to LP is normally used if the flow model 

contains only linear features (see “Solution of Ground-Water Management Problems 

with GWM” for a discussion of linear and nonlinear features). Alternatively, the use of 

LP allows the user to force the problem to be solved as a linear problem even if nonlinear 

features are present. In that case, the management formulation will be solved with a 

single response matrix. This option should be used carefully and may lead to inaccurate 

results if the management problem has a significant nonlinear response.

SLP (with or without binary variables): the optimization formulation is solved using 

sequential (iterative) linear programming, or, for problems with binary variables, 

sequential linear programming and the branch and bound method.

If SOLNTYP is NS:

DELTA—is a real variable equal to the perturbation parameter (eq. 63) used to determine 

the response matrix. DELTA is multiplied by FVMAX for each flow-rate decision variable 

to determine each perturbation value. A positive value of DELTA implies a forwarddifference calculation of the response coefficient (that is, an increase in flow rate), whereas 

a negative value implies a backward-difference calculation (that is, a decrease in flow rate).

NSIGDIG—is an integer variable equal to a lower limit on the number of significant digits 

in response-matrix entries. For each entry in a column of the response matrix, the ratio of 

the difference in observed state to the HCLOSE variable is computed. If the largest ratio in 

the column has fewer than NSIGDIG significant digits, the perturbation is considered to 

have failed.

NPGNMX—is an integer variable equal to the maximum number of attempts to achieve a 

successful flow-process run. Failure may occur during either base or perturbation flowprocess runs. When failure is detected, automatic resetting of flow-rate decision-variable 

values may produce a successful solution. NPGNMX controls the maximum number of 

attempts for perturbation adjustments (controlled by PGFACT) or for base adjustments 

(controlled by AFACT). 

δ

0
48 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

PGFACT—is a real variable equal to the perturbation step-length adjustment factor used 

during perturbation failure. PGFACT must be greater than 0 and less than 1. A value of 0.5 

is suggested. 

RMNAME—is the file name (or pathname) to which the response matrix will be written.

If SOLNTYP is MPS:

DELTA—as defined for item 2a.

NSIGDIG, NPGNMX, PGFACT—as defined for item 2b.

MPSNAME—is the file name (or pathname) to which the the formulation will be written 

in MPS format.

If SOLNTYP is LP:

IRM—is an integer variable equal to 0, 1, or 2. Its value specifies whether or not the 

response matrix will be calculated or read from an input file and whether or not the response 

matrix will be saved. A value of 0 indicates that the response matrix has been generated in 

a previous run and will be read from the file specified by RMNAME in record 4f; this feature can greatly reduce the time required to solve GWM if the response matrix has been 

previously calculated. A value of 1 indicates that the response matrix will be calculated by 

GWM and written to the file specified by RMNAME in record 4f. A value of 2 indicates 

that the response matrix will be calculated by GWM but not written to a file.

LPITMAX—is an integer variable. Its value is the maximum number of iterations 

allowed for the linear program solver; this limit prevents the solver from iterating indefinitely if it does not converge to a solution. If the linear solver is being used and the value 

of LPITMAX is reached, the program will be terminated, and the output file will indicate 

that the maximum number of iterations has been exceeded. A typical value for LPITMAX 

is ten times the number of constraints.

BBITMAX—is an integer variable. BBITMAX is only relevant if the management problem contains one or more binary variables, otherwise, its value is ignored. BBITMAX is 

the maximum number of iterations allowed for the branch and bound program solver. 

Each iteration consists of one solution of the linear program. If the value of BBITMAX is 

reached, the program will be terminated, and the output file will indicate that the maximum 

number of iterations has been exceeded. 

DELTA—as defined for item 2a.

NSIGDIG, NPGNMX, PGFACT—as defined for item 2b.

BBITPRT—is an integer variable that specifies whether output describing the details of the 

branch and bound algorithm for solving mixed binary problems will be written to the 

GLOBAL output file. A value of 1 indicates that this output will be created and a value of 

0 indicates that it will not. For problems with many binary variables, this file can be very 

large. If the management problem contains no binary variables, this value will be ignored.

RANGE—is an integer variable that indicates the status of the range analysis. A value of 1 

indicates that range analysis should be done and the results written to the GLOBAL file. A 

value of 0 indicates no range analysis. Range analysis is described in Appendix 2. Range 

analysis is based on the assumption that the optimization problem is strictly linear with continuous variables. If binary variables or nonlinear responses are significant in the problem, 

then the range analysis may be inaccurate.
Input Instructions and Output Files 49

RMNAME—is the file name (or pathname) from which the response matrix will be read if 

IRM equals 0 and is the file name (or pathname) to which the response matrix will be 

written if IRM equals 1.

If SOLNTYP is SLP:

SLPITMAX—is an integer variable. Its value is the maximum number of iterations allowed 

for the sequential linear-programming algorithm. If the value of SLPITMAX is reached, the 

program will be terminated, and the output file will indicate that the maximum number of 

iterations has been exceeded.

LPITMAX—as defined for item 4b.

BBITMAX—as defined for item 4b.

SLPVCRIT—is a real variable. Its value is the convergence criterion (eq. 69A), which 

is the first of two termination rules when the sequential linear programming algorithm is 

used. This rule is satisfied when the change in the values of all flow-rate decision variables 

from the previous iteration to the current iteration is less than a fraction, , of the magnitude of the flow-rate decision variables at the current iteration.

SLPZCRIT—is a real variable. Its value is the convergence criterion (eq. 69B), which is 

the second of two termination rules when the sequential linear programming algorithm is 

used. This rule is satisfied when the change in the value of the objective function, Z, is less 

than a specified fraction, , of the value of the objective function.

DINIT, DMIN, DSC—are real variables that control the value of the perturbation variable 

(eq. 72) used to compute response coefficients. DINIT is the perturbation variable used for 

the first iteration, DMIN is the minimum perturbation variable used, and DSC is a parameter that controls the rate of change of the perturbation parameter. DINIT and DMIN must 

have the same sign. Positive values of DINIT and DMIN imply a forward-difference calculation of the response coefficient (that is, an increase in flow rate), whereas negative values 

imply a backward-difference calculation (that is, a decrease in flow rate). DSC must always 

be positive.

NSIGDIG, NPGNMX, PGFACT—as defined for item 2b.

AFACT—is a real variable equal to the relaxation parameter ( in eq. 73) used to determine 

a temporary base solution when a base run fails. AFACT controls the interpolation between 

the current base solution and the most recent successful base solution. AFACT must be 

greater than 0 and less than 1. A value close to 0 implies that the temporary base solution 

will be close to the current base solution, whereas a value close to 1.0 implies that the temporary base solution will be close to the prior base solution. A value of 0.5 is suggested.

NINFMX—is an integer variable that specifies the maximum number of consecutive infeasible iterations that will be accepted by the SLP algorithm before the algorithm terminates.

SLPITPRT—is an integer variable that specifies whether output describing the details of 

the sequential-iteration algorithm will be written to the GLOBAL output file. A value of 1 

indicates that this output will be created and a value of 0 indicates that it will not. 

BBITPRT—as defined in item 4e.

RANGE—as defined in item 4e.

ε1

ε1

ε2

ε2

α
50 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

IBASE—is an integer variable equal to 0 or 1 that indicates the source for the values of the 

flow-rate decision variables that will be used as the base run. For problems solved using the 

SLP algorithm, these values are the starting point for the iterative algorithm. A value of 

IBASE equal to 0 indicates that the reference flow rates (FVREF) specified for each flowrate decision variable in file VARCON will be used in the base run (and record 6b is not 

necessary). A value of IBASE equal to 1 indicates that the flow rates specified for each 

decision variable by FVBASE in record 6b will be used to calculate the base run.

FVNAME—is a character variable (maximum length of 10 characters) that is one of the 

flow-rate decision-variable names. Each of the FVNAME variables listed must be defined 

in the DECVAR file. A flow-rate decision-variable name can only be listed once in the 

SOLN file.

FVBASE—is a real variable equal to the rate for the flow-rate decision variable. These 

values are used by GWM to calculate the base run. If the SLP solution algorithm is used, 

these values are the starting point for the iterative algorithm.

Output Files

Two output files are always produced by a GWM run—a GLOBAL file and a LIST 

file. Two additional files, one to hold response matrixes and one to hold formulations written in MPS format, also are produced if specified by the user (see instructions for variables 

RMNAME and MPSNAME in SOLN file). The GLOBAL and LIST files are used in a 

manner that is similar to their use in MODFLOW applications that include the parameterestimation or sensitivity processes (Harbaugh and others, 2000, p. 6–7). Most of the output 

from a particular run of the GWF Process is written to the LIST file. Input data that are written to the LIST file include array-allocation information, hydraulic properties of simulated 

aquifers, and initial and boundary conditions. Output that results from the GWF Process 

includes calculated heads, drawdowns, and volumetric-budget terms. The LIST file is 

erased and generated anew each time a ground-water flow run is required.

The GLOBAL file contains information that applies to the GWM run as a whole. The 

first information that is written includes the names and types of MODFLOW files that are 

opened from the NAME file, information about the spatial and temporal discretization of 

the model that is read from the MODFLOW discretization (DIS) file, and input for the 

selected MODFLOW solver package. The remainder of the GLOBAL file consists of input 

to, information about, and results of the GWM run. The first part of the GWM section consists of data that are echoed from the input files; the amount of input data that is printed 

depends on the values of IPRN specified by the user in the GWM input files. Information 

about the size of the management problem, which consists of the number of decision 

variables and the number of constraints, also is printed.

The next part of the GWM section ("Solution Algorithm") provides information about 

the solution process for the management problem. First, results of the call to the GWF Process for the base-condition run are reported, including the status of each of the constraints. 

The status of each constraint is indicated as either "satisfied," "not met," or "near-binding." 

For output purposes, GWM considers a constraint to be near-binding if the left and right 

sides of the constraint are in agreement to at least five digits. Next, information about the 

flow-process runs that are required to generate the response matrix is reported, including 

the perturbation value used at each candidate well and information about the number of significant digits in the response matrix. If the problem is nonlinear or includes binary variables (and if either of the input variables SLPITPRT or BBITPRT has been set to 1), then 

this section of the output also will provide detailed information on the progress of iterations 

necessary to solve the problem. If the problem contains binary variables, then the output 
Input Instructions and Output Files 51

will include the solution to the management problem calculated for each iteration (subproblem) of the branch and bound algorithm. The solutions are designated as relaxed, current 

best solution, feasible, or infeasible. For each solution, the name of the flow-rate or external 

decision variable and its value for the current solution are reported, as well as the status of 

each binary variable.

The last part of the GWM section ("Ground-Water Management Solution") provides 

information on the solution of the management problem, including the value of the objective function, the optimal values for each of the decision variables, and the binding constraints. Although linear programs often have many constraints, only a subset of them typically control the optimal solution. Those constraints that restrict the value of the objective 

function are said to be binding, because they prevent the decision variables from taking on 

values that further improve the objective function, and, therefore, bind the solution. At the 

optimal solution, the inequality constraints that are met as equalities are binding. Conversely, nonbinding constraints do not affect the optimal values of the decision variables 

and could be removed from the formulation without changing the solution.

Each constraint has associated with it a shadow price (or dual variable) that provides 

information on the effect of small changes in the value of the right-hand-side constraint 

coefficient on the value of the objective function. Shadow prices represent the marginal 

increase in benefit that results from relaxing each constraint coefficient. When the righthand-side value of a binding constraint increases by a unit amount, the objective function 

will change by an amount given by the shadow price. This implies that the optimal solution 

is quite sensitive to constraints that have large shadow prices. Conversely, non-binding 

constraints have zero shadow prices, because small changes in the right-hand side of a nonbinding constraint will have no effect on the objective function. Shadow prices are reported 

in the GLOBAL file along with the list of binding constraints.

Shadow prices represent the local sensitivity of the optimal solution to the right-hand 

side of the constraints and are valid as long as changes in the constraints do not alter the 

basis at the optimal solution. For management problems solved by the RMS Package of 

GWM, the basis consists of the set of flow-rate and external decision variables that have 

values between their upper and lower bounds at the optimal solution, as well as the set of 

slack variables associated with the nonbinding constraints. Range analysis is used to determine the range of values over which a shadow price is valid. The range analysis is reported 

after the list of binding constraints if the input variable RANGE in the SOLN file is set to 

1. The range analysis calculates the interval over which the right-hand-side value of the 

constraint can vary without changing the basis of the optimal solution, if all other parameters of the management model remain unchanged. The analysis also determines the decision variables and (or) constraints that enter and leave the basis when the right-hand-side 

value is increased or decreased beyond its range.

A range analysis on the objective-function coefficients is also done and reported if 

requested by the user. In this analysis, the interval over which the coefficients can vary 

without changing the basis of the optimal solution is determined. Although the optimal 

basis remains unchanged, if the objective-function coefficients stay within their specified 

ranges, the values of the objective function and shadow prices are functions of the coefficients and therefore will be affected by any changes. The coefficient range analysis also 

determines the decision variables and (or) constraints that enter and leave the basis when 

the objective-function coefficient of interest is increased or decreased beyond its range. In 

addition to reporting the objective-function coefficient range-analysis information, the 

reduced cost associated with each decision variable is provided. The reduced cost of each 

decision variable that is not part of the basis at the optimal solution (nonbasic variables) is 

the amount by which the value of the objective function would be penalized if that decision 
52 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

variable were to become basic. For example, if a nonbasic decision variable has a value of 

zero and were to be brought into the optimal solution at a nonzero value, the objective function would be worsened—increased for a minimization problem and decreased for a maximization problem—by an amount equal to the reduced cost of the variable. The reduced 

cost of each of the decision variables in the optimal basis is zero.

Because range analysis assumes continuous variables and linear constraints, it is 

strictly valid for linear, non-binary problems only. Range analysis may be reasonably accurate for mildly nonlinear problems. If binary variables are present, then range analysis can 

be reported; however, it should be used with caution. The reported range analysis is based 

on the final linear program solved during the branch and bound algorithm and does not 

include the changes in binary variables that might result when right-hand-side and objective-function coefficients are changed. The most complete range analysis, when binary 

variables are present, can be obtained by resolving the problem as a non-binary, linear program with the inactive flow variables removed from the problem.

Additional details concerning sensitivity and range analyses can be found in Appendix 

2 and Ahlfeld and Mulligan (2000).

The last step in a successful GWM run is a final run of the GWF Process using the optimal flow rates determined by the solution algorithm. The LIST file will contain output from 

the GWF Process using these flow rates. The GLOBAL output file will contain the final 

status of the constraints. A successful run should indicate that all constraints are either satisfied or near-binding. A distinction is made between binding constraints and near-binding 

constraints. When solving the linear program, which contains an approximation to the 

GWF Process, certain constraints will be binding, that is, satisfied as strict equalities. However, when the GWF Process is run with the optimal flow rates, these same constraints may 

not be exactly binding. This can result from nonlinear responses in the GWF Process and 

precision limitations in the computation of heads. In addition, some constraints that are not 

strictly binding in the linear program, may be very close to binding in the GWF Process. 

As a result, there may be a difference in the set of constraints listed as binding in the linear 

program output and those listed as near-binding in the final GWF Process output.

Sample Problems

Three sample problems are provided to demonstrate how to set up files for GWM, 

some of the output generated by GWM, and how GWM may be used to solve typical 

ground-water management problems. It is assumed that the reader is familiar with the 

MODFLOW code (Harbaugh and others, 2000), including terminology that is used in 

MODFLOW. Selected input and output files are listed at the end of each sample problem.

Sample Problem 1: DEWATER

This sample problem represents a steady-state dewatering problem for a construction 

site. The objective of the ground-water management problem is to minimize the cost of 

withdrawing ground water to lower heads to an elevation of 50 ft so that footings can be 

installed in the area shown in figure 3. Wells at the construction site will be pumped for 

1,000 days. The aquifer at the site is confined and is simulated by a single model layer that 

is 3,000 ft long and 2,000 ft wide. The model grid consists of 20 rows and 30 columns, and 

each grid cell is 100 ft by 100 ft (fig. 3). The model uses no-flow boundary conditions along 

the north and south boundaries of the aquifer and constant heads of 80 ft and 60 ft along the 

east and west boundaries of the aquifer, respectively. The transmissivity of the aquifer is 

50 ft2/d.
Sample Problems 53

N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

1

20

10

2

3

4

5

6

7

8

9

11

12

13

14

15

16

17

18

19

CANDIDATE WITHDRAWAL WELL

HEAD CONSTRAINT LOCATION

CONSTRUCTION AREA

* * *

* *

* *

*

**

*

1

3

5

2

4

6

7

5

EXPLANATION

COLUMNS

ROWS

Constant head of 80 feetConstant head of 60 feet

Figure 3. Model grid for DEWATER sample problem.

The MODFLOW input files consist of a NAME file, a DIS-Package file, a BAS6-

Package file, a BCF6-Package file, and a PCG-Package file. In the DIS file, the number of 

stress periods (NPER) is set to 1, the time units are specified as days (ITMUNI = 4), the 

length of the single stress period is 1,000 days (PERLEN = 1.000E+3), and a single time 

step is used for the stress period (NSTP = 1). A very small head-change convergence 

criterion of 1.0E-8 ft is specified in the PCG solution package. 

The management problem is solved in two ways. In the first formulation, only the 

operational costs of pumping are considered; this formulation results in a linear problem. 

The GWM input files that are necessary for the linear formulation are: DECVAR, 

OBJFNC, VARCON, HEDCON, and SOLN. In the second formulation, both the operational and construction costs of the wells are considered; in addition, there is a restriction 

on the minimum number of wells that must be constructed. The construction costs and constraint on the minimum number of wells that must be constructed add binary variables to 

the problem; the formulation is therefore a mixed-binary linear problem. The GWM input 

files that are necessary for the mixed-binary linear problem include all of those necessary 

for the linear formulation, as well as a SUMCON file to specify the minimum number of 

wells that must be constructed. Input files for these two formulations are listed at the end 

of the sample problem.
54 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Linear Formulation

Seven candidate well locations are selected as possible locations of withdrawal. These 

seven flow-rate decision variables (named Q1, Q2, and so forth) are specified in the DECVAR file with input variable NFVAR, and each decision variable extends over a single cell 

(NC=1). The minimum and maximum pumping rates at each well, which are specified in 

the VARCON file, are 0 and 20,000 ft3/d, respectively. 

The linear problem is solved using two approaches. In the first approach, the operational costs of pumping are considered to be directly proportional to the amount of water 

withdrawn (in units of cubic feet per day), so the objective is to minimize total withdrawal 

from the seven wells, with unitless objective-function coefficients of 1.0:

Minimize . (74)

Note that equation 74 is multiplied by the length of the stress period (1,000 days) by GWM 

to determine the value of the objective function in units of cubic feet.

The 50-ft head criterion can be imposed by defining 10 upper-bound head constraints 

(constraint type NHB) in the HEDCON file at the cell locations shown in figure 3. The 

problem solution type specified in the SOLN file is LP (linear program). A perturbationparameter (DELTA) value of 0.5 is specified. DELTA is multiplied by the maximum withdrawal rate at each well (input variable FVMAX in file VARCON) to determine the initial 

perturbation rate at each well. In this example, because the maximum withdrawal rate for 

each well is 20,000 ft3/day, the initial perturbation value for each well is -10,000 ft3/d. 

Variable IBASE in the SOLN file is specified as 0, which means that the base-condition 

withdrawal rates for each well are specified by FVREF in the VARCON file; these rates 

are specified as 0 ft3/d for each well.

The value of the objective function at the optimal solution for this approach is 

2.8657x106 ft3 of water withdrawn. Four wells were selected for pumping in the optimal 

solution: well Q1 pumps at a rate of 1,077 ft3/d, well Q2 at 78.2 ft3/d, well Q4 at 769.0 ft3/d, 

and well Q7 at 941.1 ft3/d. Wells Q3, Q5, and Q6 are inactive in the optimal solution (that 

is, they have pumping rates of 0 ft3/d). Four of the head constraints are binding at the optimal solution, those located in cells (1, 6, 13), (1, 6, 17), (1, 10, 13), and (1, 14, 17). The 

GLOBAL output file for this GWM run is listed at the end of the sample problem.

In the second approach, the actual costs of pumping during construction are included. 

These costs, which are independent of the depth to water (that is, the lift), are $0.02/ft3 of 

water withdrawn:

Minimize 

.

(75)

Again, GWM multiplies each term in the objective function by the length of the stress 

period (1,000 days).

Q1 ++++++ Q2 Q3 Q4 Q5 Q6 Q7

$0.02

ft3 ------------ ⎝ ⎠ ⎛ ⎞Q1 $0.02

ft3 ------------ ⎝ ⎠ ⎛ ⎞Q2 $0.02ft3 ------------ ⎝ ⎠ ⎛ ⎞Q3 $0.02ft3 ------------ ⎝ ⎠ ⎛ ⎞Q4 

$0.02

ft3 ------------ ⎝ ⎠ ⎛ ⎞Q5 $0.02

ft3 ------------ ⎝ ⎠ ⎛ ⎞Q6 $0.02ft3 ------------ ⎝ ⎠ ⎛ ⎞Q7

++++

+ +
Sample Problems 55

The value of the objective function at the optimal solution for the second approach is 

$57,313 (output from this formulation is not included in this report). The same four wells 

are selected for pumping as were selected in the first approach, and the pumping rates at 

each of the wells are the same as for the first approach. The reason that the solutions for the 

two approaches are the same is that the solution to minimization (or maximization) problems that differ only by a constant multiple in the objective function will be the same. In 

this case, the two approaches differ only by the multiple $0.02/ft3 in the objective function, 

so the solutions are the same except for the value (and the units) of the objective function 

(2.8657x106 ft3 compared to $57,313). Note that this would not be the case if the cost of 

withdrawing water differed among the seven wells.

Mixed-Binary Linear Formulation

In this second formulation, the costs of both installing and operating the wells will be 

considered. This change requires that seven binary variables be defined for the problem 

(named BV1, BV2, and so forth), one of which is associated with each flow-rate decision 

variable. The operational costs of the system are $0.02/ft3; installation costs are $2,000 per 

well. The objective function of the mixed-binary linear problem is

Minimize

(76)

To ensure reliability of the system, a minimum of at least three wells must be installed, 

and the minimum withdrawal rate that is allowed at any well is 100 ft3/d (maximum pumping rates are still 20,000 ft3/d). A linear-summation constraint is added to the formulation 

that states that a minimum of at least three wells must be installed (see eq. 26):

. (77)

This constraint is specified in the SUMCON file. No changes to the HEDCON or SOLN 

files used in the linear formulation are required for this formulation.

The optimal solution was determined to be $63,598 and only three wells were selected 

for pumping: well Q1 pumps at a rate of 1,242 ft3/d, well Q4 at 694.1 ft3/d, and well Q7 at 

943.3 ft3/d. Note that well Q2, which pumped at a rate of only 78.2 ft3/d in the linear formulations, was not selected for pumping in this formulation because of the relatively high 

installation cost for the well in comparison to the amount of water produced. The installation cost for the three wells is $6,000 and the operational cost for the 1,000-day period is 

$57,598, which is close to the value of $57,313 determined for the linear formulation.

$0.02

ft3 ------------ ⎝ ⎠ ⎛ ⎞ Q1 $0.02

ft3 ------------ ⎝ ⎠ ⎛ ⎞Q2 $0.02ft3 ------------ ⎝ ⎠ ⎛ ⎞Q3 $0.02ft3 ------------ ⎝ ⎠ ⎛ ⎞Q4 $0.02ft3 ------------ ⎝ ⎠ ⎛ ⎞Q5 $0.02ft3 ------------ ⎝ ⎠ ⎛ ⎞Q6 

 $0.02

ft3 ------------- ⎝ ⎠ ⎛ ⎞ Q7 $2,000 ( )BV1 $2,000 ( )BV2 $2,000 ( )BV3 

 $2,000 ( )BV4 $2,000 ( )BV5 $2,000 ( )BV6 $2,000 ( )BV7 .

+++++

++ + +

++++

BV1 ++++++ BV2 BV3 BV4 BV5 BV6 BV7 ≥ 3
56 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Selected Input and Output Files for Sample Problem

Name File (dewater.nam):

GLOBAL 9 dewater.glo

LIST 10 dewater.lst

DIS 11 ..\data\dewater.dis

BAS6 12 ..\data\dewater.ba6

BCF6 13 ..\data\dewater.bc6

PCG 14 ..\data\dewater.pcg

GWM 15 ..\data\dewater.gwm

GWM File (dewater.gwm) for linear formulation:

#DEWATER Sample Problem, GWM file

#February 20, 2005

DECVAR ..\data\dewater.decvar

OBJFNC ..\data\dewater.objfnc

VARCON ..\data\dewater.varcon

HEDCON ..\data\dewater.hedcon

SOLN ..\data\dewater.soln

Decision variable (DECVAR) File (dewater.decvar) for linear formulation:

#DEWATER Sample Problem, DECVAR file

#February 20, 2005

 1 #1-IPRN

 7 0 0 #2-NFVAR NEVAR NBVAR

 Q1 1 1 7 14 W Y 1 #3a-FVNAME NC LAY ROW COL FTYPE ...

 Q2 1 1 7 16 W Y 1 

 Q3 1 1 8 15 W Y 1

 Q4 1 1 9 14 W Y 1

 Q5 1 1 9 16 W Y 1

 Q6 1 1 11 17 W Y 1

 Q7 1 1 13 16 W Y 1

Objective Function (OBJFNC) File (dewater.objfnc) for linear formulation:

#DEWATER Sample Problem, OBJFNC file

#February 20, 2005

 1 #1-IPRN

 MIN WSDV #2-OBJTYP FNTYP

 7 0 0 #3-NFVOBJ NEVOBJ NBVOBJ

 Q1 1.0 #4-FVNAME FVOBJC

 Q2 1.0 

 Q3 1.0 

 Q4 1.0 

 Q5 1.0 

 Q6 1.0 

 Q7 1.0 
Sample Problems 57

Decision-variable constraints (VARCON) File (dewater.varcon) for linear 

formulation:

#DEWATER Sample Problem, VARCON file

#February 20, 2005

 1 #1-IPRN

 Q1 0.0d2 2.0d4 0.0d2 #2-FVNAME FVMIN FVMAX FVREF

 Q2 0.0d2 2.0d4 0.0d2

 Q3 0.0d2 2.0d4 0.0d2

 Q4 0.0d2 2.0d4 0.0d2

 Q5 0.0d2 2.0d4 0.0d2

 Q6 0.0d2 2.0d4 0.0d2

 Q7 0.0d2 2.0d4 0.0d2

Head constraints (HEDCON) File (dewater.hedcon) for linear formulation:

#DEWATER Sample Problem, HEDCON file

#February 20, 2005

 1 #1-IPRN

 10 0 0 0 #2-NHB NDD NDF NGD

 b-01 1 6 13 le 50.0 1 #3-HBNAME LAYH ROWH COLH TYPH BND NSP

 b-02 1 6 15 le 50.0 1

 b-03 1 6 17 le 50.0 1

 b-04 1 8 13 le 50.0 1

 b-05 1 8 17 le 50.0 1

 b-06 1 10 13 le 50.0 1

 b-07 1 10 15 le 50.0 1

 b-08 1 11 16 le 50.0 1

 b-09 1 14 16 le 50.0 1

 b-10 1 14 17 le 50.0 1

Solution and output control file (SOLN) File (dewater.soln) for linear formulation:

#DEWATER Sample Problem, SOLN file

#February 20, 2005

 LP #1-SOLNTYP

2 #4a-IRM

 1000 2000 #4b-LPITMAX BBITMAX

 0.5 #4c-DELTA

 1 10 0.5 #4d-NSIGDIG NPGNMX PGFACT

 1 1 #4e-BBITPRT RANGE

 0 #6a-IBASE

GWM File (dewatermb.gwm) for mixed-binary linear formulation:

#DEWATER Sample Problem, GWM file

#February 20, 2005

DECVAR ..\data\dewatermb.decvar

OBJFNC ..\data\dewatermb.objfnc

VARCON ..\data\dewatermb.varcon

HEDCON ..\data\dewater.hedcon

SUMCON ..\data\dewatermb.sumcon

SOLN ..\data\dewatermb.soln
58 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Decision variable (DECVAR) File (dewatermb.decvar) for mixed-binary linear 

formulation:

#DEWATER Sample Problem, DECVAR file

#February 20, 2005

 1 #1-IPRN

 7 0 7 #2-NFVAR NEVAR NBVAR

 Q1 1 1 7 14 W Y 1 #3a-FVNAME NC LAY ROW COL FTYPE ...

 Q2 1 1 7 16 W Y 1 

 Q3 1 1 8 15 W Y 1

 Q4 1 1 9 14 W Y 1

 Q5 1 1 9 16 W Y 1

 Q6 1 1 11 17 W Y 1

 Q7 1 1 13 16 W Y 1

 BV1 1 Q1 #5-BVNAME NDV BVLIST

 BV2 1 Q2

 BV3 1 Q3

 BV4 1 Q4

 BV5 1 Q5

 BV6 1 Q6

 BV7 1 Q7

Objective Function (OBJFNC) File (dewatermb.objfnc) for mixed-binary linear 

formulation:

#DEWATER Sample Problem, OBJFNC file

#February 20, 2005

 1 #1-IPRN

 MIN WSDV #2-OBJTYP FNTYP

 7 0 7 #3-NFVOBJ NEVOBJ NBVOBJ

 Q1 0.02 #4-FVNAME FVOBJC

 Q2 0.02 

 Q3 0.02 

 Q4 0.02 

 Q5 0.02 

 Q6 0.02 

 Q7 0.02 

 BV1 2000.0 #6-BVNAME BVOBJC

 BV2 2000.0 

 BV3 2000.0 

 BV4 2000.0 

 BV5 2000.0 

 BV6 2000.0 

 BV7 2000.0 
Sample Problems 59

Decision-variable constraints (VARCON) File (dewatermb.varcon) for mixed-binary 

linear formulation:

#DEWATER Sample Problem, VARCON file

#February 20, 2005

 1 #1-IPRN

 Q1 1.0d2 2.0d4 0.0d2 #2-FVNAME FVMIN FVMAX FVREF

 Q2 1.0d2 2.0d4 0.0d2

 Q3 1.0d2 2.0d4 0.0d2

 Q4 1.0d2 2.0d4 0.0d2

 Q5 1.0d2 2.0d4 0.0d2

 Q6 1.0d2 2.0d4 0.0d2

 Q7 1.0d2 2.0d4 0.0d2

Linear-summation constraint file (SUMCON) File (dewatermb.sumcon) for mixedbinary linear formulation:

#DEWATER Sample Problem, SUMCON file

#February 20, 2005

 1 #1-IPRN

 1 #2-SMCNUM

 Binlower 7 ge 3.0 #3a-SMCNAME NTERMS TYPE RHS 

 BV1 1 #3b-GVNAME GVCOEFF

 BV2 1

 BV3 1

 BV4 1 

 BV5 1 

 BV6 1 

 BV7 1 

GLOBAL File (dewater.glo) for linear formulation with objective-function 

coefficients equal to 1.0:

MODFLOW-2000

U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUND-WATER FLOW MODEL

VERSION MF2K_GWM 1.0.0 022005,FROM MF2K V1.13.00

This model run produced both GLOBAL and LIST files. This is the GLOBAL file.

 GLOBAL LISTING FILE: dewater.glo

 UNIT 9

 OPENING dewater.lst

 FILE TYPE:LIST UNIT 10 STATUS:REPLACE

 FORMAT:FORMATTED ACCESS:SEQUENTIAL 

 OPENING ..\data\dewater.dis

 FILE TYPE:DIS UNIT 11 STATUS:OLD 

 FORMAT:FORMATTED ACCESS:SEQUENTIAL 

 OPENING ..\data\dewater.ba6

 FILE TYPE:BAS6 UNIT 12 STATUS:OLD 

 FORMAT:FORMATTED ACCESS:SEQUENTIAL 

 OPENING ..\data\dewater.bc6

 FILE TYPE:BCF6 UNIT 13 STATUS:OLD 

 FORMAT:FORMATTED ACCESS:SEQUENTIAL 
60 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

 OPENING ..\data\dewater.pcg

 FILE TYPE:PCG UNIT 14 STATUS:OLD 

 FORMAT:FORMATTED ACCESS:SEQUENTIAL 

 OPENING ..\data\dewater.gwm

 FILE TYPE:GWM UNIT 15 STATUS:OLD 

 FORMAT:FORMATTED ACCESS:SEQUENTIAL 

 THE FREE FORMAT OPTION HAS BEEN SELECTED

 DISCRETIZATION INPUT DATA READ FROM UNIT 11

 #DEWATER Sample Problem

 #February 20, 2005

 1 LAYERS 20 ROWS 30 COLUMNS

 1 STRESS PERIOD(S) IN SIMULATION

 MODEL TIME UNIT IS DAYS

 MODEL LENGTH UNIT IS UNDEFINED

 THE GROUND-WATER TRANSPORT PROCESS IS INACTIVE

 THE OBSERVATION PROCESS IS INACTIVE

 THE SENSITIVITY PROCESS IS INACTIVE

 THE PARAMETER-ESTIMATION PROCESS IS INACTIVE

 MODE: FORWARD

 Confining bed flag for each layer:

 0

 5450 ELEMENTS OF GX ARRAY USED OUT OF 5450

 600 ELEMENTS OF GZ ARRAY USED OUT OF 600

 600 ELEMENTS OF IG ARRAY USED OUT OF 600

 DELR = 100.000 

 DELC = 100.000 

 TOP ELEVATION OF LAYER 1 = 100.000 

 MODEL LAYER BOTTOM EL. = 100.000 FOR LAYER 1

STRESS PERIOD LENGTH TIME STEPS MULTIPLIER FOR DELT SS FLAG

-------------------------------------------------------------------------

 1 1000.000 1 1.000 SS

 STEADY-STATE SIMULATION

 PCG2 -- CONJUGATE GRADIENT SOLUTION PACKAGE, VERSION 2.4, 12/29/98

 MAXIMUM OF 50 CALLS OF SOLUTION ROUTINE

 MAXIMUM OF 5 INTERNAL ITERATIONS PER CALL TO SOLUTION ROUTINE

 MATRIX PRECONDITIONING TYPE : 1

 1700 ELEMENTS IN X ARRAY ARE USED BY PCG

 1750 ELEMENTS IN IX ARRAY ARE USED BY PCG

 2400 ELEMENTS IN Z ARRAY ARE USED BY PCG

 1700 ELEMENTS OF X ARRAY USED OUT OF 1700

 2400 ELEMENTS OF Z ARRAY USED OUT OF 2400

 1750 ELEMENTS OF IX ARRAY USED OUT OF 1750
Sample Problems 61

 SOLUTION BY THE CONJUGATE-GRADIENT METHOD

 -------------------------------------------

 MAXIMUM NUMBER OF CALLS TO PCG ROUTINE = 50

 MAXIMUM ITERATIONS PER CALL TO PCG = 5

 MATRIX PRECONDITIONING TYPE = 1

 RELAXATION FACTOR (ONLY USED WITH PRECOND. TYPE 1) = 0.10000E+01

 PARAMETER OF POLYNOMIAL PRECOND. = 2 (2) OR IS CALCULATED : 0

 HEAD CHANGE CRITERION FOR CLOSURE = 0.10000E-07

 RESIDUAL CHANGE CRITERION FOR CLOSURE = 0.10000E+01

 PCG HEAD AND RESIDUAL CHANGE PRINTOUT INTERVAL = 1

 PRINTING FROM SOLVER IS LIMITED(1) OR SUPPRESSED (>1) = 0

 DAMPING PARAMETER = 0.10000E+01

GWM1 -- GROUND-WATER MANAGEMENT PROCESS, VERSION 1.0.0 022005

INPUT READ FROM UNIT 15

----------------------------------------------------------------------

 Reading GWM Input

----------------------------------------------------------------------

 #DEWATER Sample Problem, GWM file

 #February 20, 2005

 OPENING DECISION-VARIABLE FILE ON UNIT 99:

..\data\dewater.decvar 

 #DEWATER Sample Problem, DECVAR file

 #February 20, 2005

 NO. OF FLOW-RATE DECISION VARIABLES (NFVAR) 7

 NO. OF EXTERNAL DECISION VARIABLES (NEVAR): 0

 BINARY VARIABLES ARE NOT ACTIVE.

 FLOW-RATE VARIABLES:

 FRACTION

 NUMBER NAME TYPE LAY ROW COL OF FLOW

 ----------------------------------------------------------

 1 Q1 WITHDRAWAL 1 7 14 1.0000

AVAILABLE IN STRESS PERIODS: 1 

 2 Q2 WITHDRAWAL 1 7 16 1.0000

AVAILABLE IN STRESS PERIODS: 1 

 3 Q3 WITHDRAWAL 1 8 15 1.0000

AVAILABLE IN STRESS PERIODS: 1 

 4 Q4 WITHDRAWAL 1 9 14 1.0000

AVAILABLE IN STRESS PERIODS: 1 

 5 Q5 WITHDRAWAL 1 9 16 1.0000

AVAILABLE IN STRESS PERIODS: 1 

 6 Q6 WITHDRAWAL 1 11 17 1.0000

AVAILABLE IN STRESS PERIODS: 1 

 7 Q7 WITHDRAWAL 1 13 16 1.0000

AVAILABLE IN STRESS PERIODS: 1 
62 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

 490 BYTES OF MEMORY ALLOCATED TO STORE DATA FOR DECISION VARIABLES

 CLOSING DECISION-VARIABLE FILE

 OPENING OBJECTIVE-FUNCTION FILE ON UNIT 99:

..\data\dewater.objfnc 

 #DEWATER Sample Problem, OBJFNC file

 #February 20, 2005

 OBJECTIVE TYPE: MIN FUNCTION TYPE: WSDV

 NO. OF FLOW-RATE DECISION VARIABLES IN OBJECTIVE FUNCTION (NFVOBJ): 7

 NO. OF EXTERNAL DECISION VARIABLES IN OBJECTIVE FUNCTION (NEVOBJ): 0

 NO. OF BINARY DECISION VARIABLES IN OBJECTIVE FUNCTION (NBVOBJ): 0

 OBJECTIVE FUNCTION:

 MIN + 1.00E+00 Q1 + 1.00E+00 Q2 + 1.00E+00 Q3

+ 1.00E+00 Q4 + 1.00E+00 Q5 + 1.00E+00 Q6

+ 1.00E+00 Q7

 28 BYTES OF MEMORY ALLOCATED TO STORE DATA FOR OBJECTIVE-FUNCTION

 CLOSING OBJECTIVE-FUNCTION FILE

 OPENING DECISION-VARIABLE CONSTRAINTS FILE

 ON UNIT 99:

..\data\dewater.varcon 

 #DEWATER Sample Problem, VARCON file

 #February 20, 2005

 FLOW RATE VARIABLES:

 MINIMUM MAXIMUM REFERENCE

 NUMBER NAME FLOW RATE FLOW RATE FLOW RATE

----------------------------------------------------------------

 1 Q1 0.000E+00 2.000E+04 0.000E+00

 2 Q2 0.000E+00 2.000E+04 0.000E+00

 3 Q3 0.000E+00 2.000E+04 0.000E+00

 4 Q4 0.000E+00 2.000E+04 0.000E+00

 5 Q5 0.000E+00 2.000E+04 0.000E+00

 6 Q6 0.000E+00 2.000E+04 0.000E+00

 7 Q7 0.000E+00 2.000E+04 0.000E+00

 CLOSING DECISION-VARIABLE CONSTRAINTS FILE

 OPENING HEAD CONSTRAINTS FILE

 ON UNIT 99:

..\data\dewater.hedcon 

 #DEWATER Sample Problem, HEDCON file

 #February 20, 2005
Sample Problems 63

 HEAD CONSTRAINTS:

 RIGHT-HAND STRESS

 NUMBER NAME LAY ROW COL TYPE SIDE PERIOD

 ---------------------------------------------------------------

 1 b-01 1 6 13 < 5.0000E+01 1

 2 b-02 1 6 15 < 5.0000E+01 1

 3 b-03 1 6 17 < 5.0000E+01 1

 4 b-04 1 8 13 < 5.0000E+01 1

 5 b-05 1 8 17 < 5.0000E+01 1

 6 b-06 1 10 13 < 5.0000E+01 1

 7 b-07 1 10 15 < 5.0000E+01 1

 8 b-08 1 11 16 < 5.0000E+01 1

 9 b-09 1 14 16 < 5.0000E+01 1

 10 b-10 1 14 17 < 5.0000E+01 1

 660 BYTES OF MEMORY ALLOCATED TO STORE DATA FOR HEAD CONSTRAINTS

 CLOSING HEAD CONSTRAINTS FILE

 OPENING SOLUTION FILE ON UNIT 99:

..\data\dewater.soln 

 #DEWATER Sample Problem, SOLN file

 #February 20, 2005

 SOLNTYP IS LP: GWM WILL COMPLETE A SINGLE ITERATION OF THE LINEAR PROBLEM.

 IRM EQUALS 2: RESPONSE MATRIX WILL BE CALCULATED BY GWM

 BUT NOT WRITTEN TO THE RESPONSE FILE.

 MAXIMUM NUMBER OF LP ITERATIONS: 1000

 MAXIMUM NUMBER OF BRANCH AND BOUND ITER: 2000

 PERTURBATION VALUE: 0.50D+00

 MAXIMUM NUMBER OF PERTURBATION ATTEMPTS: 10

 PERTURBATION ADJUSTMENT FACTOR (PGFACT): 0.50000

 OUTPUT FROM BRANCH-AND-BOUND ALGORITHM WILL NOT BE PRINTED.

 BASE PUMPING RATES TAKEN FROM FVREF SPECIFIED IN VARCON INPUT FILE

 PROBLEM SIZE

 NUMBER OF VARIABLES (INCLUDING SLACKS) 17

 NUMBER OF CONSTRAINT EQUATIONS 10

 6268 BYTES OF MEMORY ALLOCATED FOR RESPONSE MATRIX ALGORITHM

 CLOSING SOLUTION AND OUTPUT FILE
64 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

---------------------------------------------------------------

 Solution Algorithm

---------------------------------------------------------------

 Begin Solution Algorithm

 Running Base Flow Process Simulation

 Status of Simulation-Based Constraints 

 Constraint Type Name Status Distance To RHS

 --------------- ---- ------ ---------------

 Head Bound b-01 Not Met 2.1724E+01

 Head Bound b-02 Not Met 2.0345E+01

 Head Bound b-03 Not Met 1.8966E+01

 Head Bound b-04 Not Met 2.1724E+01

 Head Bound b-05 Not Met 1.8966E+01

 Head Bound b-06 Not Met 2.1724E+01

 Head Bound b-07 Not Met 2.0345E+01

 Head Bound b-08 Not Met 1.9655E+01

 Head Bound b-09 Not Met 1.9655E+01

 Head Bound b-10 Not Met 1.8966E+01

 

 Calculating Response Matrix

 Perturb Flow Variable 1

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 2

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 3

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 4

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 5

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 6

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 7

 By Perturbation Value: -1.000000E+04

 Average Number of Significant Digits in Matrix 1.011429E+01

 

 Solving Linear Program

 Optimal Solution Found
Sample Problems 65

----------------------------------------------------------------------

 Ground-Water Management Solution

----------------------------------------------------------------------

 OPTIMAL SOLUTION FOUND 

 OPTIMAL RATES FOR EACH FLOW VARIABLE 

 ---------------------------------------

Variable Withdrawal Injection Contribution

Name Rate Rate To Objective

---------- -------------- ------------ ------------

 Q1 1.077390E+03 1.077390E+06

 Q2 7.823877E+01 7.823877E+04

 Q3 0.000000E+00 0.000000E+00

 Q4 7.689506E+02 7.689506E+05

 Q5 0.000000E+00 0.000000E+00

 Q6 0.000000E+00 0.000000E+00

 Q7 9.410751E+02 9.410751E+05

 ------------ ------------ ------------

TOTALS 2.865655E+03 0.000000E+00 2.865655E+06

 OBJECTIVE FUNCTION VALUE 2.865655E+06

 BINDING CONSTRAINTS 

Constraint Type Name Status Shadow Price

--------------- ---- ------ ------------

Head Bound b-01 Binding -2.7273E+04

Head Bound b-03 Binding -3.2593E+04

Head Bound b-06 Binding -3.1185E+04

Head Bound b-10 Binding -5.1544E+04

Binding constraint and range analysis values are determined from the linear 

program and based on the response matrix approximation of the flow-process.

 RANGE ANALYSIS 

 Constraint Ranges

Lower/Upper Bound are the values of the RHS beyond which basis will change.

 Leaving is the variable which will leave the basis. 

 Entering is the variable which will enter the basis.

 If the entering or leaving variable is a constraint name,

 then the constraint slack variable is active
66 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Constraint Original Lower/Upper

Name Slack RHS Bound Entering Leaving

---------- ---------- ---------- ---------- ---------- -------

b-01 0.0000E+00 5.0000E+01 4.9477E+01 b-03 Q2

5.3228E+01 b-01 b-04

b-02 2.0745E+00 5.0000E+01 4.7926E+01 b-01 b-02

Infinity ----- No Change -----

b-03 0.0000E+00 5.0000E+01 4.3065E+01 b-01 Q1

5.0317E+01 b-03 Q2

b-04 2.0528E+00 5.0000E+01 4.7947E+01 b-01 b-04

Infinity ----- No Change -----

b-05 1.1167E+00 5.0000E+01 4.8883E+01 Q3 b-05

Infinity ----- No Change -----

b-06 0.0000E+00 5.0000E+01 4.7939E+01 b-03 Q2

5.2635E+01 b-06 Q4

b-07 2.6182E+00 5.0000E+01 4.7382E+01 Q3 b-07

Infinity ----- No Change -----

b-08 1.8584E+00 5.0000E+01 4.8142E+01 Q6 b-08

Infinity ----- No Change -----

b-09 1.0158E+00 5.0000E+01 4.8984E+01 b-10 b-09

Infinity ----- No Change -----

b-10 0.0000E+00 5.0000E+01 4.7205E+01 b-03 Q2

5.0850E+01 b-10 b-09

 Objective-Function Coefficient Ranges 

Lower/Upper Bound are the values of the coefficients beyond which basis 

will change.

 Leaving is the variable which will leave the basis. 

 Entering is the variable which will enter the basis.

 If the entering or leaving variable is a constraint name,

 then the constraint slack variable is active

 Basic variables are shown with zero reduced cost

Variable Reduced Original Lower/Upper

Name Cost Coefficient Bound Entering Leaving

---------- ---------- ---------- ---------- ---------- -------

Q1 0.0000E+00 1.0000E+03 9.1368E+02 b-01 Q2

1.0669E+03 Q3 Q2

Q2 0.0000E+00 1.0000E+03 8.6811E+02 b-03 Q1

1.0438E+03 Q3 Q2

Q3 1.5770E+01 1.0000E+03 9.8423E+02 Q3 Q2

Infinity ----- No Change -----

Q4 0.0000E+00 1.0000E+03 8.9312E+02 b-06 Q2

1.0471E+03 Q3 Q2
Sample Problems 67

Q5 4.4085E+01 1.0000E+03 9.5592E+02 Q5 Q2

Infinity ----- No Change -----

Q6 7.4018E+01 1.0000E+03 9.2598E+02 Q6 Q2

Infinity ----- No Change -----

Q7 0.0000E+00 1.0000E+03 6.7387E+02 b-10 Q2

1.1286E+03 Q6 Q2

----------------------------------------------------------------------

 Final Flow Process Simulation

----------------------------------------------------------------------

 Status of Simulation-Based Constraints 

 Using Optimal Flow Rate Variable Values

 Constraint Type Name Status Distance To RHS

 --------------- ---- ------ ---------------

 Head Bound b-01 Near-Binding 3.4739E-07

 Head Bound b-02 Satisfied 2.0745E+00

 Head Bound b-03 Near-Binding 2.8216E-07

 Head Bound b-04 Satisfied 2.0528E+00

 Head Bound b-05 Satisfied 1.1167E+00

 Head Bound b-06 Near-Binding 2.6732E-07

 Head Bound b-07 Satisfied 2.6182E+00

 Head Bound b-08 Satisfied 1.8584E+00

 Head Bound b-09 Satisfied 1.0158E+00

 Head Bound b-10 Near-Binding 1.8573E-07

 Because of precision limitations and possible nonlinear behavior, 

 the status of binding constraints computed directly by the flow process 

 may differ slightly from those computed using the linear program. 

Sample Problem 2: SEAWATER

In this sample problem, seawater intrusion is to be controlled in a coastal area where 

ground water is pumped for water supply. Because of the contrast in water density between 

freshwater and saltwater, a density-dependent ground-water flow model is needed to rigorously address the problem of seawater intrusion caused by ground-water pumping; specifically, to adequately protect the coastal aquifers from the intrusion of salty water. As a first 

approximation, however, a constant-density, MODFLOW ground-water flow model has 

been developed that uses hydraulic controls at the coast as a surrogate for density effects. 

Such an approach has been used by Reichard (1995), Nishikawa (1998), and Reichard and 

others (2003).

The ground-water flow system consists of an upper unconfined aquifer with a uniform 

bottom elevation of 10 ft below local sea level and a lower confined aquifer with a bottom 

elevation of 20 ft below local sea level. The flow system extends 6,000 ft landward from 

the coast to a mountainous area underlain by impermeable rocks. The flow system is 

simulated by two model layers separated by a confining unit with a very low hydraulic 

conductivity. The width of the modeled area parallel to the coast is 4,000 ft. The model grid 

consists of 20 rows and 30 columns, and each grid cell is 200 ft by 200 ft (fig. 4). The 

hydraulic conductivity of each aquifer is homogeneous and isotropic. The hydraulic conductivity of the unconfined aquifer is 5 ft/d; the transmissivity of the confined aquifer 

is 800 ft2/d. The confining unit is modeled implicitly by use of a vertical conductance 

between the layers equal to 0.05 d-1. An existing production well is placed in the lower 

aquifer at model cell (2, 13, 20) and withdraws at a rate of 5,000 ft3/d.
68 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

1

20

10

2

3

4

5

6

7

8

9

11

12

13

14

15

16

17

18

19

PRODUCTION WELL

HORIZONTAL-GRADIENT 

CONSTRAINT LOCATION

CANDIDATE WITHDRAWAL WELL

CANDIDATE INJECTION WELL

EXPLANATION

COLUMNS

ROWS

No-flow Boundary

Constant head of 0 feet

VERTICAL HEAD-DIFFERENCE CONSTRAINT LOCATION

27 28 29 30

COLUMNS

Candidate

injection

well

LAYERS

1

2

Seawater

W5

I1

I2

W4

W2 W3

W1

I1

W1

I3

Figure 4. Model grid for SEAWATER sample problem.

Seawater at the coast is modeled by constant heads of 0 ft in column 30 of the top layer 

of the model and in columns 29 and 30 of the bottom layer of the model (fig. 4 inset). The 

model uses no-flow boundary conditions at the contact between the bottom of the lower 

aquifer and the impermeable rocks and along the western, northern, and southern boundaries of each layer. The upper aquifer is recharged at a uniform rate of 0.002 ft/d. As an 

example of the use of the Sensitivity Process to input model parameters, recharge is 

specified in the model by use of the RECH and SEN input files.
Sample Problems 69

In the absence of pumping at the existing production well, there is a uniform gradient 

toward the coast of 2.67 ft/400 ft (0.00668 ft/ft) from column 28 to column 30 in the upper 

aquifer and of 2.62 ft/400 ft (0.00655 ft/ft) from column 28 to column 30 in the lower aquifer. There is also a downward head difference everywhere along the coast of 0.11 ft from 

the unconfined aquifer to the lower confined aquifer in column 29. With pumping at the 

production well, however, the gradients are reduced to a minimum value of 2.33 ft/400 ft 

(0.00583 ft/ft) from column 28 to column 30 in the upper aquifer and 2.28 ft/400 ft (0.0057 

ft/ft) from column 28 to column 30 in the lower aquifer. Also, the downward head difference everywhere along the coast is lowered to about 0.10 ft.

The MODFLOW input files consist of a NAME file, a DIS-Package file, a BAS6-

Package file, a BCF6-Package file, a WEL-Package file, a RCH-Package file, a SENProcess file, and a PCG-Package file. Because the model simulates steady-state flow, 

the simulation time of the management problem is arbitrarily set to 1.0 day (NPER = 1, 

ITMUNI=4, and PERLEN = 1.0) in the discretization file. A very small head-change 

convergence criterion of 1.0E-8 ft is used.

Nonlinear Formulation

Water managers would like to develop additional ground-water supplies from 

the two aquifers at the five candidate withdrawal locations shown in figure 4 (W1, W2, 

and so forth), and simultaneously maintain a seaward hydraulic gradient of 1.50 ft/400 ft 

(0.00375 ft/ft) along the coast and a minimum downward gradient from the upper to the 

lower aquifer of 0.05 ft in column 29 of the modeled area. Three candidate injection wells 

(I1, I2, and I3) near the coastline have been proposed to help control the resulting seawater 

intrusion. Withdrawals from wells W1 through W4 would be from the unconfined aquifer 

(layer 1), whereas withdrawals from well W5 would be distributed evenly between the 

upper and lower aquifers (layers 1 and 2). Injection at each of the three candidate injection 

wells also would be distributed evenly between the upper and lower aquifers. The maximum rate of withdrawal or injection at any of the wells is 10,000 ft3/d and the minimum 

rate of withdrawal or injection is 0.

The objective of the ground-water management problem is to maximize the sum of 

withdrawals from the five candidate withdrawal wells minus the amount of water that is 

reinjected to control seawater intrusion at the three candidate injection wells:

Maximize . (78)

GWM multiplies each term in the objective function by the length of the single, steady-state 

stress period (1 day).

A linear-summation constraint is added to ensure (1) that the sum of withdrawals is 

greater than or equal to the sum of injections and (2) that the value of the objective function 

is positive:

. (79)

Three types of head constraints are used to control hydraulic conditions along the 

coast. First, water-level mounding at the three candidate injection wells must be less than 

or equal to land surface at the three locations (5 ft above local sea level). Second, 10 vertical 

head-difference constraints are used at locations shown on figure 4 to maintain a downward 

hydraulic gradient from the upper aquifer to the lower aquifer along column 29 of the 

model. Finally, 20 horizontal head-gradient constraints are used to maintain a seaward 

W1 ++++ W2 W3 W4 W5 – I1 – I2 – I3

W1 ++++ W2 W3 W4 W5 ≥ I1 + + I2 I3
70 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

hydraulic gradient of at least 1.50 ft/400 ft (0.00375 ft/ft) from column 28 to column 30 in 

each aquifer; these constraints are located in the odd-numbered rows (1, 3, 5, and so forth) 

and are shown schematically by the inset illustration on figure 4.

Although the objective function and many of the constraints of the problem are linear, 

the head constraints in the water-table aquifer are nonlinear. Therefore, sequential linear 

programming is selected in the SOLN file to solve the nonlinear formulation (SOLNTYP 

is set to SLP). An initial perturbation value of -5 percent (DINIT equal to -0.05) of the maximum withdrawal or injection rate of each candidate well is specified, which results in an 

initial perturbation withdrawal rate of 500 ft3/d at each withdrawal well and -500 ft3/d at 

each injection well. A minimum perturbation value of -0.005 percent (DMIN equal to 

-5.0 x 10-5) is selected, which results in an asymptotic perturbation withdrawal rate of 

0.5 ft3/d at each withdrawal well and -0.5 ft3/d at each injection well. Convergence criteria 

on the withdrawal/injection rates, , of 1.0 ft3/d (SLPVCRIT equal to 1.0) and on the 

value of the objective function, , of 0.001 ft3 (SLPZCRIT equal to 0.001) are specified. 

Base withdrawal and injection rates are set to 0 by specifying IBASE = 0 in the SOLN file 

and FVREF = 0.0d0 in the VARCON file for each withdrawal and injection well.

The GWM input files necessary for the problem formulation are DECVAR, OBJFNC, 

VARCON, SUMCON, HEDCON, and SOLN. These files are listed at the end of the sample problem.

Three iterations of the sequential linear-programming algorithm were required to satisfy the two convergence criteria and (see GLOBAL file output at end of sample 

problem). The value of the objective function at the optimal solution (that is, net withdrawal) is 14,088 ft3 of water withdrawn each day. Nonzero withdrawals were calculated 

for three of the candidate withdrawal wells: W1 (10,000 ft3/d), W2 (3,301 ft3/d), and W5 

(1,678 ft3/d); therefore, total withdrawals are 14,978 ft3/d. All three injection wells were 

selected, with a total daily injection rate into the three wells of 890 ft3/d. Five of the hydraulic-gradient constraints were binding, and all of them were located in the bottom layer 

of the model (rows 5, 7, 13, 17, and 19). The “Solution Algorithm,” “Ground-Water 

Management Solution,” and “Final Flow Process Simulation” sections of the GLOBAL 

output file for this GWM run are listed at the end of the sample problem.

The optimal solution is highly dependent on the value of the hydraulic-gradient 

constraints specified at the coast; as the required gradient constraint is increased, the 

amount of net withdrawal decreases. For example, at a specified gradient of 1.7 ft/400 ft 

(0.00425 ft/ft), net withdrawal decreases to 10,666 ft3/d, and for a specified gradient of 

2.0 ft/400 ft (0.005 ft/ft), net withdrawal decreases to 544 ft3/d.

ε1

ε2

ε1 ε2
Sample Problems 71

Selected Input and Output Files for Sample Problem

Name File (seawater.nam):

GLOBAL 9 seawater.glo

LIST 10 seawater.lst

DIS 11 ..\data\seawater.dis

BAS6 12 ..\data\seawater.ba6

BCF6 13 ..\data\seawater.bc6

WEL 14 ..\data\seawater.wel

RCH 15 ..\data\seawater.rch

PCG 16 ..\data\seawater.pcg

SEN 17 ..\data\seawater.sen

GWM 18 ..\data\seawater.gwm

GWM File (seawater.gwm):

#SEAWATER Sample Problem, GWM file

#February 20, 2005

DECVAR ..\data\seawater.decvar

OBJFNC ..\data\seawater.objfnc

VARCON ..\data\seawater.varcon

SUMCON ..\data\seawater.sumcon

HEDCON ..\data\seawater.hedcon

SOLN ..\data\seawater.soln

Decision variable (DECVAR) File (seawater.decvar):

#SEAWATER Sample Problem, DECVAR file

#February 20, 2005

 1 #1-IPRN

 8 0 0 #2-NFVAR NEVAR NBVAR

 W1 1 1 7 12 W Y 1 #3a-FVNAME NC LAY ROW COL FTYPE ...

 W2 1 1 10 13 W Y 1

 W3 1 1 10 16 W Y 1

 W4 1 1 12 15 W Y 1

 W5 2 0 0 0 W Y 1

 0.500 1 15 11 #3b-RATIO LAY ROW COL

 0.500 2 15 11 

 I1 2 0 0 0 I Y 1

 0.500 1 7 27

 0.500 2 7 27

 I2 2 0 0 0 I Y 1

 0.500 1 11 27

 0.500 2 11 27

 I3 2 0 0 0 I Y 1

 0.500 1 15 27

 0.500 2 15 27
72 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Objective Function (OBJFNC) File (seawater.objfnc):

#SEAWATER Sample Problem, OBJFNC file

#February 20, 2005

 1 #1-IPRN

 MAX WSDV #2-OBJTYP FNTYP

 8 0 0 #3-NFVOBJ NEVOBJ NBVOBJ

 W1 1.00 #4-FVNAME FVOBJC

 W2 1.00

 W3 1.00

 W4 1.00

 W5 1.00

 I1 -1.00

 I2 -1.00

 I3 -1.00

Decision-variable constraints (VARCON) File (seawater.varcon):

#SEAWATER Sample Problem, VARCON file

#February 20, 2005

 1 #1-IPRN

 W1 0.0d0 1.0d4 0.0d0 #2-FVNAME FVMIN FVMAX FVREF

 W2 0.0d0 1.0d4 0.0d0

 W3 0.0d0 1.0d4 0.0d0

 W4 0.0d0 1.0d4 0.0d0

 W5 0.0d0 1.0d4 0.0d0

 I1 0.0d0 1.0d4 0.0d0

 I2 0.0d0 1.0d4 0.0d0

 I3 0.0d0 1.0d4 0.0d0

Linear-summation constraint (SUMCON) File (seawater.sumcon):

#SEAWATER Sample Problem, SUMCON file

#February 20, 2005

 1 #1-IPRN

 1 #2-SMCNUM

 rech-a 8 ge 0. #3a-SMCNAME NTERMS TYPE RHS 

 W1 1.0 #3b-GVNAME GVCOEFF

 W2 1.0

 W3 1.0

 W4 1.0

 W5 1.0

 I1 -1.0

 I2 -1.0

 I3 -1.0
Sample Problems 73

Head constraints (HEDCON) File (seawater.hedcon):

#SEAWATER Sample Problem, HEDCON file

#February 20, 2005

 1 #1-IPRN

 3 0 10 20 #2-NHB NDD NDF NGD

 m-01 1 7 27 le 5. 1 #3-HBNAME LAYH ROWH COLH TYPH ...

 m-02 1 11 27 le 5. 1

 m-03 1 15 27 le 5. 1

 hd-01 1 1 29 2 1 29 0.05 1 #5-HDIFNAME LAY1 ROW1 COL1 ...

 hd-02 1 3 29 2 3 29 0.05 1

 hd-03 1 5 29 2 5 29 0.05 1

 hd-04 1 7 29 2 7 29 0.05 1

 hd-05 1 9 29 2 9 29 0.05 1

 hd-06 1 11 29 2 11 29 0.05 1

 hd-07 1 13 29 2 13 29 0.05 1

 hd-08 1 15 29 2 15 29 0.05 1

 hd-09 1 17 29 2 17 29 0.05 1

 hd-10 1 19 29 2 19 29 0.05 1

 gd-01 1 1 28 1 1 30 400. 0.00375 1 #6-GRADNAME (L,R,C)1 ...

 gd-02 2 1 28 2 1 30 400. 0.00375 1

 gd-03 1 3 28 1 3 30 400. 0.00375 1

 gd-04 2 3 28 2 3 30 400. 0.00375 1

 gd-05 1 5 28 1 5 30 400. 0.00375 1

 gd-06 2 5 28 2 5 30 400. 0.00375 1

 gd-07 1 7 28 1 7 30 400. 0.00375 1

 gd-08 2 7 28 2 7 30 400. 0.00375 1

 gd-09 1 9 28 1 9 30 400. 0.00375 1

 gd-10 2 9 28 2 9 30 400. 0.00375 1

 gd-11 1 11 28 1 11 30 400. 0.00375 1

 gd-12 2 11 28 2 11 30 400. 0.00375 1

 gd-13 1 13 28 1 13 30 400. 0.00375 1

 gd-14 2 13 28 2 13 30 400. 0.00375 1

 gd-15 1 15 28 1 15 30 400. 0.00375 1

 gd-16 2 15 28 2 15 30 400. 0.00375 1

 gd-17 1 17 28 1 17 30 400. 0.00375 1

 gd-18 2 17 28 2 17 30 400. 0.00375 1

 gd-19 1 19 28 1 19 30 400. 0.00375 1

 gd-20 2 19 28 2 19 30 400. 0.00375 1

Solution and output control file (SOLN) File (seawater.soln):

#SEAWATER Sample Problem, SOLN file

#February 20, 2005

 SLP #1-SOLNTYP

 10 1000 0 #5a-SLPITMAX LPITMAX BBITMAX

1.0 0.001 -0.05 -0.00005 5 #5b-SLPVCRIT SLPZCRIT DINIT DMIN DSC

1 10 0.5 0.5 5 #5c-NSIGDIG NPGNMX PGFACT AFACT NINFMX

1 0 1 #5d-SLPITPRT BBITPRT RANGE

 0 #6a-IBASE
74 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Selected output from the GLOBAL File (seawater.glo):

---------------------------------------------------------------

 Solution Algorithm

---------------------------------------------------------------

 Begin Solution Algorithm

 Running Base Flow Process Simulation

 0 Well parameters

 1 Recharge parameters

 PARAMETER NAME:REC-SP1 TYPE:RCH CLUSTERS: 1

 Parameter value from package file is: 0.0000 

 This value has been changed to: 2.00000E-03, as read from

 the Sensitivity Process file

 MULTIPLIER ARRAY: NONE ZONE ARRAY: ALL

 1 PARAMETER HAS BEEN DEFINED IN ALL PACKAGES.

 (SPACE IS ALLOCATED FOR 500 PARAMETERS.)

 Status of Simulation-Based Constraints 

 Constraint Type Name Status Distance To RHS

 --------------- ---- ------ ---------------

 Head Bound m-01 Satisfied 3.7103E-01

 Head Bound m-02 Satisfied 4.8707E-01

 Head Bound m-03 Satisfied 5.1690E-01

 Head Difference hd-01 Satisfied 5.3533E-02

 Head Difference hd-02 Satisfied 5.3335E-02

 Head Difference hd-03 Satisfied 5.2875E-02

 Head Difference hd-04 Satisfied 5.2172E-02

 Head Difference hd-05 Satisfied 5.1312E-02

 Head Difference hd-06 Satisfied 5.0499E-02

 Head Difference hd-07 Satisfied 5.0027E-02

 Head Difference hd-08 Satisfied 5.0048E-02

 Head Difference hd-09 Satisfied 5.0374E-02

 Head Difference hd-10 Satisfied 5.0681E-02

 Head Gradient gd-01 Satisfied 9.4952E-01

 Head Gradient gd-02 Satisfied 8.9861E-01

 Head Gradient gd-03 Satisfied 9.4284E-01

 Head Gradient gd-04 Satisfied 8.9200E-01

 Head Gradient gd-05 Satisfied 9.2728E-01

 Head Gradient gd-06 Satisfied 8.7661E-01

 Head Gradient gd-07 Satisfied 9.0348E-01

 Head Gradient gd-08 Satisfied 8.5307E-01

 Head Gradient gd-09 Satisfied 8.7427E-01

 Head Gradient gd-10 Satisfied 8.2417E-01

 Head Gradient gd-11 Satisfied 8.4655E-01

 Head Gradient gd-12 Satisfied 7.9675E-01

 Head Gradient gd-13 Satisfied 8.3044E-01

 Head Gradient gd-14 Satisfied 7.8080E-01

 Head Gradient gd-15 Satisfied 8.3121E-01

 Head Gradient gd-16 Satisfied 7.8156E-01

 Head Gradient gd-17 Satisfied 8.4240E-01

 Head Gradient gd-18 Satisfied 7.9262E-01

 Head Gradient gd-19 Satisfied 8.5287E-01

 Head Gradient gd-20 Satisfied 8.0298E-01

 
Sample Problems 75

 Calculating Response Matrix

 Perturb Flow Variable 1

 By Perturbation Value: 5.000000E+02

 Perturb Flow Variable 2

 By Perturbation Value: 5.000000E+02

 Perturb Flow Variable 3

 By Perturbation Value: 5.000000E+02

 Perturb Flow Variable 4

 By Perturbation Value: 5.000000E+02

 Perturb Flow Variable 5

 By Perturbation Value: 5.000000E+02

 Perturb Flow Variable 6

 By Perturbation Value: -5.000000E+02

 Perturb Flow Variable 7

 By Perturbation Value: -5.000000E+02

 Perturb Flow Variable 8

 By Perturbation Value: -5.000000E+02

 Average Number of Significant Digits in Matrix 6.393939E+00

 

 Solving Linear Program

 Optimal Solution Found

 Objective Value -1.412443E+04

 Maximum Relative Change in Flow Variable 9.999000E-01

 SLP Algorithm: End Iteration 1

 SLP Algorithm: Begin Iteration 2

 Running Base Flow Process Simulation

 Status of Simulation-Based Constraints 

 Constraint Type Name Status Distance To RHS

 --------------- ---- ------ ---------------

 Head Bound m-01 Satisfied 2.0549E+00

 Head Bound m-02 Satisfied 1.8048E+00

 Head Bound m-03 Satisfied 1.9387E+00

 Head Difference hd-01 Satisfied 2.7829E-02

 Head Difference hd-02 Satisfied 2.7734E-02

 Head Difference hd-03 Satisfied 2.7586E-02

 Head Difference hd-04 Satisfied 2.7594E-02

 Head Difference hd-05 Satisfied 2.7646E-02

 Head Difference hd-06 Satisfied 2.9108E-02

 Head Difference hd-07 Satisfied 2.7611E-02

 Head Difference hd-08 Satisfied 2.8253E-02

 Head Difference hd-09 Satisfied 2.7601E-02

 Head Difference hd-10 Satisfied 2.7596E-02

 Head Gradient gd-01 Satisfied 4.9407E-02

 Head Gradient gd-02 Satisfied 6.4056E-03

 Head Gradient gd-03 Satisfied 4.5939E-02

 Head Gradient gd-04 Satisfied 2.9631E-03

 Head Gradient gd-05 Satisfied 4.0503E-02

 Head Gradient gd-06 Not Met 2.4294E-03

 Head Gradient gd-07 Satisfied 4.0970E-02

 Head Gradient gd-08 Not Met 2.3047E-03

 Head Gradient gd-09 Satisfied 4.2385E-02

 Head Gradient gd-10 Not Met 5.2509E-04

 Head Gradient gd-11 Satisfied 9.7447E-02

 Head Gradient gd-12 Satisfied 5.1164E-02

 Head Gradient gd-13 Satisfied 4.0931E-02

 Head Gradient gd-14 Not Met 1.9472E-03

 Head Gradient gd-15 Satisfied 6.5654E-02

 Head Gradient gd-16 Satisfied 2.0929E-02

 Head Gradient gd-17 Satisfied 4.0929E-02

 Head Gradient gd-18 Not Met 1.9924E-03
76 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

 Head Gradient gd-19 Satisfied 4.0902E-02

 Head Gradient gd-20 Not Met 2.0380E-03

 

 Calculating Response Matrix

 Perturb Flow Variable 1

 By Perturbation Value: 1.004000E+02

 Perturb Flow Variable 2

 By Perturbation Value: 1.004000E+02

 Perturb Flow Variable 3

 By Perturbation Value: 1.004000E+02

 Perturb Flow Variable 4

 By Perturbation Value: 1.004000E+02

 Perturb Flow Variable 5

 By Perturbation Value: 1.004000E+02

 Perturb Flow Variable 6

 By Perturbation Value: -1.004000E+02

 Perturb Flow Variable 7

 By Perturbation Value: -1.004000E+02

 Perturb Flow Variable 8

 By Perturbation Value: -1.004000E+02

 Average Number of Significant Digits in Matrix 5.742424E+00

 

 Solving Linear Program

 Optimal Solution Found

 Objective Value -1.408782E+04

 Maximum Relative Change in Flow Variable 7.213245E-03

 SLP Algorithm: End Iteration 2

 SLP Algorithm: Begin Iteration 3

 Running Base Flow Process Simulation

 Status of Simulation-Based Constraints 

 Constraint Type Name Status Distance To RHS

 --------------- ---- ------ ---------------

 Head Bound m-01 Satisfied 2.0505E+00

 Head Bound m-02 Satisfied 1.8023E+00

 Head Bound m-03 Satisfied 1.9355E+00

 Head Difference hd-01 Satisfied 2.7897E-02

 Head Difference hd-02 Satisfied 2.7802E-02

 Head Difference hd-03 Satisfied 2.7653E-02

 Head Difference hd-04 Satisfied 2.7658E-02

 Head Difference hd-05 Satisfied 2.7706E-02

 Head Difference hd-06 Satisfied 2.9156E-02

 Head Difference hd-07 Satisfied 2.7665E-02

 Head Difference hd-08 Satisfied 2.8304E-02

 Head Difference hd-09 Satisfied 2.7656E-02

 Head Difference hd-10 Satisfied 2.7652E-02

 Head Gradient gd-01 Satisfied 5.1911E-02

 Head Gradient gd-02 Satisfied 8.8922E-03

 Head Gradient gd-03 Satisfied 4.8431E-02

 Head Gradient gd-04 Satisfied 5.4380E-03

 Head Gradient gd-05 Satisfied 4.2950E-02

 Head Gradient gd-06 Near-Binding 1.4160E-08

 Head Gradient gd-07 Satisfied 4.3288E-02

 Head Gradient gd-08 Near-Binding 2.4126E-08

 Head Gradient gd-09 Satisfied 4.4562E-02

 Head Gradient gd-10 Satisfied 1.6376E-03

 Head Gradient gd-11 Satisfied 9.9181E-02

 Head Gradient gd-12 Satisfied 5.2901E-02

 Head Gradient gd-13 Satisfied 4.2891E-02

 Head Gradient gd-14 Near-Binding 2.8028E-08

 Head Gradient gd-15 Satisfied 6.7484E-02
Sample Problems 77

 Head Gradient gd-16 Satisfied 2.2754E-02

 Head Gradient gd-17 Satisfied 4.2935E-02

 Head Gradient gd-18 Near-Binding 1.4009E-08

 Head Gradient gd-19 Satisfied 4.2954E-02

 Head Gradient gd-20 Near-Binding 8.0249E-09

 

 Calculating Response Matrix

 Perturb Flow Variable 1

 By Perturbation Value: 2.048000E+01

 Perturb Flow Variable 2

 By Perturbation Value: 2.048000E+01

 Perturb Flow Variable 3

 By Perturbation Value: 2.048000E+01

 Perturb Flow Variable 4

 By Perturbation Value: 2.048000E+01

 Perturb Flow Variable 5

 By Perturbation Value: 2.048000E+01

 Perturb Flow Variable 6

 By Perturbation Value: -2.048000E+01

 Perturb Flow Variable 7

 By Perturbation Value: -2.048000E+01

 Perturb Flow Variable 8

 By Perturbation Value: -2.048000E+01

 Average Number of Significant Digits in Matrix 5.234848E+00

 

 Solving Linear Program

 Optimal Solution Found

 Objective Value -1.408782E+04

 Maximum Relative Change in Flow Variable 4.443988E-08

 SLP Algorithm: End Iteration 3

 Iterations have converged

----------------------------------------------------------------------

 Ground-Water Management Solution

----------------------------------------------------------------------

 OPTIMAL SOLUTION FOUND 

 OPTIMAL RATES FOR EACH FLOW VARIABLE 

 ---------------------------------------

Variable Withdrawal Injection Contribution

Name Rate Rate To Objective

---------- -------------- ------------ ------------

 W1 1.000000E+04 1.000000E+04

 W2 3.300918E+03 3.300918E+03

 W3 0.000000E+00 0.000000E+00

 W4 0.000000E+00 0.000000E+00

 W5 1.677461E+03 1.677461E+03

 I1 6.293964E+01 -6.293964E+01

 I2 5.321156E+02 -5.321156E+02

 I3 2.955083E+02 -2.955083E+02

 ------------ ------------ ------------

TOTALS 1.497838E+04 8.905636E+02 1.408782E+04
78 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

 OBJECTIVE FUNCTION VALUE 1.408782E+04

 BINDING CONSTRAINTS 

Constraint Type Name Status Shadow Price

--------------- ---- ------ ------------

Head Gradient gd-06 Binding -5.7711E+03

Head Gradient gd-08 Binding -9.5887E+02

Head Gradient gd-14 Binding -8.0729E+03

Head Gradient gd-18 Binding -1.4635E+02

Head Gradient gd-20 Binding -2.1473E+03

Maximum Flow Rate W1 Binding Not Available

Binding constraint and range analysis values are determined from the linear

program and based on the response matrix approximation of the flow-process.

 RANGE ANALYSIS 

 Constraint Ranges

Lower/Upper Bound are the values of the RHS beyond which basis will change.

 Leaving is the variable which will leave the basis. 

 Entering is the variable which will enter the basis.

 If the entering or leaving variable is a constraint name,

 then the constraint slack variable is active

Constraint Original Lower/Upper

Name Slack RHS Bound Entering Leaving

---------- ---------- ---------- ---------- ---------- --------

m-01 2.0505E+00 5.0000E+00 2.9495E+00 gd-06 m-01

Infinity ----- No Change -----

m-02 1.8023E+00 5.0000E+00 3.1977E+00 gd-18 m-02

Infinity ----- No Change -----

m-03 1.9355E+00 5.0000E+00 3.0645E+00 gd-20 m-03

Infinity ----- No Change -----

hd-01 2.7897E-02 5.0000E-02 -Infinity ----- No Change -----

7.7897E-02 gd-06 hd-01

hd-02 2.7802E-02 5.0000E-02 -Infinity ----- No Change -----

7.7802E-02 gd-06 hd-02

hd-03 2.7653E-02 5.0000E-02 -Infinity ----- No Change -----

7.7653E-02 gd-06 hd-03

hd-04 2.7658E-02 5.0000E-02 -Infinity ----- No Change -----

7.7658E-02 gd-08 hd-04

hd-05 2.7706E-02 5.0000E-02 -Infinity ----- No Change -----

7.7706E-02 gd-20 hd-05

hd-06 2.9156E-02 5.0000E-02 -Infinity ----- No Change -----

7.9156E-02 gd-20 hd-06

hd-07 2.7665E-02 5.0000E-02 -Infinity ----- No Change -----

7.7665E-02 gd-14 hd-07

hd-08 2.8304E-02 5.0000E-02 -Infinity ----- No Change -----

7.8304E-02 gd-18 hd-08
Sample Problems 79

hd-09 2.7656E-02 5.0000E-02 -Infinity ----- No Change -----

7.7656E-02 gd-18 hd-09

hd-10 2.7652E-02 5.0000E-02 -Infinity ----- No Change -----

7.7652E-02 gd-20 hd-10

gd-01 5.1911E-02 1.5000E+00 -Infinity ----- No Change -----

1.5519E+00 gd-06 gd-01

gd-02 8.8922E-03 1.5000E+00 -Infinity ----- No Change -----

1.5089E+00 gd-06 gd-02

gd-03 4.8431E-02 1.5000E+00 -Infinity ----- No Change -----

1.5484E+00 gd-06 gd-03

gd-04 5.4380E-03 1.5000E+00 -Infinity ----- No Change -----

1.5054E+00 gd-06 gd-04

gd-05 4.2950E-02 1.5000E+00 -Infinity ----- No Change -----

1.5429E+00 gd-06 gd-05

gd-06 0.0000E+00 1.5000E+00 1.4962E+00 gd-06 gd-04

1.5100E+00 gd-18 I1

gd-07 4.3288E-02 1.5000E+00 -Infinity ----- No Change -----

1.5433E+00 gd-08 gd-07

gd-08 0.0000E+00 1.5000E+00 1.4935E+00 gd-20 gd-10

1.5133E+00 gd-06 gd-04

gd-09 4.4562E-02 1.5000E+00 -Infinity ----- No Change -----

1.5446E+00 gd-20 gd-09

gd-10 1.6375E-03 1.5000E+00 -Infinity ----- No Change -----

1.5016E+00 gd-20 gd-10

gd-11 9.9181E-02 1.5000E+00 -Infinity ----- No Change -----

1.5992E+00 gd-20 gd-11

gd-12 5.2901E-02 1.5000E+00 -Infinity ----- No Change -----

1.5529E+00 gd-20 gd-12

gd-13 4.2891E-02 1.5000E+00 -Infinity ----- No Change -----

1.5429E+00 gd-14 gd-13

gd-14 0.0000E+00 1.5000E+00 1.4984E+00 gd-20 gd-10

1.5301E+00 gd-18 I1

gd-15 6.7484E-02 1.5000E+00 -Infinity ----- No Change -----

1.5675E+00 gd-18 gd-15

gd-16 2.2754E-02 1.5000E+00 -Infinity ----- No Change -----

1.5228E+00 gd-1 gd-16

gd-17 4.2935E-02 1.5000E+00 -Infinity ----- No Change -----

1.5429E+00 gd-18 gd-17

gd-18 0.0000E+00 1.5000E+00 1.4935E+00 gd-18 gd-16

1.5015E+00 gd-20 gd-10
80 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

gd-19 4.2954E-02 1.5000E+00 -Infinity ----- No Change -----

1.5430E+00 gd-20 gd-19

gd-20 0.0000E+00 1.5000E+00 1.4975E+00 gd-20 gd-10

1.5094E+00 gd-18 gd-16

rech-a 1.4088E+04 0.0000E+00 -Infinity ----- No Change -----

1.4088E+04 -- Degenerate Basis -

 Objective-Function Coefficient Ranges 

Lower/Upper Bound are the values of the coefficients beyond which basis 

will change.

 Leaving is the variable which will leave the basis. 

 Entering is the variable which will enter the basis.

 If the entering or leaving variable is a constraint name,

 then the constraint slack variable is active

 Basic variables are shown with zero reduced cost

Variable Reduced Original Lower/Upper

Name Cost Coefficient Bound Entering Leaving

--------- ---------- ---------- ---------- ---------- -------

W1 3.1892E-03 1.0000E+00 9.9681E-01 W1 W5 

Infinity ----- No Change -----

W2 0.0000E+00 1.0000E+00 9.9801E-01 gd-18 gd-10 

1.0020E+00 W1 W5 

W3 -3.4953E-03 1.0000E+00 Infinity ----- No Change -----

1.0035E+00 W3 W2 

W4 -3.5766E-03 1.0000E+00 Infinity ----- No Change -----

1.0036E+00 W4 W5 

W5 0.0000E+00 1.0000E+00 9.9466E-01 W1 W5 

1.0018E+00 gd-18 gd-10 

I1 0.0000E+00 -1.0000E+00 -1.9137E+00 gd-06 I1 

-9.4052E-01 gd-18 gd-10 

I2 0.0000E+00 -1.0000E+00 -1.0083E+00 gd-18 gd-10 

-7.8592E-01 W3 W2 

I3 0.0000E+00 -1.0000E+00 -1.1114E+00 gd-20 gd-16 

-9.9345E-01 gd-18 gd-10 
Sample Problems 81

----------------------------------------------------------------------

 Final Flow Process Simulation

----------------------------------------------------------------------

 Status of Simulation-Based Constraints 

 Using Optimal Flow Rate Variable Values

 Constraint Type Name Status Distance To RHS

 --------------- ---- ------ ---------------

 Head Bound m-01 Satisfied 2.0505E+00

 Head Bound m-02 Satisfied 1.8023E+00

 Head Bound m-03 Satisfied 1.9355E+00

 Head Difference hd-01 Satisfied 2.7897E-02

 Head Difference hd-02 Satisfied 2.7802E-02

 Head Difference hd-03 Satisfied 2.7653E-02

 Head Difference hd-04 Satisfied 2.7658E-02

 Head Difference hd-05 Satisfied 2.7706E-02

 Head Difference hd-06 Satisfied 2.9156E-02

 Head Difference hd-07 Satisfied 2.7665E-02

 Head Difference hd-08 Satisfied 2.8304E-02

 Head Difference hd-09 Satisfied 2.7656E-02

 Head Difference hd-10 Satisfied 2.7652E-02

 Head Gradient gd-01 Satisfied 5.1911E-02

 Head Gradient gd-02 Satisfied 8.8922E-03

 Head Gradient gd-03 Satisfied 4.8431E-02

 Head Gradient gd-04 Satisfied 5.4380E-03

 Head Gradient gd-05 Satisfied 4.2950E-02

 Head Gradient gd-06 Near-Binding 1.4530E-08

 Head Gradient gd-07 Satisfied 4.3288E-02

 Head Gradient gd-08 Near-Binding 1.3689E-08

 Head Gradient gd-09 Satisfied 4.4562E-02

 Head Gradient gd-10 Satisfied 1.6376E-03

 Head Gradient gd-11 Satisfied 9.9181E-02

 Head Gradient gd-12 Satisfied 5.2901E-02

 Head Gradient gd-13 Satisfied 4.2891E-02

 Head Gradient gd-14 Near-Binding 1.2701E-08

 Head Gradient gd-15 Satisfied 6.7484E-02

 Head Gradient gd-16 Satisfied 2.2754E-02

 Head Gradient gd-17 Satisfied 4.2935E-02

 Head Gradient gd-18 Near-Binding 1.2746E-08

 Head Gradient gd-19 Satisfied 4.2954E-02

 Head Gradient gd-20 Near-Binding 1.2570E-08

 Because of precision limitations and possible nonlinear behavior, 

 the status of binding constraints computed directly by the flow process 

 may differ slightly from those computed using the linear program. 
82 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Sample Problem 3: SUPPLY

This sample problem represents a transient water-supply problem in which total 

ground-water withdrawals over a 3-year period are limited by the amount of streamflow 

depletion allowed in two streams that are in hydraulic connection with the aquifer (fig. 5). 

The 3-year period is divided into 12 seasons (winter, spring, summer, and fall of each year), 

each of which is represented by a single stress period. The aquifer is confined and the area 

of interest is 6,000 ft long by 5,000 ft wide. The model consists of a single layer with 25 

rows and 30 columns; each model cell is 200 ft by 200 ft. The aquifer is homogenous and 

isotropic with a transmissivity of 5,000 ft2/d and a storage coefficient of 0.05. The modeled 

area is bounded on the east and west by no-flow conditions and on the north and south by 

constant heads that decrease in elevation from west to east (fig. 5). The aquifer is recharged 

at a rate of 0.0005 ft/d in the winter, 0.002 ft/d in the spring, 0 ft/d in the summer, and 

0.001 ft/d in the fall. The two streams are simulated with the MODFLOW Stream (STR1) 

Package (Prudic, 1989) by use of three stream segments. Both streams are 20 ft wide and 

have a streambed conductance of 20,000 ft2/d. The mainstem has a slope of 0.0025, 

whereas the tributary stream has a slope of 0.0010.

The MODFLOW input files consist of a NAME file, a DIS-Package file, a BAS6-

Package file, a BCF6-Package file, a RCH-Package file, a STR1-Package file, an 

OC-Package file, and a PCG-Package file.
Sample Problems 83

N

1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 21 20 22 23 24 25 26 27 28 29 30

1

20

10

2

3

4

5

6

7

8

9

11

12

13

14

15

16

17

18

19

1

3

4

1 CANDIDATE WITHDRAWAL WELL

RIVER

STREAMFLOW-DEPLETION CONSTRAINT LOCATION

2

21

22

23

24

25

EXPLANATION

COLUMNS

ROWS

No-Flow BoundaryNo-Flow Boundary

Constant heads that range from 100.5 to 86.0 feet

Constant heads that range from 80.5 to 66.0 feet

1

2

3

4

Figure 5. Model grid for SUPPLY sample problem.
84 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Nonlinear Formulation

The management objective is to maximize the value of ground water withdrawn from 

four wells over the 3-year period, while simultaneously limiting streamflow depletions at 

four streamflow-constraint locations along the two streams (fig. 5). Because it may be 

advantageous to have some of the wells pump at different rates during each season, multiple decision variables are defined at well sites 2 and 4. Well 2 is allowed to have one withdrawal rate during the winter (decision variable Q2a, stress periods 1, 5, and 9), another 

during spring (Q2b, stress periods 2, 6, and 10), another during summer (Q2c, stress periods 

3, 7, and 11), and another during fall (Q2d, stress periods 4, 8, and 12). Well 4 only withdraws during the spring and fall, and is allowed to have one withdrawal rate during the 

spring (Q4a, stress periods 2, 6, and 10) and a second withdrawal rate during the fall (Q4b, 

stress periods 4, 8, and 12). Well 1 must have a constant withdrawal rate for all 12 stress 

periods and well 3 is allowed to withdraw only during the second year (stress periods 5-8). 

The lower and upper bounds on withdrawals at each well are 0 and 50,000 ft3/d, respectively. There is a net benefit from pumping at each well of $1 per 1,000 ft3 of water 

withdrawn (that is, $0.001/ft3). 

Importation of water from an external source is allowed in each of the last four stress 

periods. These external variables are named Im9, Im10, Im11, and Im12. The maximum 

amount of water that can be imported each season is 1,000,000 ft3/d. The cost to import the 

water is $0.0012/ft3. Because this is a cost rather than a benefit, it appears in the objective 

function with a negative coefficient. 

The objective function of the problem is, therefore,

Maximize

(80)

and the units of the objective function are dollars (that is, dollars per cubic foot multiplied 

by cubic feet per day multiplied by days of pumping at each well).

There are constraints that specify the upper and lower bounds on the total amounts of 

water that must be withdrawn each season to meet water-supply demands. The upper bound 

on the total withdrawals needed each season is 80,000 ft3/d, whereas the lower bound on 

total withdrawals is 30,000 ft3/d during the first year, 25,000 ft3/d during the second year, 

and 45,000 ft3/d during the third year. These water-supply demands are specified using 

linear-summation constraints that include both the flow-rate and external variables.

Allowed streamflow depletions at constraint locations 1 and 3 (fig. 5) are 15,000 ft3/d 

during each season, those at constraint location 2 are 20,000 ft3/d during the summer and 

fall, and those at constraint location 4 are 30,000 ft3/d during the summer and fall.

Although the aquifer is confined and the objective function is linear, the streamflowdepletion constraints are nonlinear because of the presence of the head-dependent boundary conditions at the streams. Nonlinearities arise at these boundaries in two ways. First, 

there may be a nonlinear relation between pumping at the four well sites and seepage rates 

across the streambeds (see fig. 1); second, the iterative method that is used to calculate the 

$0.001

ft3 --------------- ⎝ ⎠ ⎛ ⎞Q1 $0.001

ft3 --------------- ⎝ ⎠ ⎛ ⎞ Q2a $0.001ft3 --------------- ⎝ ⎠ ⎛ ⎞Q2b $0.001ft3 --------------- ⎝ ⎠ ⎛ ⎞Q2c $0.001ft3 --------------- ⎝ ⎠ ⎛ ⎞ Q2d 

 $0.001

ft3 --------------- ⎝ ⎠ ⎛ ⎞Q3 $0.001

ft3 --------------- ⎝ ⎠ ⎛ ⎞ Q4a $0.001ft3 --------------- ⎝ ⎠ ⎛ ⎞Q4b $0.0012ft3 ------------------ ⎝ ⎠ ⎛ ⎞ – Im9 

 $0.0012

ft3 ------------------ ⎝ ⎠ ⎛ ⎞Im10 $0.0012

ft3 ------------------ ⎝ ⎠ ⎛ ⎞ – Im11 $0.0012ft3 ------------------ ⎝ ⎠ ⎛ ⎞ – – Im12 ,

++++

++ +
Sample Problems 85

stage of each stream reach can produce nonlinear relations between pumping rates and 

streamflow (Prudic, 1989, p. 10-11). Therefore, sequential linear programming is selected 

in the SOLN file to solve the nonlinear formulation (SOLNTYP is set to SLP). An initial 

perturbation value of 20 percent (DINIT equal to 0.2) of the maximum withdrawal rate of 

each candidate well is specified, which results in an initial perturbation withdrawal rate of 

10,000 ft3/d at each well. A minimum perturbation value of 0.002 percent (DMIN equal 

to 2.0 x 10-5) is selected, which results in an asymptotic perturbation withdrawal rate of 

1.0 ft3/d. Convergence criteria on the withdrawal/injection rates, , of 1.0 x 10-5 ft3/d 

(SLPVCRIT equal to 1.0 x 10-5) and on the value of the objective function, , of 

$1.0 x 10-4 (SLPZCRIT equal to 1.0 x 10-4) are specified. Base withdrawal rates are set to 

0 ft3/d by specifying IBASE = 0 in the SOLN file and FVREF = 0.0d0 in the VARCON file 

for each well. 

The GWM input files necessary for the problem formulation are: DECVAR, OBJFNC, 

VARCON, SUMCON, STRMCON, and SOLN. These files are listed at the end of the 

sample problem.

Three iterations of the sequential linear programming algorithm were required to satisfy the two convergence criteria and . The value of the objective function at the optimal solution is $53,022, which is the value of the water withdrawn over the 3-year period. 

All eight decision variables are active in the solution, although the rate calculated for well 

3 (decision variable Q3) is inconsequential. Withdrawal rates at the other three well sites 

range from 3,483 ft3/d at well 1 (decision variable Q1) to a maximum of 50,000 ft3/d at well 

2 during the winter (decision variable Q2a). Two of the four external variables are active in 

the solution. Six of the streamflow-depletion constraints are binding: at constraint site 1 in 

winter of the third year, at constraint site 2 in summer and fall of the third year, at constraint 

site 3 in spring of the third year, and at constraint site 4 in fall of the third and fourth years. 

The “Solution Algorithm,” “Ground-Water Management Solution,” and “Final Flow 

Process Simulation” sections of the GLOBAL output file for this GWM run are listed at the 

end of the sample problem.

Selected Input and Output Files for Sample Problem

Name File (supply.nam):

GLOBAL 9 supply.glo

LIST 10 supply.lst

DIS 11 ..\data\supply.dis

BAS6 12 ..\data\supply.ba6

BCF6 13 ..\data\supply.bc6

RCH 15 ..\data\supply.rch

STR 16 ..\data\supply.str

OC 17 ..\data\supply.oc

PCG 18 ..\data\supply.pcg

GWM 19 ..\data\supply.gwm

GWM File (supply.gwm):

#SUPPLY Sample Problem, GWM file

#February 20, 2005

DECVAR ..\data\supply.decvar

OBJFNC ..\data\supply.objfnc

VARCON ..\data\supply.varcon

SUMCON ..\data\supply.sumcon

STRMCON ..\data\supply.strmcon

SOLN ..\data\supply.soln

ε1

ε2

ε1 ε2
86 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Decision Variable (DECVAR) File (supply.decvar):

#SUPPLY Sample Problem, DECVAR file

#February 20, 2005

 1 #1-IPRN

 8 4 0 #2-NFVAR NEVAR NBVAR

 Q1 1 1 12 11 W Y 1-12 #3a-FVNAME NC LAY ROW COL ...

 Q2a 1 1 16 17 W Y 1:5:9

 Q2b 1 1 16 17 W Y 2:6:10

 Q2c 1 1 16 17 W Y 3:7:11

 Q2d 1 1 16 17 W Y 4:8:12

 Q3 1 1 11 22 W Y 5-8

 Q4a 1 1 14 25 W Y 2:6:10

 Q4b 1 1 14 25 W Y 4:8:12

 Im9 IM 9 #4-EVNAME ETYPE ESP 

 Im10 IM 10

 Im11 IM 11

 Im12 IM 12

Objective Function (OBJFNC) File (supply.objfnc):

#SUPPLY Sample Problem, OBJFNC file

#February 20, 2005

 1 #1-IPRN

 MAX WSDV #2-OBJTYP FNTYP

 8 4 0 #3-NFVOBJ NEVOBJ NBVOBJ

 Q1 0.001 #4-FVNAME FVOBJC

 Q2a 0.001 

 Q2b 0.001

 Q2c 0.001

 Q2d 0.001

 Q3 0.001

 Q4a 0.001

 Q4b 0.001

 Im9 -0.0012 #5-EVNAME EVOBJC

 Im10 -0.0012

 Im11 -0.0012

 Im12 -0.0012

Decision-variable constraints (VARCON) File (supply.varcon):

#SUPPLY Sample Problem, VARCON file

#February 20, 2005

 1 #1-IPRN

 Q1 0.0 5.0d4 0.0D2 #2-FVNAME FVMIN FVMAX FVREF

 Q2a 0.0 5.0d4 0.0D2

 Q2b 0.0 5.0d4 0.0D2

 Q2c 0.0 5.0d4 0.0D2

 Q2d 0.0 5.0d4 0.0D2

 Q3 0.0 5.0d4 0.0D2

 Q4a 0.0 5.0d4 0.0D2

 Q4b 0.0 5.0d4 0.0D2

 Im9 0.0 1.0d6 #3-EVNAME EVMIN EVMAX 

 Im10 0.0 1.0d6

 Im11 0.0 1.0d6

 Im12 0.0 1.0d6
Sample Problems 87

Part of the linear-summation constraints (SUMCON) File (supply.sumcon):

#SUPPLY Sample Problem, SUMCON file

#February 20, 2005

 1 #1-IPRN

 24 #2-SMCNUM

 p01u 2 le 80000. #3a-SMCNAME NTERMS TYPE RHS

 Q1 1.0 #3b-GVNAME GVCOEFF

 Q2a 1.0

 p01l 2 ge 30000.

 Q1 1.0

 Q2a 1.0

 p02u 3 le 80000.

 Q1 1.0

 Q2b 1.0

 Q4a 1.0

 p02l 3 ge 30000.

 Q1 1.0

 Q2b 1.0

 Q4a 1.0

**12 constraints deleted here**

 p09u 3 le 80000.

 Q1 1.0

 Q2a 1.0

 Im9 1.0

 p09l 3 ge 45000.

 Q1 1.0

 Q2a 1.0

 Im9 1.0

 p10u 4 le 80000.

 Q1 1.0

 Q2b 1.0

 Q4a 1.0

 Im10 1.0

 p10l 4 ge 45000.

 Q1 1.0

 Q2b 1.0

 Q4a 1.0

 Im10 1.0

 p11u 3 le 80000.

 Q1 1.0

 Q2c 1.0

 Im11 1.0

 p11l 3 ge 45000.

 Q1 1.0

 Q2c 1.0

 Im11 1.0

 p12u 4 le 80000.

 Q1 1.0

 Q2d 1.0

 Q4b 1.0

 Im12 1.0

 p12l 4 ge 45000.

 Q1 1.0

 Q2d 1.0

 Q4b 1.0

 Im12 1.0
88 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Streamflow constraints (STRMCON) File (supply.strmcon):

#SUPPLY Sample Problem, STRMCON file

#February 20, 2005

 1 #1-IPRN

 0 36 #2-NSF NSD

 S01.1 1 14 le 15000. 1 #4-SDNAME SEG REACH TYPSD BND NSP

 S01.2 1 14 le 15000. 2

 S01.3 1 14 le 15000. 3

 S01.4 1 14 le 15000. 4

 S01.5 1 14 le 15000. 5

 S01.6 1 14 le 15000. 6

 S01.7 1 14 le 15000. 7

 S01.8 1 14 le 15000. 8

 S01.9 1 14 le 15000. 9

 S01.10 1 14 le 15000. 10

 S01.11 1 14 le 15000. 11

 S01.12 1 14 le 15000. 12

 S02.3 1 21 le 20000. 3

 S02.4 1 21 le 20000. 4

 S02.7 1 21 le 20000. 7

 S02.8 1 21 le 20000. 8

 S02.11 1 21 le 20000. 11

 S02.12 1 21 le 20000. 12

 S03.1 2 8 le 15000. 1

 S03.2 2 8 le 15000. 2

 S03.3 2 8 le 15000. 3

 S03.4 2 8 le 15000. 4

 S03.5 2 8 le 15000. 5

 S03.6 2 8 le 15000. 6

 S03.7 2 8 le 15000. 7

 S03.8 2 8 le 15000. 8

 S03.9 2 8 le 15000. 9

 S03.10 2 8 le 15000. 10

 S03.11 2 8 le 15000. 11

 S03.12 2 8 le 15000. 12

 S04.3 3 5 le 30000. 3

 S04.4 3 5 le 30000. 4

 S04.7 3 5 le 30000. 7

 S04.8 3 5 le 30000. 8

 S04.11 3 5 le 30000. 11

 S04.12 3 5 le 30000. 12

Solution and output control file (SOLN) File (supply.soln):

#SUPPLY Sample Problem, SOLN file

#February 20, 2005

 SLP #1-SOLNTYP

 50 10000 2000 #5a-SLPITMAX LPITMAX BBITMAX

0.00001 0.0001 0.2 0.00002 2 #5b-SLPVCRIT SLPZCRIT DINIT DMIN DSC

1 10 0.5 0.5 5 #5c-NSIGDIG NPGNMX PGFACT AFACT NINFMX

 1 1 0 #5d-SLPITPRT BBITPRT RANGE

 0 #6a-IBASE
Sample Problems 89

Selected output from the GLOBAL File (supply.glo):

---------------------------------------------------------------

 Solution Algorithm

---------------------------------------------------------------

 Begin Solution Algorithm

 Running Flow Process Simulation

 for both Reference and Base 

 0 Recharge parameters

 0 Stream parameters

 Status of Simulation-Based Constraints 

 Constraint Type Name Status Distance To RHS

 --------------- ---- ------ ---------------

 Stream Depletion S01.1 Satisfied 1.5000E+04

 Stream Depletion S01.2 Satisfied 1.5000E+04

 Stream Depletion S01.3 Satisfied 1.5000E+04

 Stream Depletion S01.4 Satisfied 1.5000E+04

 Stream Depletion S01.5 Satisfied 1.5000E+04

 Stream Depletion S01.6 Satisfied 1.5000E+04

 Stream Depletion S01.7 Satisfied 1.5000E+04

 Stream Depletion S01.8 Satisfied 1.5000E+04

 Stream Depletion S01.9 Satisfied 1.5000E+04

 Stream Depletion S01.10 Satisfied 1.5000E+04

 Stream Depletion S01.11 Satisfied 1.5000E+04

 Stream Depletion S01.12 Satisfied 1.5000E+04

 Stream Depletion S02.3 Satisfied 2.0000E+04

 Stream Depletion S02.4 Satisfied 2.0000E+04

 Stream Depletion S02.7 Satisfied 2.0000E+04

 Stream Depletion S02.8 Satisfied 2.0000E+04

 Stream Depletion S02.11 Satisfied 2.0000E+04

 Stream Depletion S02.12 Satisfied 2.0000E+04

 Stream Depletion S03.1 Satisfied 1.5000E+04

 Stream Depletion S03.2 Satisfied 1.5000E+04

 Stream Depletion S03.3 Satisfied 1.5000E+04

 Stream Depletion S03.4 Satisfied 1.5000E+04

 Stream Depletion S03.5 Satisfied 1.5000E+04

 Stream Depletion S03.6 Satisfied 1.5000E+04

 Stream Depletion S03.7 Satisfied 1.5000E+04

 Stream Depletion S03.8 Satisfied 1.5000E+04

 Stream Depletion S03.9 Satisfied 1.5000E+04

 Stream Depletion S03.10 Satisfied 1.5000E+04

 Stream Depletion S03.11 Satisfied 1.5000E+04

 Stream Depletion S03.12 Satisfied 1.5000E+04

 Stream Depletion S04.3 Satisfied 3.0000E+04

 Stream Depletion S04.4 Satisfied 3.0000E+04

 Stream Depletion S04.7 Satisfied 3.0000E+04

 Stream Depletion S04.8 Satisfied 3.0000E+04

 Stream Depletion S04.11 Satisfied 3.0000E+04

 Stream Depletion S04.12 Satisfied 3.0000E+04

 

 Calculating Response Matrix

 Perturb Flow Variable 1

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 2

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 3

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 4
90 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 5

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 6

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 7

 By Perturbation Value: -1.000000E+04

 Perturb Flow Variable 8

 By Perturbation Value: -1.000000E+04

 Average Number of Significant Digits in Matrix 7.019841E+00

 

 Solving Linear Program

 Optimal Solution Found

 Objective Value -5.302102E+04

 Maximum Relative Change in Flow Variable 9.999800E-01

 Max Relative External Variable Change 9.999316E-01

 SLP Algorithm: End Iteration 1

 SLP Algorithm: Begin Iteration 2

 Running Base Flow Process Simulation

 Status of Simulation-Based Constraints 

 Constraint Type Name Status Distance To RHS

 --------------- ---- ------ ---------------

 Stream Depletion S01.1 Satisfied 4.4820E+02

 Stream Depletion S01.2 Satisfied 5.3103E+03

 Stream Depletion S01.3 Satisfied 5.8133E+03

 Stream Depletion S01.4 Satisfied 5.8671E+03

 Stream Depletion S01.5 Near-Binding 5.7007E-02

 Stream Depletion S01.6 Satisfied 5.2871E+03

 Stream Depletion S01.7 Satisfied 5.8121E+03

 Stream Depletion S01.8 Satisfied 5.8670E+03

 Stream Depletion S01.9 Near-Binding 6.2400E-02

 Stream Depletion S01.10 Satisfied 5.2871E+03

 Stream Depletion S01.11 Satisfied 5.8121E+03

 Stream Depletion S01.12 Satisfied 5.8670E+03

 Stream Depletion S02.3 Satisfied 3.2657E+00

 Stream Depletion S02.4 Near-Binding 2.5023E-01

 Stream Depletion S02.7 Satisfied 1.3145E+00

 Stream Depletion S02.8 Near-Binding 1.0742E-01

 Stream Depletion S02.11 Satisfied 1.3145E+00

 Stream Depletion S02.12 Near-Binding 1.0741E-01

 Stream Depletion S03.1 Satisfied 1.4120E+04

 Stream Depletion S03.2 Satisfied 2.5641E+00

 Stream Depletion S03.3 Satisfied 1.4209E+04

 Stream Depletion S03.4 Satisfied 1.2081E+04

 Stream Depletion S03.5 Satisfied 1.4032E+04

 Stream Depletion S03.6 Not Met 1.1239E+00

 Stream Depletion S03.7 Satisfied 1.4209E+04

 Stream Depletion S03.8 Satisfied 1.2081E+04

 Stream Depletion S03.9 Satisfied 1.4032E+04

 Stream Depletion S03.10 Not Met 1.1240E+00

 Stream Depletion S03.11 Satisfied 1.4209E+04

 Stream Depletion S03.12 Satisfied 1.2081E+04

 Stream Depletion S04.3 Satisfied 6.4478E+03

 Stream Depletion S04.4 Near-Binding 8.1515E-02

 Stream Depletion S04.7 Satisfied 6.4453E+03

 Stream Depletion S04.8 Near-Binding 1.2603E-01

 Stream Depletion S04.11 Satisfied 6.4453E+03

 Stream Depletion S04.12 Near-Binding 1.2603E-01

 
Sample Problems 91

 Calculating Response Matrix

 Perturb Flow Variable 1

 By Perturbation Value: -5.000500E+03

 Perturb Flow Variable 2

 By Perturbation Value: -5.000500E+03

 Perturb Flow Variable 3

 By Perturbation Value: -5.000500E+03

 Perturb Flow Variable 4

 By Perturbation Value: -5.000500E+03

 Perturb Flow Variable 5

 By Perturbation Value: -5.000500E+03

 Perturb Flow Variable 6

 By Perturbation Value: -5.000500E+03

 Perturb Flow Variable 7

 By Perturbation Value: -5.000500E+03

 Perturb Flow Variable 8

 By Perturbation Value: -5.000500E+03

 Average Number of Significant Digits in Matrix 6.698413E+00

 

 Solving Linear Program

 Optimal Solution Found

 Objective Value -5.302168E+04

 Maximum Relative Change in Flow Variable 8.326215E-05

 Max Relative External Variable Change 1.369594E-04

 SLP Algorithm: End Iteration 2

 SLP Algorithm: Begin Iteration 3

 Running Base Flow Process Simulation

 Status of Simulation-Based Constraints 

 Constraint Type Name Status Distance To RHS

 --------------- ---- ------ ---------------

 Stream Depletion S01.1 Satisfied 4.4826E+02

 Stream Depletion S01.2 Satisfied 5.3093E+03

 Stream Depletion S01.3 Satisfied 5.8128E+03

 Stream Depletion S01.4 Satisfied 5.8671E+03

 Stream Depletion S01.5 Near-Binding 5.3802E-03

 Stream Depletion S01.6 Satisfied 5.2861E+03

 Stream Depletion S01.7 Satisfied 5.8116E+03

 Stream Depletion S01.8 Satisfied 5.8670E+03

 Stream Depletion S01.9 Near-Binding 1.2965E-05

 Stream Depletion S01.10 Satisfied 5.2861E+03

 Stream Depletion S01.11 Satisfied 5.8116E+03

 Stream Depletion S01.12 Satisfied 5.8670E+03

 Stream Depletion S02.3 Satisfied 1.9513E+00

 Stream Depletion S02.4 Near-Binding 1.4319E-01

 Stream Depletion S02.7 Near-Binding 1.5495E-04

 Stream Depletion S02.8 Near-Binding 3.8380E-04

 Stream Depletion S02.11 Near-Binding 1.2510E-04

 Stream Depletion S02.12 Near-Binding 3.8262E-04

 Stream Depletion S03.1 Satisfied 1.4120E+04

 Stream Depletion S03.2 Satisfied 3.6890E+00

 Stream Depletion S03.3 Satisfied 1.4209E+04

 Stream Depletion S03.4 Satisfied 1.2081E+04

 Stream Depletion S03.5 Satisfied 1.4032E+04

 Stream Depletion S03.6 Near-Binding 1.0418E-03

 Stream Depletion S03.7 Satisfied 1.4209E+04

 Stream Depletion S03.8 Satisfied 1.2081E+04

 Stream Depletion S03.9 Satisfied 1.4032E+04

 Stream Depletion S03.10 Near-Binding 9.9406E-04

 Stream Depletion S03.11 Satisfied 1.4209E+04

 Stream Depletion S03.12 Satisfied 1.2081E+04
92 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

 Stream Depletion S04.3 Satisfied 6.4462E+03

 Stream Depletion S04.4 Near-Binding 2.0796E-01

 Stream Depletion S04.7 Satisfied 6.4438E+03

 Stream Depletion S04.8 Near-Binding 4.1881E-04

 Stream Depletion S04.11 Satisfied 6.4438E+03

 Stream Depletion S04.12 Near-Binding 4.1881E-04

 

 Calculating Response Matrix

 Perturb Flow Variable 1

 By Perturbation Value: -2.500750E+03

 Perturb Flow Variable 2

 By Perturbation Value: -2.500750E+03

 Perturb Flow Variable 3

 By Perturbation Value: -2.500750E+03

 Perturb Flow Variable 4

 By Perturbation Value: -2.500750E+03

 Perturb Flow Variable 5

 By Perturbation Value: -2.500750E+03

 Perturb Flow Variable 6

 By Perturbation Value: -2.500750E+03

 Perturb Flow Variable 7

 By Perturbation Value: -2.500750E+03

 Perturb Flow Variable 8

 By Perturbation Value: -2.500750E+03

 Average Number of Significant Digits in Matrix 6.400794E+00

 

 Solving Linear Program

 Optimal Solution Found

 Objective Value -5.302168E+04

 Maximum Relative Change in Flow Variable 7.157547E-08

 Max Relative External Variable Change 3.980410E-08

 SLP Algorithm: End Iteration 3

 Iterations have converged

----------------------------------------------------------------------

 Ground-Water Management Solution

----------------------------------------------------------------------

 OPTIMAL SOLUTION FOUND 

 OPTIMAL RATES FOR EACH FLOW VARIABLE 

 ---------------------------------------

Variable Withdrawal Injection Contribution

Name Rate Rate To Objective

---------- -------------- ------------ ------------

 Q1 3.483296E+03 3.814209E+03

 Q2a 5.000000E+04 1.380000E+04

 Q2b 2.755393E+04 7.522222E+03

 Q2c 2.690226E+04 7.344317E+03

 Q2d 2.673937E+04 7.299849E+03

 Q3 2.829550E-06 1.032786E-06

 Q4a 4.896278E+04 1.336684E+04

 Q4b 8.068633E+03 2.202737E+03

 ------------ ------------ ------------

TOTALS 1.917103E+05 0.000000E+00 5.535017E+04
Sample Problems 93

 OPTIMAL RATES FOR EACH EXTERNAL VARIABLE 

 ---------------------------------------

Variable Export Import Contribution

Name Rate Rate To Objective

---------- -------------- ------------ ------------

 Im9 0.000000E+00 0.000000E+00

 Im10 0.000000E+00 0.000000E+00

 Im11 1.461444E+04 -1.595897E+03

 Im12 6.708699E+03 -7.325900E+02

 ------------ ------------ ------------

TOTALS 0.000000E+00 2.132314E+04 -2.328487E+03

 OBJECTIVE FUNCTION VALUE 5.302168E+04

 BINDING CONSTRAINTS 

Constraint Type Name Status Shadow Price

--------------- ---- ------ ------------

Stream Depletion S01.9 Binding -3.2407E-01

Stream Depletion S02.11 Binding -5.9427E-01

Stream Depletion S02.12 Binding -1.1097E-01

Stream Depletion S03.10 Binding -5.3048E-02

Stream Depletion S04.8 Binding -1.2527E-01

Stream Depletion S04.12 Binding -3.0472E-01

Summation p06u Binding 2.5658E-01

Summation p11l Binding -1.0920E-01

Summation p12l Binding -1.0920E-01

Maximum Flow Rate Q2a Binding Not Available

 Binding constraint values are determined from the linear program

 and based on the response matrix approximation of the flow-process.

 RANGE ANALYSIS NOT REPORTED 

----------------------------------------------------------------------

 Final Flow Process Simulation

----------------------------------------------------------------------

 Status of Simulation-Based Constraints 

 Using Optimal Flow Rate Variable Values

 Constraint Type Name Status Distance To RHS

 --------------- ---- ------ ---------------

 Stream Depletion S01.1 Satisfied 4.4826E+02

 Stream Depletion S01.2 Satisfied 5.3093E+03

 Stream Depletion S01.3 Satisfied 5.8128E+03

 Stream Depletion S01.4 Satisfied 5.8671E+03

 Stream Depletion S01.5 Near-Binding 5.3803E-03

 Stream Depletion S01.6 Satisfied 5.2861E+03

 Stream Depletion S01.7 Satisfied 5.8116E+03

 Stream Depletion S01.8 Satisfied 5.8670E+03

 Stream Depletion S01.9 Near-Binding 1.2851E-05

 Stream Depletion S01.10 Satisfied 5.2861E+03

 Stream Depletion S01.11 Satisfied 5.8116E+03

 Stream Depletion S01.12 Satisfied 5.8670E+03

 Stream Depletion S02.3 Satisfied 1.9514E+00

 Stream Depletion S02.4 Near-Binding 1.4320E-01

 Stream Depletion S02.7 Near-Binding 2.4777E-04

 Stream Depletion S02.8 Near-Binding 3.8769E-04

 Stream Depletion S02.11 Near-Binding 2.1793E-04

 Stream Depletion S02.12 Near-Binding 3.8650E-04

 Stream Depletion S03.1 Satisfied 1.4120E+04

 Stream Depletion S03.2 Satisfied 3.6880E+00
94 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

 Stream Depletion S03.3 Satisfied 1.4209E+04

 Stream Depletion S03.4 Satisfied 1.2081E+04

 Stream Depletion S03.5 Satisfied 1.4032E+04

 Stream Depletion S03.6 Near-Binding 4.3092E-05

 Stream Depletion S03.7 Satisfied 1.4209E+04

 Stream Depletion S03.8 Satisfied 1.2081E+04

 Stream Depletion S03.9 Satisfied 1.4032E+04

 Stream Depletion S03.10 Near-Binding 9.0869E-05

 Stream Depletion S03.11 Satisfied 1.4209E+04

 Stream Depletion S03.12 Satisfied 1.2081E+04

 Stream Depletion S04.3 Satisfied 6.4462E+03

 Stream Depletion S04.4 Near-Binding 2.0796E-01

 Stream Depletion S04.7 Satisfied 6.4438E+03

 Stream Depletion S04.8 Near-Binding 4.2327E-04

 Stream Depletion S04.11 Satisfied 6.4438E+03

 Stream Depletion S04.12 Near-Binding 4.2327E-04

 Because of precision limitations and possible nonlinear behavior, 

 the status of binding constraints computed directly by the flow process 

 may differ slightly from those computed using the linear program. 

Acknowledgments

The authors thank Arlen Harbaugh of the USGS for his advice and assistance during 

the development of GWM. We also thank Eric Reichard and Mary Ashman of the USGS 

for their helpful review comments on a draft of this report. Several colleagues made 

contributions to GWM and its precursor code, MODOFC, while at the University of 

Connecticut and University of Massachusetts: R. Guy Riefler wrote early versions of the 

branch and bound algorithm, procedures for testing the precision of the response-matrix 

coefficients, and methods for handling multiple stress periods; Michele Minihane contributed to the development of the sequential linear programming algorithm; Yi-Chang Ethan 

Lin contributed to the coding and testing of GWM; and Gemma Baro-Montes did extensive 

testing of GWM. 

References Cited

Ahlfeld, D.P., and Heidari, Manoutchehr, 1994, Application of optimal hydraulic control 

to ground-water systems: Journal of Water Resources Planning and Management, v. 120, 

no. 3, p. 350–365.

Ahlfeld, D.P., and Mulligan, A.E., 2000, Optimal management of flow in groundwater 

systems: San Diego, CA, Academic Press, 185 p.

Ahlfeld, D.P., and Riefler, Guy, 2003, Documentation for MODOFC—A program for 

solving optimal flow control problems based on MODFLOW simulation, version 2.3: 

Amherst, MA, University of Massachusetts, Department of Civil and Environmental 

Engineering, 48 p.

Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Du Croz, J., 

Greenbaum, A., Hammarling, S., McKenney, A., and Sorensen, D., 1999, LAPACK 

users’ guide (3d ed.): Philadelpha, PA, Society for Industrial and Applied Mathematics, 

407 p.

Bradley, S.P., Hax, A.C., and Magnanti, T.L., 1977, Applied mathematical programming: 

Reading, MA, Addison-Wesley, 539 p.

Dantzig, G. B., 1963, Linear programming and extensions: Princeton, NJ, Princeton 

University Press, 627 p.
References Cited 95

Gass, Saul, 1985, Linear programming—Methods and applications: New York, McGrawHill, 532 p.

Gill, P.E., Murray, Walter, and Wright, M.H., 1981, Practical optimization: New York, 

Academic Press, 401 p.

Gorelick, S.M, 1983, A review of distributed parameter groundwater management 

modeling methods: Water Resources Research, v. 19, no. 2, p. 305–319.

Gorelick, S.M., Freeze, R.A., Donohue, David, and Keely, J.F., 1993, Groundwater 

contamination—Optimal capture and containment: Boca Raton, FL, Lewis Publishers, 

385 p.

Greenwald, R.M., 1998, Documentation and user’s guide—MODMAN, an optimization 

module for MODFLOW, version 4.0: Freehold, NJ, HSI GeoTrans, 112 p.

Halford, K.J., and Hanson, R.T., 2002, User guide for the drawdown-limited, multi-node 

well (MNW) package for the U.S. Geological Survey’s modular three-dimensional 

finite-difference ground-water flow model, versions MODFLOW-96 and MODFLOW2000: U.S. Geological Survey Open-File Report 02-293, 33 p.

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

the U.S. Geological Survey modular ground-water model—User guide to modularization 

concepts and the ground-water flow process: U.S. Geological Survey Open-File Report 

00-92, 121 p.

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.

Hillier, F.S., and Lieberman, G.J., 2001, Introduction to operations research (5th ed.): New 

York, McGraw-Hill Publishing, 954 p. 

Lefkoff, L.J., and Gorelick, S.M., 1987, AQMAN-Linear and quadratic programming 

matrix generator using two-dimensional ground-water flow simulation for aquifer 

management modeling: U.S. Geological Survey Water-Resources Investigations Report 

87-4061, 164 p.

McDonald, M.G., and Harbaugh, A.W., 1988, A modular three-dimensional finitedifference ground-water flow model: U.S. Geological Survey Techniques of WaterResources Investigations, book 6, chap. A1, 586 p.

Minihane, M.R., 2002, Stream-aquifer interactions and a water supply management model: 

Amherst, MA, University of Massachusetts Department of Civil and Environmental 

Engineering, Masters Thesis, 161 p.

Nemhauser, G. L., and Wolsey, L.A., 1988, Integer and combinatorial optimization: New 

York, Wiley, 763 p. 

Nishikawa, Tracy, 1998, Water-resources optimization model for Santa Barbara, 

California: Journal of Water Resources Planning and Management, v. 124, no. 5, 

p. 252–263.

Peralta, R.C., 2004, Simulation/Optimization MOdeling System (SOMOS): Logan, UT, 

Utah State University Department of Biological and Irrigation Engineering, accessed 

July 28, 2004, at http://www.engineering.usu.edu/bie/

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.

Prudic, D.E., Konikow, L.F., and Banta, E.R., 2004, A new streamflow-routing (SFR1) 

package to simulate stream-aquifer interaction with MODFLOW-2000: U.S. Geological 

Survey Open-File Report 2004-1042, 95 p.

Reichard, E.G., 1995, Groundwater-surface water management with stochastic surface 

water supplies—A simulation optimization approach: Water Resources Research, v. 31, 

no. 11, p. 2845–2865.
96 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Reichard, E.G., Land, Michael, Crawford, S.M., Johnson, Tyler, Everett, R.R., Kulshan, 

T.V., Ponti, D.J., Halford, K.J., Johnson, T.A., Paybins, K.S., and Nishikawa, Tracy, 

2003, Geohydrology, geochemistry, and ground-water simulation-optimization of the 

Central and West Coast Basins, Los Angeles County, California: U.S. Geological Survey 

Water-Resources Investigations Report 03-4065, 184 p.

Riefler, R.G., and Ahlfeld, D.P., 1996, The impact of numerical precision on the solution 

of confined and unconfined optimal hydraulic control problems: Hazardous Waste and 

Hazardous Materials, v. 13, no. 2, p. 167–176.

Wagner, B.J., 1995, Recent advances in simulation-optimization groundwater management 

modeling: Reviews of Geophysics, Supplement, U.S. National Report to International 

Union of Geodesy and Geophysics 1991–1994, p. 1021–1028.

Willis, R., and Yeh, W.W.-G., 1987, Groundwater systems planning and management: 

Englewood Cliffs, NJ, Prentice-Hall, 416 p.

Yeh, W.W.-G., 1992, Systems analysis in ground-water planning and management: Journal 

of Water Resources Management and Planning, v. 118, no. 3, p. 224–237.

Zheng, Chunmiao, and Wang, P.P., 2002, MGO—A Modular Groundwater Optimizer 

incorporating MODFLOW and MT3DMS, Documentation and user’s guide: 

Tuscaloosa, AL, The University of Alabama and Groundwater Systems Research Ltd., 

118 p.
Appendix Contents

Appendix 1. Examples of Decision Variables that can be Defined in GWM . . . . . . . . . . . . . . . . . . . . . . . 99

Appendix 2. GWM Solution Algorithms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

Linear Program Solution Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

Simplex Algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

Algorithm Specifics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

Branch and Bound Solution Algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

Appendix 3. Programmers’ Guide to GWM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

Procedures Used by GWM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

GWM Packages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

Adding New GWM Packages. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

Some Programming Details. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

Figures

1-1. Hypothetical model grid illustrating some of the types of flow-rate and external 

decision variables that can be defined in GWM: A, four flow-rate decision 

variables; and B, four flow-rate decision variables and an external decision 

variable that is a source of water to the management model . . . . . . . . . . . . . . . . . . . . . . . . . . .102

3-1. Simplified flow chart showing Ground-Water Flow and Ground-Water Management 

Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .120

Tables

3-1. Outline of packages, procedures, and modules for GWM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .122

APPENDIX 1: Examples of Decision 

Variables that can be Defined in GWM

Appendix 1 101

This section provides examples of how decision variables can be defined in a GWM 

problem. Included in the examples are discussions of how characteristics of the decision 

variables are then specified by input variables in the decision variable (DECVAR) input file 

to GWM.

Example 1: Flow-rate decision variables for a steady-state model: Figure 1-1A

shows a small model grid that consists of 3 layers, 6 rows, and 4 columns; the model consists of a single stress period. Four flow-rate decision variables have been defined within 

the model grid. Decision variable Q1 is a withdrawal well located at model cell (1, 6, 4). 

Decision variable Q2 is a withdrawal well that extends over two adjacent vertical cells (1, 

6, 2) and (2, 6, 2); transmissivity of the aquifer is such that it can be assumed that each cell 

contributes 50 percent to the well’s discharge. Decision variable Q3 is a withdrawal well 

that extends over two vertical cells that are not adjacent, (1, 4, 4) and (3, 4, 4); it can be 

assumed that 75 percent of the well’s discharge is from the upper layer and 25 percent of 

the well’s discharge is from the lower layer. Decision variable Q4 is a recharge basin 

extending over five cells in the top layer of the model (row 1, columns 2-4 and row 2, 

columns 3-4); recharge is distributed evenly to the five cells.

Records 3a and 3b of the DECVAR file would have the following information for 

these four decision variables:

3a. FVNAME NC LAY ROW COL FTYPE FSTAT WSP

3b. RATIO LAY ROW COL

Q1 1 1 6 4 W Y 1

Q2 2 0 0 0 W Y 1

0.5 1 6 2

0.5 2 6 2

Q3 2 0 0 0 W Y 1

0.75 1 4 4

0.25 3 4 4

Q4 5 0 0 0 I Y 1

0.20 1 1 2

0.20 1 1 3

0.20 1 1 4

0.20 1 2 3

0.20 1 2 4

Example 2: Flow-rate decision variables for a transient model: The model grid 

and locations of the wells and recharge basin are the same as for example 1, but the model 

is now transient and consists of four stress periods that represent the four seasons of the year 

(Winter, Spring, Summer, and Fall). The withdrawal well at model cell (1, 6, 4) can have 

different pumping rates during each of the four stress periods, so four decision variables 

now must be used for the well site; these are specified as decision variables Q1-W, Q1-Sp, 

Q1-Su, and Q1-F. The withdrawal rates at the second and third well sites (decision variables Q2 and Q3) must be constant throughout the year. The recharge basin is used only 

during the winter and spring, and the recharge rates to the basin during the two seasons need 

not be the same. Two decision variables for recharge therefore are necessary for the basin, 

and these are specified as Q4-W and Q4-Sp.
102 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

COLUMNS

ROWS

LAYERS

1

2

3

12 3 4

1

2

3

4

5

6

A.

Q4

Q3

Q2 Q1

COLUMNS

ROWS

LAYERS

1

2

3

12 3 4

1

2

3

4

5

6

B.

Q4

Q3

Q2 Q1

IM1

RESERVOIR

Figure 1-1. Hypothetical model grid illustrating some of the types of flow-rate and external 

decision variables that can be defined in GWM: A, four flow-rate decision variables (Q1, Q2, Q3, 

and Q4); and B, four flow-rate decision variables and an external decision variable (a surfacewater reservoir, IM1) that is a source of water to the management model.
Appendix 1 103

There are now a total of eight decision variables for the management model, and 

records 3a and 3b of the DECVAR file would have the following information for the 8 

variables:

3a. FVNAME NC LAY ROW COL FTYPE FSTAT WSP

3b. RATIO LAY ROW COL

Q1-W 1 1 6 4 W Y 1

Q1-Sp 1 1 6 4 W Y 2

Q1-Su 1 1 6 4 W Y 3

Q1-F 1 1 6 4 W Y 4

Q2 2 0 0 0 W Y 1-4

0.5 1 6 2

0.5 2 6 2

Q3 2 0 0 0 W Y 1-4

0.75 1 4 4

0.25 3 4 4

Q4-W 5 0 0 0 I Y 1

0.20 1 1 2

0.20 1 1 3

0.20 1 1 4

0.20 1 2 3

0.20 1 2 4

Q4-Sp 5 0 0 0 I Y 2

0.20 1 1 2

0.20 1 1 3

0.20 1 1 4

0.20 1 2 3

0.20 1 2 4
104 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Example 3: Flow-rate and binary decision variables for a transient model: The 

model grid and all of the flow-rate decision variables for this example are the same as for 

example 2, but binary variables now are associated with the flow-rate decision variables. 

These binary variables are used to define the flow-rate decision variable as either active or 

inactive; this might be necessary, for example, when there is a minimum amount of withdrawal or injection that is required at each of the decision variables before the well or 

recharge basin will be constructed. The decision variables at well site (1, 6, 4) are associated with a single binary variable BV1, because if withdrawals occur at the well in at least 

one of the seasons, the well site must be constructed. The same situation holds for the 

recharge basin, so binary variable BV4 is associated with flow-rate variables Q4-W and 

Q4-Sp.

Records 3a and 3b of the DECVAR file would have the same information for the four 

flow-rate decision variables as given in example 2. Characteristics of the binary variables 

are specified in record 5 of the DECVAR file:

5. BVNAME NDV BVLIST

BV1 4 Q1-W Q1-Sp Q1-Su Q1-F

BV2 1 Q2

BV3 1 Q3

BV4 2 Q4-W Q4-Sp

Example 4: Flow-rate and external decision variables for a transient 

model: The model grid and all of the flow-rate decision variables for this example are the 

same as for the previous example, but there is now an external source of water that is available to the management problem from a surface-water reservoir shown schematically in 

figure 1-1B. The rate at which this source of water can be imported can vary during each 

of the four stress periods; therefore, four external decision variables are needed to model 

the source: IM1-W, IM1-Sp, IM1-Su, and IM1-F. Note that these external decision variables are not associated with any of the model grid cells. Characteristics of the external 

decision variable are specified in record 4 of the DECVAR file: 

4. EVNAME ETYPE ESP

IM1-W IM 1

IM1-Sp IM 2

IM1-Su IM 3

IM1-F IM 4
APPENDIX 2: GWM Solution Algorithms

Appendix 2 107

GWM includes the RMS (Response Matrix Solution) Package for solving optimization 

problems. This package contains three solution algorithms: linear programming, branch 

and bound, and sequential linear programming. Linear programming is used to solve problems with linear objective and constraint functions. When binary variables are present, the 

branch and bound algorithm is used. This algorithm is based on solving a series of linear 

programs. If the objective or constraints include nonlinear responses of system state to 

flow-rate variables, then the sequential linear programming (SLP) algorithm is used. SLP 

operates by solving a series of linear programs. Because solution of a linear program forms 

the core of all three solution algorithms, it is the primary focus of this section.

Linear Program Solution Algorithm

Simplex Algorithm

GWM1SIMPLEX1 is a collection of FORTRAN subroutines that solve a linear program (LP) by using a 2-phase revised simplex method (Dantzig, 1963; Gass, 1985). Sensitivity (or range) analysis is also done by GWM1SIMPLEX1 if an optimal solution to the 

LP exists. GWM1SIMPLEX1 assumes that the problem takes the following form:

Minimize (2-1)

subject to (2-2)

, (2-3)

where is the vector of decision variables, is the transposed column vector of objectivefunction coefficients, is an m x n matrix where m is the number of constraints and n is 

the number of decision variables (including slack variables), is the total amount of 

resource available, and is the upper-bound vector on the decision variables. Other LP formulations can be transformed into equations 2-1 through 2-3 through fairly simple means. 

For example, a maximization objective can be transformed into an equivalent minimization 

objective by multiplying all objective-function coefficients by (-1). The transformation 

from inequality constraint to equality constraint is accomplished by adding one variable to 

each inequality constraint. These extra variables indicate the difference between the amount 

of resource used by a particular solution and the amount of resource available. In that context, the variables represent the slack or surplus of a resource, and are called slack or surplus 

variables.

The original LP in equations 2-1 through 2-3 represents a system of n unknowns in m

equations, where n is greater than m. The simplex method proceeds by fixing (n-m) decision 

variables and using the equality constraints to solve for the remaining m variables. At any 

iteration, the variables with fixed values are termed nonbasic variables, , while the 

Z c

t = x

Ax = b

0 ≤ ≤ x u

x c

t

A

b

u

xN
108 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

remaining m variables are called basic variables, . By using these variables and defining 

the matrices and as the columns of corresponding to the basic and nonbasic variables, respectively, the original LP can be written as:

Minimize (2-4)

subject to (2-5)

(2-6)

, (2-7)

where , and , are the objective-function coefficients and upper bounds for the 

basic and nonbasic variables, respectively. Assuming that the list of basic variables is 

known and that the nonbasic variables are assigned to specific values, then the values of the 

basic variables are determined by rearranging equation 2-5:

. (2-8)

Furthermore, substituting equation 2-8 into 2-4, and rearranging, the objective function is 

calculated as

Minimize . (2-9)

The rules for identifying the set of basic variables and assigning values to the nonbasics are 

a key component of the simplex method.

The simplex method proceeds by finding a sequence of feasible solutions to the system 

of equations in equations 2-1 through 2-3 until an optimal solution is found. At each iteration, the nonbasic variables are considered fixed at either their upper or lower bound. The 

basic variables are then computed from equation 2-8 to ensure feasibility. Information from 

equation 2-9 is used to determine if a given solution is optimal. If the current solution is not 

optimal, then a new simplex iteration begins and one basic variable becomes nonbasic 

while one nonbasic variable becomes basic.

The nonbasic variable that becomes basic (that is, enters the basis) is chosen on the 

basis of its impact on the objective function. The objective function, as written in equation 

2-9, includes a term that only depends on and a constant term. The coefficient on in 

equation 2-9 is called the reduced cost. The nonbasic variable that will have the most 

advantageous impact on the objective function, that is, the best reduced cost, is chosen as 

the entering nonbasic variable. The basic variable that first goes to either its upper or lower 

bound as the nonbasic variable is changed, leaves the basis, and becomes nonbasic.

The simplex method finds a sequence of feasible solutions by altering the set of basic 

variables one variable at a time until the reduced costs indicate that no further improvement 

is possible. At this point, the optimal solution has been found. In many LPs, it is a nontrivial 

problem to find an initial feasible solution from which to begin the search for an optimal 

solution. This problem is addressed in GWM1SIMPLEX1 by use of the 2-phase revised 

simplex method, which proceeds in Phase I by altering the original problem (eqs. 2-1 

through 2-3) in such a way that an initial feasible solution to the altered problem is trivial 

and the optimal solution of the Phase I problem is a feasible solution to the original 

problem.

xB

B N A

Z cB

t

xB cN

t

x = + N

BxB NxN + = b

0 xB ≤ ≤ uB

0 xN ≤ ≤ uN

cB uB cN uN

xB B–1 b NxN = ( ) –

Z cB

t

B–1b cN

t

 cB

t

B–1 ( ) + N x = + N

xN xN
Appendix 2 109

Algorithm Specifics

Input to GWM1SIMPLEX1: GWM1SIMPLEX1 is a self-contained implementation of the 2-phase revised simplex method. Augmentation of the original LP with slack 

variables so that the problem takes the form in equations 2-1 through 2-3 is required prior 

to calling GWM1SIMPLEX1. Variables related to solving the LP are defined in the module 

GWM1RMS1, whereas variables used exclusively by GWM1SIMPLEX1 and its subroutines are defined in module GWM1RMS1LP. Calls to GWM1SIMPLEX1 should be of the 

form

CALL GWM1SIMPLEX1(M, NV, NDV, AMAT, COST, BNDS, RHS, OBJ, IFLG,

LPITMAX),

where

M is the number of constraints;

NV is the number of variables (decision variables + slacks);

NDV is the number of decision variables + 1;

AMAT(M,NDV) is the original matrix augmented with slack-variable coefficients in 

the last column;

COST(NV) is the vector of objective-function coefficients for each variable, 

including slacks;

BNDS(NV) is the array of upper bounds on each variable, including slacks;

RHS(M) is the original right-hand side values of the constraints; and

LPITMAX is the maximum number of iterations allowed.

If an optimal solution to the original LP is found, GWM1SIMPLEX1 stores the optimal solution, shadow prices, reduced costs, and range-analysis information in space that 

contained the original LP information, so that output from GWM1SIMPLEX1 includes

IFLG the output-status flag indicating optimal solution found (0), problem is 

infeasible (1), or problem is unbounded (2);

COST(NV) the values of the decision variables at the optimal solution;

BNDS(NV) the reduced costs of the decision variables;

RHS(M) the shadow prices of the constraints (that is, the dual variables); and

OBJ the value of the objective function at the optimal solution.

If specified by the user, range analysis is done on the objective-function coefficients 

( ) and the right-hand-side values (RHS, ). During range analysis, one parameter at a 

time is perturbed about its original value until the optimal basis changes. When this occurs, 

the variables that would enter and leave the basis are determined by means of a ratio test. 

For the objective-function coefficient ranges, the optimal values of the decision variables 

will not change provided the coefficients remain within the stated ranges. Conversely, 

when an RHS value is modified, the value of the decision variables will change (see eq. 2-

8), but the optimal basis and shadow prices will not be affected. Upon completion of range 

analysis, program control is then passed back to the calling program.

A

c b
110 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Phase I and II: The Phase I procedure is needed to obtain an initial feasible solution 

to the LP. Phase I begins by augmenting the original problem (eqs. 2-1 through 2-3) with 

artificial variables, , such that the artificial LP becomes

Minimize (2-10)

subject to (2-11)

(2-12)

, (2-13)

where is an (m x 1) vector of ones and is an (m x m) identity matrix. The Phase I problem has the form of a linear program. An initial feasible solution to this artificial problem 

is obtained by setting all original variables and using the system of equality equations 

in the constraints to solve for the values of (see Bradley and others, 1977, for more detail 

about the 2-phase simplex method). Because the objective-function coefficients on the original variables have been set to zero in the artificial problem, application of the simplex algorithm to equations 2-10 through 2-13 will tend to drive the artificial variables from the 

problem.

If a feasible solution to the original problem exists, then all artificial variables will be 

driven to zero and the optimal solution in Phase I represents a feasible solution to the original problem. The original LP is then solved during Phase II from the initial feasible solution obtained from the artificial LP. If the solution to equations 2-10 through 2-13 retains 

any artificial variables, the original LP must be infeasible and the algorithm terminates.

Pricing: As discussed above, each iteration of the simplex method proceeds by substituting one basic variable with one nonbasic variable. The decision process in this step 

requires first determining which nonbasic variable will enter the basis. This decision is 

made during the pricing step. The reduced cost of all nonbasic variables is determined and 

the entering variable is chosen as the one that results in the largest marginal improvement 

in the objective function:

 

. (2-14)

When the nonbasic variable is at its lower bound, the most negative reduced cost is advantageous. Conversely, when the nonbasic variable is at its upper bound, the most positive 

reduced cost is advantageous.

The Ratio Test: Once the entering variable is identified in the pricing step, the leaving variable must be determined. At each iteration, the set of equality constraints in equations 2-1 through 2-3 must be satisfied to maintain feasibility. Thus, the constraints are used 

to determine which variable will leave the basis. When the entering variable increases 

from its lower bound (that is, it is currently equal to zero and will increase in value upon 

entering the basis), the ratio test

x'

Z 1

t

x' 0

t = + x

Ax + Ix' = b

0 ≤ ≤ x u

0 ≤ x'

1 I

x = 0

x'

Min

cN cB

t

B–1 ( ) – N xN ∀ = 0

cN cB

t

B–1 –( ) – N xN ⎩ ∀ = u

⎨

⎧

xj
Appendix 2 111

(2-15)

is done, where i and j refer to elements in the associated vector or matrix. If one of the first 

two ratios provides the limiting value, then the variable leaves the basis. If the third argument, , is the minimum, however, then the entering variable hits its upper bound prior to 

a basic variable hitting one of its bounds. Consequently, the variable remains nonbasic but 

its value is set to its upper bound.

Similarly, when the entering variable enters from its upper bound and decreases in 

value upon entering the basis, the next iteration is determined by

 

.

(2-16)

Once the entering and leaving variables are determined, a pivot operation is done in 

which the basis pointer array is updated and the basis matrix is updated and factorized. The 

process of substituting the appropriate columns of into the matrices and is termed 

pivoting. After each pivot operation, the basis matrix is factorized by the method of LU 

decomposition, by use of routines from LAPACK (Anderson and others, 1999).

The steps for determining the entering and leaving variables followed by pivoting and 

factorization are done iteratively until one of three stopping criteria is met:

1. An optimal solution is found. This is determined during the pricing operation: if all 

prices are greater than or equal to zero for nonbasics at their lower bounds or less than 

or equal to zero for nonbasics at their upper bounds, then the solution cannot be 

improved and the optimal solution has been found. 

2. The iteration limit, LPITMAX, is exceeded. When this occurs, an error message is 

printed and GWM is terminated.

3. The problem is determined to be unbounded. If a leaving variable cannot be found in 

the ratio test, then the problem is unbounded.

Range analysis: When an optimal solution is found, range analysis can be performed. The goal during range analysis is to determine the range over which a specified 

parameter can vary without changing the optimal basis. At the limit of these ranges, the new 

optimal basis is also determined. Range analysis is applicable for strictly linear programs 

and should be used with caution if the problem is a mixed-binary linear program or is nonlinear but being solved by using sequential linear programming.

Min

B–1 b NxN [ ] ( ) – i

B–1 ( ) N ij

----------------------------------- B–1 ( ) N ij 0 >

uB, i

 B–1 b NxN [ ] ( ) – i –

 B–1 – ( ) N ij

---------------------------------------------------- B–1 ( ) N ij 0 <

uj ⎩

⎪

⎪

⎪

⎨

⎪

⎪

⎪

⎧

xi

uj

Min

B–1 b NxN [ ] ( ) – i

 B–1 – ( ) N ij

----------------------------------- B–1 ( ) N ij 0 <

uB, i

 B–1 b NxN [ ] ( ) – i –

B–1 ( ) N ij

---------------------------------------------------- B–1 ( ) N ij 0 >

uj ⎩

⎪

⎪

⎪

⎨

⎪

⎪

⎪

⎧

A B N
112 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

Objective-function coefficient ranges are determined by considering the value of the 

reduced cost for each variable that would cause the variable to enter the basis. For nonbasic 

variables, objective-function coefficients can vary in the range

 . (2-17)

When the range limit is exceeded, xj becomes a candidate to enter the basis. 

For basic variables, the objective-function coefficient values that would cause the variable to leave the basis are determined by imposing the optimality condition on all reduced 

costs. Letting the nonbasic variable be variable j, the objective-function coefficient range 

for the basic variable i is determined by

 

. (2-18)

The nonbasic variable producing the limiting ratio in each case enters the basis. For all 

objective-function coefficient ranges, the leaving variables are determined by the 

standard ratio test, as described above. Because the optimal solution is calculated as 

, the optimal solution is unaffected by objective-function coefficients that 

vary within the ranges determined above.

Right-hand-side range analysis is done by determining the range of RHS values for 

which the optimal basis remains unchanged. Note, however, that the optimal values of the 

basic variables (given by eq. 2-8) will change with a change in the RHS. RHS range analysis proceeds by first determining the leaving variables as those that provide the limiting 

ratios

,

(2-19)

where

is the change in the kth RHS (k = 1, ..., m);

is the optimal value of the ith basic variable;

is the ith row, kth column element of ; and

is the upper bound of the ith basic variable.

cj cB

t

B–1 ( ) N j – [ ] – Δcj ≤ ≤ ∞ xN ∀ = 0

– ∞ Δcj

 cj cB

t

B–1 ( ) N j ≤ ≤ – [ ] – xN ⎩ ∀ = u

⎨

⎧

Max

j

 cj

 cBB–1 ( ) A j ( ) –

B–1 ( ) A ij

----------------------------------- for

B–1 ( ) A ij < 0 xj , 0 =

B–1 ( ) A ij > 0 xj , = u

 Δc ≤ i

 Min 

j

cj

 cBB–1 ( ) A j ( ) –

B–1 ( ) A ij

----------------------------------- for

B–1 ( ) A ij > 0 xj , 0 =

B–1 ( ) A ij < 0 xj , = u

≤

xB B–1 b NxN = ( ) –

Max

i

xB ( )i 

βik

– -----------βik 

> 0

uB ( )i

 xB ( )i –

βik

----------------------------βik 

< 0

 Δbk Min

i

xB ( )i 

βik

– -----------βik 

< 0

uB ( )i

 xB ( )i –

βik

---------------------------- βik 

> 0

≤ ≤

Δbk

xB ( )i B–1 b NxN ( ) ( ) – = i

βik B–1

uB ( )i
Appendix 2 113

The entering variable is determined by the optimality condition, which requires that the 

reduced costs of nonbasic variables at their lower bounds remain non-negative while the 

reduced costs of nonbasics at their upper bounds stay non-positive. The nonbasic variable 

that first violates this condition becomes the candidate to enter the basis. Let the leaving 

variable identified above be basic in row i. If the leaving variable leaves at its lower bound, 

then the entering variable provides the limiting value in the following test:

. (2-20)

Likewise, if the leaving variable leaves at its upper bound, the entering variable is determined by

 . (2-21)

Output from GWM1SIMPLEX1: Output from GWM1SIMPLEX1 consists of the 

optimal solution, shadow prices, reduced costs, and range-analysis information. Shadow 

prices represent the sensitivity of the objective-function value to changes in the right-handside value of a particular constraint [that is, ] (Bradley and others, 1977; Ahlfeld 

and Mulligan, 2000; Hillier and Lieberman, 2001). Shadow prices also represent the optimal value of the dual variables. The slack, or the difference between the right-hand side ( ) 

and the constraint value at the optimum, indicates the amount of a resource not utilized at 

the optimal solution. 

Branch and Bound Solution Algorithm

The branch and bound algorithm is used by GWM to solve problems that include 

binary variables. The algorithm consists of solving a series of linear programs. Within each 

linear program, the binary variables are either assigned, a priori, a binary value (0 or 1) or 

allowed to take a noninteger value between 0 and 1. A binary variable that is not required 

to take a binary value is referred to as relaxed. By assigning values to some binary variables 

and relaxing others, a problem is created that no longer has decision variables that are 

required to take binary values. Such a problem is referred to as a subproblem and can be 

solved by using the simplex algorithm for linear programs. By solving a series of subproblems, the branch and bound algorithm can identify the optimal solution to the original 

mixed-binary problem.

The branch and bound algorithm begins by relaxing all binary variables and solving 

the resulting linear program. If the fully relaxed problem has no feasible solution, then the 

full binary problem will also be infeasible. The solution to the fully relaxed problem provides a starting point for the creation of additional subproblems by branching.

The algorithm proceeds by selecting one of the relaxed binary variables from the initial 

subproblem and branching on it to create two new subproblems. In the first new subproblem, the selected variable is forced to take a value of 0 while the remaining binary variables 

are still relaxed. In the second new subproblem, the selected variable is forced to 1. Each 

of these subproblems is solved. Additional subproblems are created by branching on other 

relaxed binary variables. The implementation of the branch and bound algorithm in GWM 

uses a modified depth-first search with backtracking (Nemhauser and Wolsey, 1988). 

Min

xN

cj

 cB

t

B–1 ( ) N j [ ] –

B–1 ( ) N ij

– ----------------------------------- B–1 ( ) N ij < 0 xN, j , = 0 B–1 ( ) N ij > 0 xN, j u = j ; ,

⎩

⎨

⎧

Min

xN

cj cB

t

B–1 ( ) N j [ ] –

B–1 ( ) N ij

----------------------------------- B–1 ( ) N ij > 0 xN, j , = 0 B–1 ( ) N ij < 0 xN, j u = j ; ,

⎩

⎨

⎧

( ) ∂Z ⁄ ( ) ∂b

b
114 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

For each parent subproblem, two children are created, branching from the first available 

nonspecified variable. Branching proceeds down the tree along the path with binary 

variables set to one. A notable step in the algorithm occurs when a feasible linear program 

solution is found in which all binary variables have binary values. This is referred to as a 

binary solution and is a feasible, although not necessarily optimal, solution to the original 

mixed-binary problem. Once a binary solution is determined, additional branching is 

accomplished by backtracking, with the most recently generated subproblem branched 

from first.

Branching on each relaxed binary variable in each generated subproblem would eventually lead to enumeration of all possible combinations of 0/1 values for each binary variable. Even for moderately sized problems, the number of linear programs that would be 

solved can be substantial. The branch and bound algorithm avoids enumeration of all possible combinations of binary values by eliminating subproblems and all their successive 

branching possibilities by using a series of fathoming tests. If a subproblem is fathomed, 

then additional branching from that subproblem will not yield the optimal solution and the 

subproblem and all it subsidiary subproblems can be dropped from further consideration.

The fathoming test has three parts:

1. A binary solution is found and thus the subproblem is fathomed, because no additional subproblems are possible. The objective function value for this subproblem is 

compared with any other subproblem that has a fully binary solution to identify the 

current best binary solution. The objective value for the current best solution provides 

an upper bound on the objective function for the original problem. That is, the optimal 

solution to the original problem will have an objective function value that is less than 

or equal to the one found from the subproblem with the best binary solution. 

2. The subproblem is fathomed if it is infeasible. If a subproblem, with some binary variables relaxed, is infeasible, then imposing additional restrictions on binary variables 

will not produce a feasible solution. Hence, no successive subproblems can be feasible, and the subproblem is not worth further branching.

3. The subproblem is fathomed if the objective-function value for the optimal solution 

of the subproblem is larger than the upper bound from the current best binary solution. 

This applies even if the subproblem has some variables that are relaxed. It follows that 

imposition of additional restrictions on binary variables will not improve the objective function, so that all successive subproblems will have objective values that are 

inferior to the current best binary solution.

The optimal solution to the problem is found when all subproblems are fathomed. The 

subproblem with the best binary solution at that point is the solution to the original problem.

One extension to the branch and bound algorithm is used by GWM for the case in 

which some binary variables are specified with values of 1 and some relaxed binary variables are determined by the linear program to have 0 optimal values.

Experience with ground-water flow management problems has shown that it is common for candidate-well pumping rates to take 0 values at the optimal solution even when 

binary variables are not used. This characteristic results from the geometry of the groundwater simulation problem and the proximity of the flow-rate variables to constraints.
Appendix 2 115

At the completion of the solution of a subproblem, GWM checks the value of the 

relaxed binary variables and their associated non-binary variables. If the values of all associated non-binary variables are 0, then the relaxed binary variable can be treated as having 

taken a 0 value. That is, the same solution to the subproblem would be obtained if that 

binary variable were forced to take a 0 value. If this condition is true for all relaxed variables in the subproblem, then all relaxed variables are reset to redefine the subproblem as 

if the relaxed variables had been specified. Implementation of this rule will reduce the number of subproblems that require solution. In effect, redefinition of the subproblem jumps 

ahead of the branching process to a subproblem that would have been identified after many 

additional iterations.

Two GWM input parameters control the branch and bound algorithm. BBITMAX 

specifies the maximum number of subproblems that will be considered by the algorithm. 

The RMS implementation of the branch and bound algorithm stores certain information 

about each subproblem. BBITMAX serves to define the dimension of the arrays that are 

used to store this information. BBITPRT controls the output directed to the global output 

file. Full output consists of details of the progress of the algorithm at each subproblem. 

APPENDIX 3: Programmers’ Guide to GWM

Appendix 3 119

The Ground-Water Management (GWM) Process solves an optimization problem to 

minimize or maximize a specified objective function subject to specified constraints. The 

current implementation of GWM uses a numerical method for solution that is based on construction of a response matrix and its use in a linear program. GWM does not interact with 

any process, other than the Ground-Water Flow (GWF) Process, in MODFLOW. GWM is 

organized into a number of packages that define the objective function, various constraint 

types, and the solution algorithm. A number of new procedures were defined for MODFLOW that are required for completion of the GWM Process.

Procedures Used by GWM

Five new procedures have been added to MODFLOW for use by GWM. Figure 3-1 

depicts the placement of these procedures in the overall structure of MODFLOW-2000. 

The core computations of the GWF Process are carried out in the Stress, Time Step, and 

Iteration Loops. The GWM Process requires repeated runs of the GWF Process; these runs 

are controlled with two levels of loops that surround the GWF Process—the GWM Iteration 

Loop and the Flow-Process Loop. For the Response Matrix Solution (RMS) Package, the 

Flow-Process Loop performs the perturbations required to construct a response matrix, 

whereas the GWM Iteration Loop performs iterations required by the sequential linear

programming algorithm.

AR Procedure—Allocate and Read: This procedure combines the functions of the AL 

(Allocate) and RP (Read and Prepare) Procedures found elsewhere in MODFLOW. 

Because GWM packages use FORTRAN allocatable arrays, it is not necessary to 

have a separate allocation procedure. The GWM AR Procedure is located after GWF 

AL and RP Procedures in figure 3-1. 

PL Procedure—Prepare Loop: This procedure processes information from the GWM Iteration Loop and prepares information for the Flow-Process Loop. For the RMS Package, the PL Procedure sets indices and counters, sets perturbation size, and sets base 

flow-rate variable values.

PP Procedure—Prepare Process: This procedure prepares information for a single run of 

the flow process. For the RMS Package, the PP Procedure determines if the flow-process run is for a base run or a perturbation run and sets the flow-rate variables accordingly. The GWM PP Procedure is followed by a GWM RW Procedure (fig. 3-1). The 

RW Procedure rewinds appropriate GWF Process files that will be reread during the 

Flow-Process Loop.

OS Procedure—Observe State: This procedure obtains the state of the system, as simulated 

by the GWF Process, at selected locations at the completion of a stress period. The 

specific locations at which the state is observed are determined by the needs of the 

objective and constraint packages.

FP Procedure—Formulate Process: This procedure processes the results of a completed 

flow-process run. It can be used to formulate the portions of the optimization problem 

that depend upon the flow process. For the RMS Procedure, the response coefficients 

are computed as part of the FP Procedure.

Figure 3-1 concludes with three additional GWM Procedures. The FM procedure formulates the equations that will be solved. For the RMS Package, the FM procedure assembles all the elements of the linear program. The AP Procedure solves the optimization problem at an iteration. Finally, the OT Procedure writes relevant output at the solution.
120 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

GWF RP

GWM ITERATION LOOP

FLOW-PROCESS LOOP

STRESS LOOP

TIME STEP LOOP

ITERATION LOOP

GWF and GWM Procedures:

DF — Define

AL — Allocate

RP — Read and Prepare

RW — Rewind

ST — Stress

AD — Advance a time step

FM — Formulate equations

AP — Solve equations

OC — Output Control

BD — Calculate volumetric budget

OT — Write output

GWM-Specific Procedures:

AR — Allocate and read

PL — Prepare Loop

PP — Prepare Process

OS — Observe State

FP — Formulate Process

GLOBAL DF, AL, and RP

GWF AL

GWM AR

GWM PL

GWM PP

GWM RW

GWF AL and RP

GWF ST and RP

GWF AD

GWF FM

GWF AP

GWF OC, BD, OT

GWM OS

GWM FP

GWM FM

GWM AP

GWM OT

Figure 3-1. Simplified flow chart showing Ground-Water Flow (GWF) and Ground-Water Management (GWM) 

Processes.
Appendix 3 121

GWM Packages

The Ground-Water Management (GWM) Process is organized into a series of packages. Table 3-1 summarizes the packages used by GWM, the procedures used by each 

package, and primary module names assigned to each package/procedure combination.

DCV Package—Decision Variables: This package processes information regarding the 

flow-rate, external, and binary decision variables of the optimization problem. 

Because the flow-rate decision variables affect the flow process, two new GWF modules are needed; the GWF1DCV1FM and GWF1DCV1BD modules add managed 

flows to the GWF formulation and budget calculations, respectively.

OBJ Package—Objective Function: This package processes information describing the 

objective function for the optimization problem.

HDC (Head Constraints), STC (Streamflow Constraints), SMC (Summation Constraints), 

DCC (Decision-Variable Constraints) Packages—These packages all have the same 

structure. Information about the constraints is read and processed in the AR procedure. For the head and streamflow constraints, the state of the flow-process simulated 

system is recorded in the OS procedures. All constraints are placed in the GWM equations in FM procedures and constraint status is described in the OT procedures. 

RMS Package—Response Matrix Solution: This package uses the response-matrix method 

to solve the optimization problem. RMS incorporates the simplex algorithm, the 

branch and bound algorithm, and the sequential linear programming algorithm as the 

numerical solver for the response-matrix method. The response matrix consists of the 

coefficients of the responses of head and streamflow at the constraint locations to 

changes in the flow-rate decision variables. RMS controls creation of the response 

matrix for these constraints through the PL, PP, and FP procedures. The full linear 

program coefficient matrix includes the response matrix. Additional columns are 

added for the external and binary decision variables. Additional constraint rows are 

added for summation constraints and constraints controlling binary variables. This 

assembly of the linear program is accomplished in procedure FM. The problem is 

solved, assuming linearity, in the AP procedure. If the problem is nonlinear, then it is 

solved iteratively.

BAS Package—Basic Package for GWM: This package contains the AR procedure for 

reading the GWM names file and several utility subroutines for file handling and 

input/output. 
122 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

1Underscore indicates three-character designation of constraint type (HDC, STC, SMC, or DCC).

Adding New GWM Packages

GWM is organized to make the addition of new objective or constraint packages or 

solver packages relatively easy. Adding new decision-variable types, however, is difficult. 

The response-matrix method and the associated use of a linear-program formulation 

for the optimization problem inherently leads to mingling of decision variables and constraints within the solver. GWM is organized so that each constraint package operates independently of other constraint packages. This independence means that new constraint 

packages can be added without making any modifications to other constraint packages; 

however, all constraint packages are restricted to using the existing set of decision-variable 

types—the flow-rate, external, and binary decision variables. Adding new decision-variable types would require modification of every constraint package.

The same reasoning applies to the addition of new objective-function packages; 

namely, the objective-function package operates independently of the constraints but is 

closely tied to the three types of decision variables available in GWM.

Table 3-1. Outline of packages, procedures, and modules for GWM.

Package Procedures Module names

DCV—Decision Variables AR—read, allocate, and echo GWM1DCV1AR

FM—place flow rates in GWF equations GWF1DCV1FM

BD—add flow rates to GWF budget calculations GWF1DCV1BD

OT—write optimal solution of decision variables GWM1DCV1OT

OBJ—Objective Function AR—read, allocate, and echo GWM1OBJ1AR

FM—place into GWM equations GWM1OBJ1FM

OT—write final status of objective GWM1OBJ1OT

HDC, STC, SMC, 

DCC—Constraints

AR—read, allocate, and echo GWM1__1AR1

OS—observe state GWM1__1OS1

FM—place into GWM equations GWM1__1FM1

OT—write final status of constraints GWM1__1OT1

RMS—Response Matrix Solution 

Package

AR—read, allocate, and echo GWM1RMS1AR

PL—prepare perturbation loop GWM1RMS1PL 

PP—prepare flow-process run by setting 

perturbation of flow rate

GWM1RMS1PP

FP—formulate flow process by calculating 

response-matrix entries

GWM1RMS1FP

FM—formulate objective and all constraints 

through GWM1___1FM1 modules 

GWM1RMS1FM

AP—solve the optimization problem at an iteration 

and test convergence 

GWM1RMS1AP

OT—write optimization output GWM1RMS1OT

BAS—GWM Basic AR—read GWM name file GWM1BAS1AR

Various input/output utility routines
Appendix 3 123

Adding new solver packages is also possible. Presumably, any solver package will 

require repeated runs of the flow process followed by an assessment of results of the run, 

rules for the subsequent run, and assembly of information for ultimate solution of the optimization problem. By using the two loops depicted in figure 3-1 and appropriate content of 

the PL, PP, FP, FM, and AP procedures, many different optimization solver methods can 

be accomodated in the GWM structure.

Some Programming Details

Dimensioning of GWM Arrays: GWM uses FORTRAN-90 allocatable arrays for 

all internal storage. Arrays are defined in FORTRAN Modules at the beginning of each 

package file. Arrays and variables are accessed within each subroutine with the USE command. During execution, the size required for the arrays is determined from the relevant 

input files, and the arrays are dimensioned with the ALLOCATE command. In some cases, 

memory is allocated for arrays within a package and then later deallocated when the information is no longer needed. GWM reports to the output file the number of kilobytes of 

memory allocated for each of the packages.

Integration of GWM into MODFLOW MAIN program: GWM controls the 

repeated runs of the GWF Process through two loops depicted in figure 3-1—the flowprocess loop and the GWM iteration loop. Both loops use a FORTRAN DO-WHILE structure. In both cases, the number of executions of the loop that will be required cannot be 

determined a priori. For example, when the flow-process loop is used to calculate perturbations, the number of successfully completed flow-process runs is known, but if any run 

should fail (perhaps because of a lack of flow-process convergence), then the loop will have 

to be executed an additional time.

GWM makes use of three logical variables in the MODFLOW MAIN program to track 

the current flow-process run. These are FIRSTSIM, LASTSIM, and FINISH. FIRSTSIM 

is set to indicate that the flow process is being run for the first time, whereas LASTSIM 

indicates that the flow process is being run for the final time. If GWM is not activated, then 

both FIRSTSIM and LASTSIM are set to TRUE so that only one flow-process run is executed. If GWM is active, FIRSTSIM is used to initialize variables and perform other GWM 

setup functions, whereas LASTSIM is set once a solution to the optimization problem has 

been found. FINISH is set at the completion of the last flow-process loop to terminate the 

GWM iteration loop.

Repeated Execution of the GWF Process: The GWM Process requires repeated 

execution of the GWF Process. At each execution of the GWF Process, GWM performs a 

rewind of all GWF files by using module GWM1BAS1RW. For many situations, each subsequent run will use a set of managed flow rates that are very close to those used in a prior 

run. For steady-state problems, this means that the solution to the prior run is an excellent 

first estimate of the solution to the subsequent problem. To facilitate this, a new module, 

GWM1BAS1RPP, is called if GWM is active and it is not the first run. If the problem is 

not steady-state, the head information is read from the original files; the IBOUND array is 

read and the initial heads are read into array STRT and assigned to HNEW. If the problem 
124 GWM—A Ground-Water Management Process for the USGS Modular Ground-Water Model

is steady-state, however, neither STRT nor IBOUND are read. HNEW retains its value 

from the prior flow-process run and automatically becomes the initial head for the next run. 

GWM dependence on MODFLOW variables: GWM uses several MODFLOW 

variables to control progress of the solution algorithm. When modifying MODFLOW, consideration should be given to the definition of these variables. ICNVG is used to determine 

if the flow process has converged for a given flow-process run. If it has not converged, then 

GWM takes steps to perform the run again. HCLOSE is used as a measure of the precision 

of computed head values. It is assumed that HCLOSE has a value that reflects the head 

convergence criterion, regardless of the flow-process solver used.

Sign conventions on flow variables: The GWF Process treats withdrawals of water 

as negative values and injections of water as positive values. The GWM linear-program 

solver requires that all decision variables be lower bounded by zero, so that variables that 

represent withdrawals must be redefined as positive-value variables.

The following conventions are used to keep track of sign conventions. From the input 

perspective, all variables are treated as positive. This means that the bounds on withdrawals 

and injections are treated as positive values, and the variables are assumed to take nonnegative values when they appear in the objective function or summation constraints. 

Within the code, the flow-rate decision variables are stored by using the flow-process convention; for example, the array FVBASE stores the current base flow rates. The values in 

this array will be nonpositive for withdrawal variables and nonnegative for injection variables. Calculation of the response coefficients proceeds by adding to the FVBASE value a 

number stored in DELINC. This perturbation value may either be positive or negative. For 

the withdrawal variable, a negative value of DELINC implies a forward difference and a 

positive value implies a backward difference. The opposite is true for injection variables.

The difference in sign conventions is reconciled in subroutine SGWM1RMS1AP, 

where the sign on the response-matrix coefficients for withdrawal decision variables is 

switched immediately before the call to the linear-program solver. Because all other coefficients and bounds on the variables already have the positive sign convention, the 

problem is ready for solution. After return from the optimization solver, the signs of withdrawal decision variables are switched back to be consistent with the flow-process sign 

convention.
Cover: Hypothetical model grid illustrating some of the types of flow-rate and external 

decision variables that can be defined in the Ground-Water Management Process (GWM) for 

MODFLOW-2000. Flow-rate decision variables Q1, Q2, and Q3 represent ground-water 

withdrawal wells; flow-rate decision variable Q4 represents a recharge basin; and external 

decision variable IM1 represents a surface-water reservoir that is a source of water to the 

management model. (See figure 1-1 of this report and related discussions for more details.) 
For additional information write to:

Office of Ground Water

U.S. Geological Survey

411 National Center

Reston, VA 20192

(703) 648-5001