U.S. Department of the Interior

U.S. Geological Survey

Techniques and Methods Book 6, Chapter B6

One-Dimensional Transport with Equilibrium 

Chemistry (OTEQ)—A Reactive Transport 

Model for Streams and Rivers
Cover. Upper workings of the Pennsylvania Mine in the headwaters of Peru Creek, Colorado. Photograph by Robert L. Runkel,

U.S. Geological Survey, September 2009.
One-Dimensional Transport w

Equilibrium Chemistry (OTEQ

Transport Model for Streams 

By Robert L. Runkel

U.S. Department of the Interior

U.S. Geological Survey

Techniques and Methods Book 6, Chapter B

Toxic Substances Hydrology Program
U.S. Department of the I

KEN SALAZAR, Secreta

U.S. Geological Survey

Marcia K. McNutt, Dire

U.S. Geological Survey, Res

For product and ordering information:

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

Telephone: 1–888–ASK–USGS

For more information on the USGS—the 

its natural and living resources, natural h

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

Telephone: 1–888–ASK–USGS

Any use of trade, product, or firm names 

U.S. Government.

Although this report is in the public doma

reproduce any copyrighted materials con

Suggested citation:

Runkel, R.L., 2010, One-dimensional tran

model for streams and rivers: U.S. Geolo
Contents

Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

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

1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

1.2 Applicability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

1.3 Related Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

1.4 Report Organization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

1.5 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

2.2 Conceptual Model, Governing Equations, and the Seq

2.2.1 Model Assumptions . . . . . . . . . . . . . . . . . . . . . . . 

2.2.2 Derivation of Governing Equations . . . . . . . . . . 

2.2.3 A General Solution Scheme based on Sequen

2.3 Process Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

2.3.1 pH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

2.3.2 Precipitation/Dissolution. . . . . . . . . . . . . . . . . . . 

2.3.3 Sorption . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

2.3.4 Oxidation/Reduction. . . . . . . . . . . . . . . . . . . . . . . 

2.3.5 Transient Storage . . . . . . . . . . . . . . . . . . . . . . . . . 

2.3.6 Settling of Solid Phases. . . . . . . . . . . . . . . . . . . . 

2.4 Numerical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

2.4.1 The Conceptual System — Segmentation. . . . 

2.4.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . 

2.4.3 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 

3 User’s Guide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

3.1 Conceptual System, Revisited . . . . . . . . . . . . . . . . . . . . 

3.2 Input/Output Structure . . . . . . . . . . . . . . . . . . . . . . . . . . 

3.3 Input Format . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

3.3.1 Units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

3.3.2 Internal Comments . . . . . . . . . . . . . . . . . . . . . . . . 

3.3.3 The Control File . . . . . . . . . . . . . . . . . . . . . . . . . . . 

3.3.4 The Parameter File . . . . . . . . . . . . . . . . . . . . . . . 

3.3.5 The Flow File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

3.3.6 The MINTEQ Input File. . . . . . . . . . . . . . . . . . . . . 

3.3.7 The MINTEQ Database Files. . . . . . . . . . . . . . . . 

3.4 Input File Preparation and Model Execution. . . . . . . . 

3.4.1 Preparation of the Control File . . . . . . . . . . . . . . 

3.4.2 Preparation of the Parameter and Flow Files — 

3.4.3 Preparation of the MINTEQ input file — Use of 

3.4.4 Execution of OTEQ . . . . . . . . . . . . . . . . . . . . . . . . 
iv 

3.5 Output Analysis . . . . . 

3.5.1 The Solute and 

3.5.2 Concentration-D

3.5.3 The Post-Proce

3.5.4 Plotting Alternat

4 Model Applications . . . . . . . . 

4.1 Application 1: Time-Va

4.1.1 The Control File 

4.1.2 The Parameter 

4.1.3 The Steady Flow 

4.1.4 The MINTEQ Inp

4.1.5 Simulation Resu

4.1.6 Numerical Issue

4.2 Application 2: Time-Va

4.2.1 The Parameter 

4.2.2 The MINTEQ Inp

4.2.3 Simulation Resu

4.3 Application 3: Time-Va

4.3.1 The Parameter 

4.3.2 The Unsteady F

4.3.3 The MINTEQ Inp

4.3.4 Simulation Resu

4.3.5 Numerical Issue

4.4 Application 4: Steady-

4.4.1 Quasi-Steady-S

4.4.2 Modeling Existi

4.4.3 Simulation Resu

4.5 Application 5: Steady-

4.5.1 The Parameter, 

4.5.2 Specification of 

4.5.3 Simulation Resu

5 Software Guide . . . . . . . . . . . 

5.1 Supported Platforms . 

5.2 Software Distribution 

5.3 Installation . . . . . . . . . 

5.3.1 Creating the OT

5.3.2 Updating the Us

5.3.3 Creating User W

5.4 Compilation . . . . . . . . . 

5.5 Software Overview . . 

5.5.1 Model Developm

5.5.2 Include Files . . 

5.5.3 Error Checking 

References Cited . . . . . . . . . . . . . 

Glossary . . . . . . . . . . . . . . . . . . . . . 

Appendix 1. Modifications to the 
Figures

1 Conceptual surface-water system used to develop the gove

2 The sequential iteration approach for the reactive surface-w

3 Computation of the dissolution source/sink term. . . . . . . . . . 

4 Conceptual surface-water system for sorption to a static su

5 Use of the equilibrium submodel for a static surface. . . . . . . 

6 Conceptual surface-water system for sorption to a dynamic 

7 Conceptual surface-water system for sorption to static and 

8 Use of the equilibrium submodel for static and dynamic surf

9 Iterative scheme for oxidation/reduction. . . . . . . . . . . . . . . . . 

10 Conceptual surface-water system used to develop the gov

11 Segmentation scheme used to implement the numerical so

12 Upstream boundary condition defined in terms of a fixed co

13 Downstream boundary condition defined in terms of a fixed 

14 Conceptual system that includes one or more reaches. . . . 

15 The first reach in the conceptual system and the required i

16 OTEQ Input/Output files. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

17 Upstream boundary condition options.. . . . . . . . . . . . . . . . . . 

18 Example solute output file. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

19 Example precipitate output file, for the case of PRTOPT=1. 

20 Upstream boundary condition for double pulse injection. . 

21 Control file for Application 1. . . . . . . . . . . . . . . . . . . . . . . . . . . 

22 Parameter file for Application 1, record types 1–15. . . . . . . 

23 Parameter file for Application 1, record types 16–28. . . . . . 

24 Steady flow file for Application 1. . . . . . . . . . . . . . . . . . . . . . . 

25 MINTEQ input file for Application 1. . . . . . . . . . . . . . . . . . . . . 

26 Simulated concentrations of calcium (Ca) and sulfate (SO4

27 Peak total waterborne sulfate concentration at 200 meters

28 Partial listing of the parameter file for Application 2.. . . . . . 

29 Simulated and observed concentrations of (a) total dissolv

aluminum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 

30 Segmentation scheme for the hypothetical stream with trib

31 Partial listing of the parameter file for Application 3.. . . . . . 

32 Unsteady flow file for Application 3. . . . . . . . . . . . . . . . . . . . . 

33 Copper concentrations resulting from copper sulfate treatm

34 Spatial profiles of simulated copper concentration at 4 hou

35 Time required to reach quasi-steady-state conditions.. . . . 

36 Spatial profiles: pH, iron, and aluminum (existing condition

37 Partial listing of the parameter file for Application 5.. . . . . . 

38 Spatial profiles of simulated and observed: pH, iron, and ar

39 OTEQ directory structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 
vi 

Tables 1 The OTEQ control file. . . . . . . 

2 The parameter file — record 

3 The parameter file — record 

4 The parameter file — record 

5 The parameter file — record 

6 The parameter file — record 

7 The parameter file — record 

8 The parameter file — record 

9 The parameter file — record 

10 The parameter file — record 

11 The parameter file — record 

12 The parameter file — record 

13 The parameter file — record 

14 The parameter file — record 

15 The parameter file — record 

16 The parameter file — record 

17 Steady flow file — record ty

18 Steady flow file — record ty

19 Steady flow file — record ty

20 Unsteady flow file — record 

21 Unsteady flow file — record 

22 Unsteady flow file — record 

23 The MINTEQ input file — rec

24 The MINTEQ input file — rec

25 The MINTEQ input file — rec

26 Sorption parameters — hydr

27 Stand-alone MINTEQ runs to 

28 Supported systems. . . . . . . . 

29 Files to download.. . . . . . . . . 

30 Development environments.

31 Maximum dimensions and d

32 Aqueous species for which M

those in wateq4f.dat. . . . 

33 Aqueous species for which M

those in wateq4f.dat. . . . 

34 Mineral species for which M

those in wateq4f.dat. . . . 

35 Mineral species for which M

those in wateq4f.dat. . . . 
Conversion Factors

Temperature in degrees Celsius ( ̊C) may be converted to degre

 ̊F = (1.8 x  ̊C) + 32

Temperature in degrees Fahrenheit ( ̊F) may be converted to de

 ̊C = ( ̊F - 32) / 1.8

Specific conductance is given in microsiemens per centimeter 

25 ̊C).

Concentrations of chemical constituents in water are given eith

or micrograms per liter (μg/liter).

Multiply By To o

Length

centimeter (cm) 0.3937 inch 

millimeter (mm) 0.03937 inch 

meter (m) 3.281 foot 

kilometer (km) 0.6214 mile 

Area

square meter (meter2) 0.0002471 acre

square kilometer (kilometer2) 247.1 acre

square centimeter (centimeter2) 0.001076 squa

square meter (meter2) 10.76 squa

square centimeter (centimeter2) 0.1550 squa

Volume

liter 33.82 oun

liter 0.2642 gall

cubic meter (meter3) 264.2 gall

cubic centimeter (centimeter3) 0.06102 cubi

cubic meter (meter3) 35.31 cubi

cubic meter (meter3

) 0.0008107 acre

Flow rate

cubic meter per second (meter3/s) 70.07 acre

meter per second (meter/s) 3.281 foot 

cubic meter per second (meter3/s) 35.31 cubi

liter per second (liter/s) 15.85 gall

cubic meter per day (meter3/d) 264.2 gall

Mass

gram (g) 0.03527 oun

kilogram (kg) 2.205 pou

Energy

joule (J) 0.0000002 kilo

One-Dimensional Transport w

(OTEQ): A Reactive Transpor

Rivers

By Robert L. Runkel

Abstract

OTEQ is a mathematical simulation model used to charac

rivers. The model is formed by coupling a solute transport mod

model is based on OTIS, a model that considers the physical pr

age. The equilibrium submodel is based on MINTEQ, a model t

acid-base reactions, precipitation/dissolution, and sorption.

Within OTEQ, reactions in the water column may result in 

are subject to downstream transport and settling processes. Solid 

through dissolution and sorption/desorption reactions. Consider

phases requires a unique set of governing differential equations 

ferential equations describing physical transport and the algebr

sequential iteration approach. The model’s ability to simulate p

a means of evaluating the complex interactions between instrea

This report details the development and application of OT

specifications, model applications, and installation instructions

software/OTEQ. 
2 One-Dimensional Transport with Equilibrium Chemistry (OTEQ

1 Introduction

1.1 Overview

The study of solutes in streams and rivers is inherently comp

cesses influence solute fate and transport. Study of individual proce

transport processes that act to move solutes downstream and the b

Individual processes may be studied by employing simulation mo

work.

Studies of solute fate and transport commonly employ transpo

dispersion and some specific chemical and biological reactions (fo

and Hosseinipour, 1991; Chen and others, 1996). These models de

stants and empirical partition coefficients. This general approach is 

sparse. In addition, many sorption reactions are thought to adhere t

plexation). Although these transport models provide an accurate d

degree of chemical sophistication needed to describe pH-dependen

pH-dependent reactions in batch systems, but do not consider trans

so that local equilibrium may be reasonably assumed. It is therefore 

is used to describe physical processes and a chemical equilibrium m

is used here to develop OTEQ, a solute transport model that coupl

1.2 Applicability

OTEQ is generally applicable to solutes which undergo react

(the “Local Equilibrium Assumption”; Di Toro, 1976; Rubin, 198

and application dependent, many reactions involving inorganic solu

of chemical equilibrium, inorganic solutes may be modeled using 

facilitated through the use of an existing database that describes ch

tion, solute reactions not included in the existing database may be a

associated equilibrium constants. As such, OTEQ provides a gener

of chemical equilibrium. Despite this generality, most OTEQ appl

streams and small rivers. The remainder of this document is theref

note, however, that additional applications are possible.

1.3 Related Reading

Many of the algorithms used within OTEQ are based on the O

Inflow and Storage). The reader is therefore encouraged to review 

report. Copies of the OTIS documentation are available from the a

application of OTEQ also requires considerable knowledge of equi

treatment of chemical equilibrium problems presented by Morel a

Environmental Protection Agency’s MINTEQ program (Allison a

1.4 Report Organization

The remaining sections of this report are as follows. Section 2 

the reactive transport model. This section includes descriptions of 

and the numerical methods used within the model. Section 3, a Us

Fortran computer program. Model parameters, print options, and s

4 presents several applications of the model and includes example 

tion 5), describes how to obtain the model, installation procedures
1.5 Acknowledgments

The groundwork for OTEQ was completed in the late 198

McKnight, all of the U.S. Geological Survey (USGS). These re

metal fate and transport. Analysis of these experiments led to th

(USGS) conducted many of the early OTEQ applications and w

author thanks Ming-Kuo Lee (Auburn University) and David N

manuscript. Assistance with software development and distribut

provided by the USGS Toxic Substances Hydrology Program.
4 One-Dimensional Transport with Equilibrium Chemistry (OTEQ

2 Theory

This section describes the theoretical constructs underlying th

view of how physical transport and chemical equilibrium are coup

tial equations and a description of the general solution algorithm. T

ical processes are presented in Section 2.3. Section 2.4 concludes 

numerical methods, the conceptual stream system, and the treatme

2.1 Overview

The reactive transport model is formed by coupling the OTIS 

chemical equilibrium submodel. The resultant model considers a v

dispersion, transient storage, the transport and deposition of water

tation/dissolution, and sorption. Consideration of these processes p

solute fate and transport.

Solute Transport Model. The OTIS solute transport model i

with additional terms to account for lateral inflow and transient sto

noted in many streams, where solutes are temporarily detained in e

the faster moving water near the center of the channel. In addition

within the streambed), where solutes are also detained. Lateral inf

surface inflow, overland flow, interflow, and ground-water discha

equations describing the physical transport of multiple solutes.

Equilibrium Submodel. The chemical equilibrium submode

of the MINEQL model developed by Westall and others (1976). G

MINTEQ computes the distribution of chemical species that exist 

putations include the precipitation and dissolution of solid phases a

equations describing equilibria form a set of nonlinear algebraic e

The conceptual and mathematical framework underlying MIN

ers, 1976; Morel and Hering, 1993; Allison and others, 1991). As 

are given here. Chemical “components” are defined as the fundam

derived. Chemical reactions involve two or more components that 

are selected such that (1) the components combine linearly to form 

as a combination of other components (Westall and others, 1976).

bining chemical components. The chemical equilibrium problem e

librium. This is accomplished by developing mass-action equation

equations for the chemical components.

Coupling Transport and Equilibrium Chemistry. Couplin

set of algebraic and partial differential equations. The sequential i

set of equations by dividing each time step into a “reaction” step a

submodel is executed for each segment in the stream network. Eac

rium is assumed. The equilibrium submodel thus determines the s

on this information, a transport step is taken in which the solute tr

solute. Because the transport and reaction steps neglect the couplin

specified level of convergence is achieved.

1

One exception to this general rule is the case of multiple oxidations states. For 

Fe(III), and e– are all components.
2.2 Conceptual Model, Governing Equations, a

2.2.1 Model Assumptions

The governing equations and solution algorithms used with

tions:

• Chemical Equilibria. Complexation, precipitation/dis

Under this “Local Equilibrium Assumption,” chemica

processes (Di Toro, 1976; Rubin, 1983). This assumptio

One exception to the equilibrium approach is the kineti

streambed (Section 2.3.3).

• One-Dimensional Transport. Solute mass is uniform

dimensional transport is applicable. Given this assump

consists of a series of stream segments (control volume

segment include advection, dispersion, lateral inflow, t

sorbed species resident in the water column travel at th

• Transient Storage. The physical process of transient s

(Runkel, 1998): advection and dispersion are not includ

negligible; the exchange of solute mass between the m

transfer process (Bencala and Walters, 1983). Chemica

• Physical Parameters. All model parameters describin

describing advection and lateral inflow may be tempor

main channel cross-sectional area, lateral inflow rate, an

model parameters are temporally constant.

• Chemical Parameters. All model parameters associate

stant.

• Mobile and Immobile Phases. Solute mass for each c

first three phases represent dissolved, precipitated, and 

are mobile, in that they are subject to transport. The fin

on immobile substrate (the streambed or stationary deb

layer of solute mass that interacts with the overlying w

• Precipitation. Dissolved mass in the water column ma

respect to the defined solid phases. Any precipitated ma

until it settles to the streambed or redissolution occurs. 

tation directly to the immobile bed is excluded. Precipi

precipitates are subject to the force of gravity and settle 

is unaffected by other solids in solution (Section 2.3.2)

• Dissolution. When the aqueous solution is undersatura

of the precipitate in the water column is allowed to dis

This assumption is based on the intimate contact betwee

2.3.2).

• Sorption. Dissolved species may sorb to solid phases i

versely, sorbed species may desorb from sites in the w

to sorption are presented in Section 2.3.3.

• Oxidation/Reduction. Solute mass may be transferred 

reduction reactions (Section 2.3.4).
6 One-Dimensional Transport with Equilibrium Chemistry (OTEQ

2.2.2 Derivation of Governing Equations

With these assumptions, the fundamental equations governing 

formulated in terms of the chemical components defined within th

T, is the sum of the dissolved (C), mobile precipitate (Pw), immob

(Sb) phases:

T CP 

w P

= + +

where each phase consists of one or more chemical species. For exa

Cc 

i = 

M

+= ∑

where

c concentration of the uncomplexed component species (Fe

xi concentration of the ith complexed species [moles per lite

ai stoichiometric coefficient of the component in the ith com

M number of complexed species;

and species concentrations (c and xi) are provided by equilibrium c

Pw+Pb) and total sorbed concentrations (S=Sw+Sb) are given by Y

A summary of the processes considered for each phase is pre

partments. The water column compartment contains the three mob

debris) constitutes the second compartment, containing the two imm

tified using source/sink terms (fb, fw, gb, gw; see arrows, fig. 1). Th

resented by the transport operator, L( ). The dissolved phase, C, tak

tions that occur within the water column (interactions with Pw and S

by dissolution of precipitate from the immobile substrate and by so

and Sb; fb and gb arrows, fig. 1). Finally, C may increase or decrea

exchange between the atmosphere and the water column, for exam

settle in accordance with the settling velocity, ν1 [LT–1].

Figure 1. Conceptual surface-water system used to develop the gov

consists of dissolved (C), mobile precipitate (Pw), immobile precipitate 

dissolved and mobile phases are subject to transport, as denoted by L

terms (fb, fw, gb, gw, sext) and a settling velocity (v1) as described in th

2

A full listing of chemical components is provided in the MINTEQ documentati

anions (chloride, sulfate, fluoride), cations (aluminum, ferrous iron, ferric iron), co

Water Colu

Immobile S

C

Pb 

Pw 

v1

L( ) 

fw 

fb 

Conceptual Surfac
A general mass-balance equation for each component is d

phases within a stream segment (control volume). An equation 

oped by summing the equations for the individual phases. In the 

not treated as separate control volumes, but rather as a single co

and others, 1960). Note that this approach differs from the appr

stances. These models are often developed for rivers and lakes 

column. In this instance, two or more control volumes are used t

we are concerned with streams where only a thin, immobile laye

column. As such, treatment of the system as a single control vo

Mass-balance equations for the five phases are developed 

component consists of a single species. The dissolved and sorbe

be composed of multiple species. The problem of multiple prec

Mass balances for the five phases are given by:

Dissolv

∂C

∂t

------ L C( ) f w = + +

Mobile 

Pw ∂

∂t --------- L Pw ( ) 

= 

Mobile 

Sw ∂

∂t -------- L Sw ( ) 

= 

Immobile 

dt

dPb v1

d

---- 

=

Immobi

dt

dSb v1

d

---- 

=

where3

L transport operator;

fw source/sink term for precipitation/dissolution from the 

fb source/sink term for dissolution from the immobile sub

gw source/sink term for sorption/desorption from the wate

gb source/sink term for sorption/desorption from the imm

sext source/sink term representing external gains and losses 

v1 main channel settling velocity [LT–1];

d1 effective settling depth [L] (see Section 2.3.6); and

t time [T].

Given these mass-balance equations, several comments ar

implicit functions that are dependent on the solution of the non

sink terms are not explicitly provided by the equilibrium submod

Second, the external source/sink term (sext) represents mass that 

(or sink) that is external to the system; unlike the other source/s

phases. For example, the equilibrium submodel may be used to d

CO2. This use of the equilibrium submodel results in a gain (tra

sing) of mass due to an external source/sink. Another example 

Additional details on sext are given in Section 2.2.3 and by Runk

transient storage model (Bencala and Walters, 1983; Runkel, 1

3

The fundamental units of Length [L] and Time [T] are used throughout this 
8 One-Dimensional Transport with Equilibrium Chemistry (OTEQ

L C ̃ ( ) Q

A

--- C

 ̃ ∂

∂x – ------ 1A

---∂x

∂ AD C

 ̃ ∂

∂x

( ) ------ = + + 

where

A main channel cross-sectional area [L2];

main channel concentration of an arbitrary phase [moles p

lateral inflow concentration of the arbitrary phase [moles 

storage zone concentration of an arbitrary phase [moles p

D dispersion coefficient [L2T–1];

Q volumetric flow rate [L3T–1];

qLIN lateral inflow rate [L3T–1L–1];

x distance [L]; and

α storage zone exchange coefficient [T–1].

Use of the transient storage approach introduces an additiona

tions, . The storage zone equations are discussed in Section 2.3

parameters that apply to the main channel and those that apply to t

to denote the main channel; the subscript “2” denotes the correspo

The mass-balance equation for the total component concentra

the five individual phases. This yields:

2.2.3 A General Solution Scheme based on Sequ

The basic problem is as follows. Given the mass-balance equ

model be used to determine the component concentrations in the v

and Tripathi (1989), the total component concentration consists of 

concentrations in each of these three phases are readily available a

use of the equilibrium submodel is confounded by the presence of 

itated and sorbed mass into mobile and immobile fractions. Fortuna

tiple sorptive surfaces, so that separate surfaces may be defined for 

between the mobile and immobile sorbate concentrations. For the 

and only a single value reflecting the total amount of precipitate is 

and immobile precipitate is therefore required.

The solution technique presented here uses total component c

for T is presented as equation 9. This equation is analogous to the 

1989). Here an implicit form is developed by combining equation

∂T

∂t

------ L T( ) L Sb – ( = 

where equation 10 is formulated such that the waterborne phases a

tion of Pb and Sb. Equations for these phases are given by:

dt

dPb v1

d

---- P P

( – = 

dt

dSb v1

d1

---- S S

( – = 

where P and S are the total precipitated (=Pw + Pb) and total sorbe

from eqs. 6 and 7).

The equation set governing the problem consists of three parti

and the set of algebraic equations representing chemical equilibria

mation of the governing differential equations and the sequential ite

additional nomenclature. Let n denote an initial time and n+1 den

state of the system is known, and time n+1 is the current time for 

denote the iteration number. Finally, a caret is used to indicate tha

estimate of the dissolved concentration for the current iteration at th

tion scheme and the Crank-Nicolson method are provided in Secti

C

 ̃

C

 ̃ L

C

 ̃ S

C

 ̃ S

∂T

∂t

------ LC Pw + + ( = 

1
The goal of sequential iteration is to solve the set of partial 

equation in the form of equation 10 for each chemical compone

levels are needed to solve for the total component concentratio

state variables at time level n are available from the previous ti

time level n+1. Specifically, estimates of P, Pb, S, and Sb are n

below, P and S are provided directly by the equilibrium submod

sink terms are developed algorithmically. Iteration is required b

estimates of the variables at the advanced time level. Solution of 

equilibrium calculations, transport calculations, and convergen

Tˆ n + 1

Begin

Time Step

n = n+1

k = 0

Step 1:

Estimate 

Concentrations

N

at time n+1

 Pˆ n + 1, k 

Figure 2. The sequential iteration approach for the reactive surf

Sequ
10 One-Dimensional Transport with Equilibrium Chemistry (OTE

Step 1: Initialization 

As each time step begins, the total component concentration 

( ). These concentrations are used as input to Step 2. Values f

by simply setting equal to T n. This estimation procedure is o

pleting the equilibrium and transport steps within the first iteration

described below.

Step 2: Equilibrium Calculations

Step 2 begins the iterative loop (fig. 2). For each chemical co

( ) is checked to ensure that it is a valid concentration. If 

liter), is reset to the prescribed minimum. This procedure en

passed to the equilibrium submodel. Final estimates are then input 

chemical species are computed. The concentrations of the individu

dissolved, precipitated, and sorbed phases (C, P, and S). These qua

the sorption/desorption source/sink (gb) as described in Sections 2

Gain or loss via external sources and sinks is quantified by co

the equilibrium calculations. As before, the total component conce

as . The total component concentrations after equilibration 

phases. The external source/sink term for each component is there

sˆ 

ext

n + 1, k Cˆ n + 1, k Pˆ n + 1, k 

+ 

= -------------------------------------------

where Δt is the integration time step [T].

Step 3: Transport Calculations

Given estimates of fb, gb, and sext at time n+1, the equations 

and 12 are first solved for the immobile precipitate and sorbed phas

with the variables from time n to solve equation 10 using the Cran

sents one of two states. If the solution has converged (Step 4), this 

port problem for the current time step. If convergence is not obtain

trations used as input for Step 2 in the next iteration. 

Step 4: Convergence Test

In Step 2, phase concentrations at the advanced time level (n+

calculations are based on estimates of the total component concentr

these estimates are based on the previous value of T, as described u

on the solution of the transport equations (Step 3) from the previous 

should approach T n+1

, and the solution converges. The algorithm t

convergence is performed. This convergence test is given by:

Tˆ n + 1, k + 1 Tˆ n + 

–

Tˆ n + 1, k + 1 --------------------------------------

where σe is a relative error tolerance. If equation 14 holds for all co

time step is initiated. If the left-hand side of equation 14 is greater 

converged and another iteration is required.

Tˆ n + 1 

Tˆ n + 1

Tˆ n + 1 Tˆ n + 

Tˆ n + 1

Tˆ n + 1, k
2.3 Process Formulation

A general algorithm for the solution of the reactive transpo

sents a detailed description of the solution techniques used to im

include pH, precipitation/dissolution, sorption, oxidation/reduc

2.3.1 pH

Models of chemical equilibria generally use one of two ap

approach, a charge balance equation is used to determine the aq

the PHREEQC equilibrium model (Parkhurst and Appelo, 1999

Morgan, 1972), is used within the reactive transport model. Un

“excess” hydrogen ions in solution. The proton condition appro

aqueous component that is subject to transport (eq. 10); no spec

1991). Within the equilibrium submodel, the proton condition i

T 

H H + ∑ specie

= 

where TH is the total component concentration for excess hydro

of pH are provided in Sections 4.2–4.5.

2.3.2 Precipitation/Dissolution

Step 3 of the sequential iteration procedure (Section 2.2.3

immobile precipitate equation (eq. 11). As shown here, informa

mate fb for each stream segment. To begin, consider the change i

the differential equations derived for each phase, the change in 

Estimation of fb from equation 16 is based on several simplifyin

precipitate may increase (∂P/∂t > 0). This increase may be due 

[L(Pw) > 0; eq. 8]. If precipitation is occurring, dissolution is no

increase if the gain due to transport is greater than the loss due to 

fb is also zero, as the gain due to transport indicates the presenc

from the water column, such that a nonzero Pw implies an fb of 

decreases (∂P/∂t < 0). This decrease may be due to dissolution 

At a given location, dissolution from the bed occurs (fb > 0) on

This observation is used to eliminate L(Pw) from equation 16. H

final assumption is not entirely valid, as situations may arise in 

lem, however, as the presence of precipitate in the water colum

An algorithm for determining fb is presented in figure 3. U

rium submodel is called to compute the total amount of precipit

exceeds that for time n ( > Pn), the net amount of precip

the case of precipitation, fb equals zero. If, however, the net am

occurred and two situations are possible. First, if net dissolutio

present in the water column (Pw

n), all of the dissolved mass is ta

Here again fb equals zero. Second, if net dissolution is not accou

from the immobile phase. In this case, fb is given by:

f 

ˆ b

n + 1, k Net Di

= -------------------

∂t

∂P L Pw ( ) 

– = 

Pˆ n + 1, k
12 One-Dimensional Transport with Equilibrium Chemistry (OTE

The Problem of Multiple Precipitates. The equations devel

phase consists of a single species. In this section the introduction o

ined. Note that equation 6 is developed by considering the change i

When more than one precipitate is present for a given component, 

precipitate species may settle at different velocities (hence v1 cann

immobile phase is species-specific. To consider multiple precipitat

cipitated species:

dt

dpbm v1m

d1

-------- pm 

( – = 

where

pbm immobile precipitate concentration for precipitated specie

pm total precipitate concentration for precipitated species m [

v1m settling velocity for precipitated species m [LT−1]; and

fbm source/sink term for dissolution of immobile precipitated 

Pˆ 

n + 1, k P

>

Begin Ste

input 

Tˆ n + 1, k Equilibriu

Submod

Yes

f

ˆ b

n + 1, k = 0 

Begin 

(Precipitation) 

f

ˆ b

n + 1, k 

= 

Figure 3. Computation of the dissolution source/sink term (fb).

Dissolutio
Solution of the reactive solute transport problem now requ

dure (Section 2.2.3) is modified as follows. First, rather than so

cies associated with a given component. To solve equation 18, 

the equilibrium submodel. The source/sink term for species m ( 

ogous to that for . After solving the m equations to obta

concentration of immobile precipitate is given by:

P 

ˆ b

n + 1, k 

m = 

n

= ∑

where am is the stoichiometric coefficient of the component in th

species for the current component. The solution of the governin

2.3.3 Sorption

Step 3 of the sequential iteration procedure (Section 2.2.3

solve the immobile sorbate equation (eq. 12). As shown here, i

to estimate gb. Mathematical descriptions of sorption range from 

resentations based on electro-chemical theory. Because distribu

is limited in metal-contaminated waters where the primary sorbe

tive transport model therefore estimates gb using the generalized 

ation model that explicitly considers the effects of pH and ionic 

Generalized Two-Layer Model (GTLM). This section d

1990) as implemented within the equilibrium submodel. As wi

model defines sorption reactions in terms of mass law equation

sites at equilibrium. The equilibrium constant associated with a 

resenting the chemical free energy of site binding and a second 

electrostatically charged surface. The coulombic term acts as a 

move ions from the surface layer to the bulk solution.

Because the coulombic term varies as a function of surfac

ranged and expressed in terms of intrinsic surface complexation 

SOH M 

2 + + 

↔ 

where M2+ is a divalent cation, H+ is a hydrogen ion, SOH is an 

charged surface species. The corresponding mass-action equati

K 

SOM { 

{ SOH 

= --------------

where K is the equilibrium constant and {} denotes chemical ac

terms and rearranging yields:

K 

int SO

{ 

{ } SOH M{ 

= ------------------------

where Kint is the intrinsic surface complexation constant, exp(−

[volts], F is the Faraday constant [96,485 coulomb mole−1], R i

lute temperature [K]. 

Solution of a chemical equilibrium problem that includes 

component to account for the coulombic correction factor and a 

ory, surface charge is balanced by a diffuse layer of counter ch

face potential is defined by Guoy-Chapman theory (Dzombak a

σ 8RT 

aεε0ce

= 

f

ˆ b

n + 1, k
14 One-Dimensional Transport with Equilibrium Chemistry (OTE

where σ is net surface charge density [coulomb meter−2

], ε is the d

[8.876×10−12 coulomb volt−1

 meter−1], ce is the molar electrolyte c

Given equation 23, the total component concentration for the coul

T 

CC σ

S

= ---

where CC = exp(−ΨF/RTa), SA is specific surface area [meter2 per 

liter].

Solution of the chemical equilibrium problem also requires s

A central part of the generalized two-layer model is the postulation 

ing. The first type, the high-affinity site, is generally less prevalen

A second low-affinity site is in greater abundance but has weaker bi

components are introduced for each sorptive surface. Total compo

given by:

T 

SOH

N

= ----

where NS is the site density [moles of sites per mole sorbent] and M 

sorbent].

GTLM within the Reactive Transport Model. The reactive 

static and dynamic sorptive surfaces. Static sorptive surfaces are th

change in time. Conversely, dynamic sorptive surfaces are those fo

An example of a static sorptive surface is a streambed armored wit

the number of sites available for sorption reactions is relatively co

may arise in which the number of sites changes in time, and sorpti

hydrous metal oxides form in the water column as a result of prec

defined as dynamic sorptive surfaces.

To further classify the sorption reactions, it is useful to subdi

the static sorptive surface; Pool 2 consists of the dynamic surfaces 

borne precipitates; and Pool 3 consists of dynamic surfaces that w

during the course of the simulation. The sorbed concentrations ass

respectively. Additional assumptions underlying the use of GTLM 

• Sorption reactions adhere to the generalized two-layer mo

• A static surface and(or) a dynamic surface may be defined. 

surface is distributed between Pools 2 and 3 as defined ab

• Specific surface area (SA) and sorbent molecular weight (M

for each site type on each sorptive surface. Sorbent proper

Given these assumptions, three cases are possible: (1) sorptio

sorption to static and dynamic surfaces. As shown below, these ca

equilibrium submodel is used, and how kinetic limitations are imp

Sorption to a Static Surface

A conceptual diagram depicting sorption to a static surface is 

from the dissolved phase, C, to the sorbed phase associated with P

alent to the immobile sorbed phase, Sb. During desorption, mass i

governed by the kinetic parameter, Γ. Sorption to a static surface r

• The solid concentration (SC) does not change in time. The 

or debris) and is therefore not subject to downstream trans

• The solid concentration is allowed to vary spatially on a r
Figure 4. Conceptual surface-water system for sorption to a stati

The primary task is to solve equation 12 for the concentra

surfaces, the settling term in equation 12 drops out, yielding:

dt

dSb

=

where gb is estimated using output from the equilibrium submo

equilibrium sorption and kinetically limited sorption.

Equilibrium Sorption. The first step in modeling equilib

(T n+1) used as input to the chemical equilibrium submodel. Fo

equation 10 from the previous sequential iteration. Total compo

high- and low-affinity sites on the static sorptive surface are giv

equation 24 is equal to the temporally constant, spatially variabl

of SC, SA, NS, and M are shown in table 26 (Section 3.4.3), Sec

As shown in figure 5a, the equilibrium submodel determin

sorbed phases (C n+1, P n+1, S n+1). The sorption source/sink ter

during the current time step: 

g 

b

Sn 

-----=

Given this definition of gb, equation 26 may now be solved usin

Sw equals zero when only a static surface is considered ( 

As shown by equation 28, the submodel provides the exact qua

equality given by equation 28, rather than a formal solution of 

Water Co

Immobil

Pb

Pw

v1

L( ) 

Sorption to a 

Sn =

Sb

n + 1 

=
16 One-Dimensional Transport with Equilibrium Chemistry (OTE

Figure 5. Use of the equilibrium submodel for a static surface: (a) eq

Kinetically Limited Sorption. Under the equilibrium approa

rium submodel. Here we employ a pseudo-kinetic approach in whi

sidered. This kinetic limitation is designed to model cases where o

with the static surface (the streambed). Figure 5b depicts use of th

putations during Pass 1 are very similar to the computations descr

mines the total component mass in the dissolved, precipitated, and 

equilibrium concentration in the absence of a kinetic limitation). T

change in sorbed concentration:

g 

b Γ Sn S

–

Δ

-------------

⎝ 

⎛ 

= 

where Γ is the fraction of the equilibrium quantity that is allowed t

using a forward time difference yields:

Sb

n + 1 Sbn Γ S

+= ( 

During Pass 1, output from the equilibrium submodel reflects 

rium. The phase concentrations from the submodel (C, S, P) and th

limitation. To incorporate these effects, the total sorbed concentrat

), and Pass 2 is initiated. During Pass 2, sorption reactions ar

revised to eliminate sorbed mass:

Given T ′, the equilibrium submodel provides the corrected values 

T′ Tn + 1 Sn + 1 –=

Tn + 1 

Tn + 1 

a. Equilibrium Sorption

b. Kinetically Limited Sorption

Pass 1

Pass 2

Equilibrium Submode

Sb

n + 1

T′ Tn + 1 –=
Sorption to a Dynamic Surface

Reactions between the dissolved phase, C, and the dynami

may sorb to the waterborne precipitates in Pool 2, thereby incre

equivalent to the mobile sorbed phase (Sw) and is therefore subj

sorbed mass may settle, increasing S3. Alternatively, desorption 

Sorption/desorption reactions also transfer mass between the di

ered, S3 is equivalent to the immobile sorbed phase, Sb. Assum

• In Pool 2, the concentration of the sorptive solid (SC) v

centration of a specified hydrous metal oxide (Pw

Me). T

column and are subject to transport and settling.

• Pools 2 and 3 are modeled using a single surface define

identical specific surface areas (SA), molecular weights 

S3 based on each pool’s contribution to SC.

• In Pool 3, the concentration of the sorptive solid (SC) v

concentration of a specified hydrous metal oxide (Pb

Me

attached to immobile substrate (the streambed or debri

• Sorption/desorption reactions in Pools 2 and 3 are in lo

Figure 6. Conceptual surface-water system for sorption to a dyna

Given these assumptions, the task is to solve equation 12 f

face in Pool 3, Sb. As with a static surface, a formal solution to 

be obtained directly from the equilibrium submodel. The settlin

of sorbed mass is reflected in the settling term for the hydrous m

that makes up the dynamic surface, as given by the governing e

Use of the equilibrium submodel for the dynamic surface 

conditions (fig. 5a). The only difference is computation of the to

(eq. 24) and the surface components (high- and low-affinity site

the concentration of the precipitate defined as the dynamic surf

Water Co

Immobil

Pb

Pw

v1

L( ) 

Sorption to a Dy
18 One-Dimensional Transport with Equilibrium Chemistry (OTE

S 

C MP

=

where PMe is the total precipitate concentration for the specified hy

conditions, the submodel provides values for Cn+1, Pn+1, and Sn+1. I

surface. The amount of Sn+1 on the streambed is determined by con

streambed:

S 

b

n + 1 Pb

Me

PMe ---------

=

Sorption to Static and Dynamic Surfaces

Concurrent simulation of sorption to static and dynamic surfa

solved phase and the various sorbent pools are shown. As before, s

ically limited, whereas sorption/desorption reactions for the dynam

Figure 7. Conceptual surface-water system for sorption to static and 

Equilibrium Sorption. Use of the submodel for equilibrium 

concentrations for the coulombic correction factor and the surface 

above (eqs. 24 and 25). The submodel provides values of Cn+1, Pn

Pools 2 and 3). The total sorbed concentration is the sum of the co

S 

n + 1 S1

n + 1 

=

whereas the immobile sorbed concentration is the sum of the conc

Water Column

Immobile 

C

Pb

Pw

v1

L( )

Po

Γ

S

Sorption to Static and 

Sb

n + 1 S1n + 1 

+=
Kinetically Limited Sorption. Use of the submodel for k

values of Cn+1*, Pn+1*, S1

n+1*, and S2,3n+1*. The kinetically lim

As with a static surface, phase concentrations from Pass 1 and s

C, P, S2,3, and pH are determined during Pass 2 (fig. 8b), where

In Pass 2, sorption to the static surface is not considered, as the 

equation 36. Total sorbed and immobile sorbed concentrations 

Figure 8. Use of the equilibrium submodel for static and dynamic 

S1

n + 1 S1n Γ 

+= 

T′ Tn + =

T′ Tn + 1 S1

n + 1 –=

Tn + 1 

Tn + 1 

a. Equilibrium Sorption

b. Kinetically Limited Sorption

Pass 1

Pass 2

Equilibrium Submodel: Stati
20 One-Dimensional Transport with Equilibrium Chemistry (OTE

2.3.4 Oxidation/Reduction

Oxidation/reduction reactions play an important role in determ

ric iron concentrations, for example, are often controlled by photor

Kimball and others, 1992). Previous investigators have used a fixe

solved iron (Broshears and others, 1996; Runkel and others, 1996b

(1996) in a generic manner to fix the dissolved fraction of two cou

specific oxidation/reduction reactions are not modeled. Despite thi

mechanism whereby the proper ratio of two coupled components i

reduction (θtarget; eqs. 40 and 41, below) is based on commonly av

rous and ferric iron).

To begin, consider two aqueous components that exchange so

total dissolved concentration, Ctot, is the simple sum of all dissolv

C 

tot C1 

+=

where superscripts “1” and “2” denote the first and second compo

considered by specifying a single parameter, θtarget, to dictate the 

with the first component.

A schematic diagram of the fixed fraction approach is given 

placed around the sequential iteration procedure described in Sect

iterative procedure for achieving θtarget. As each iteration begins, th

concentrations as a function of transport and chemistry. After conv

concentration is computed, as well as the fraction of the total disso

each stream segment. This fraction is then compared with that from 

The convergence test is given by:

θ 

m + 1 θm – < 

where m denotes the iteration. If equation 39 holds for all segment

If the solution has not converged, total component concentrations 

T 

1 θtargetCtot 

= 

T 

2 1 – θtarget ( )C

= 

where P and S denote total precipitate and total sorbed concentrati

the two components are computed by comparing the quantities give

[for example, =(T1,m+1 − T1,m s )/Δt]. A new iteration of the oxi

ext

1
B

Tim

n = 

m 

1. Apply θt

2. Compute

 

Figure 9. Iterative scheme for oxidation/reduction.

Oxidatio

source/si

external
22 One-Dimensional Transport with Equilibrium Chemistry (OTE

2.3.5 Transient Storage

As given by equation 8, the transport operator describes the ph

lateral inflow. By including transient storage, equation 8 introduce

Due to the introduction of the storage zone concentration, addition

transport problem that includes transient storage is illustrated sche

main channel, while the bottom part represents the storage zone. Se

are shown in the figure. First, mass in the water column of the mai

port operator, L. Mass in the storage zone, meanwhile, is not affect

described for the main channel also take place in the storage zone

with the storage zone through first-order mass transfer. This is rep

storage zone counterparts, Cs, Psw, and Ssw.

This last point is an important assumption for the model pres

approach, the storage zone represents both open water (eddies and 

(flow through porous areas within the streambed). This lumped sys

as the solutes of interest are often dissolved-phase tracers that are 

case of reactive solute transport, the modeled solutes are composed 

is primarily open water, exchange of both dissolved and solid phas

assumption. If the storage zone consists of water in the hyporheic 

questionable assumption, as solid mass entering the storage zone m

between the solid and the porous media that comprises the hyporh

The lumped nature of the storage zone therefore presents a pr

storage zone consists of open water, the correct model formulation 

solved and solid phases. Conversely, if the storage zone is primari

affect the dissolved phase only. In reality, storage zones are usuall

the development that follows, the open-water scenario is assumed; 

is a logical choice: it correctly models the open water portions of t

the solid particles are small relative to the pore space, such that so

To implement the transient storage approach, mass-balance eq

The total component concentration in the storage zone, Ts, is given 

T 

s Cs P

= + 

P 

s Psw += 

S 

s Ssw 

+= 

where Cs, Ps, and Ss are the total dissolved, precipitated, and sorbe

mobile (waterborne) precipitate, immobile (bed) precipitate, mobil

with the main channel, conservation of mass applies to the five ph

five mass-balance equations is similar to that given for the main cha

zone (Cs, Psw, Ssw) are not affected by advection, dispersion, or late

equation for Ts:

dt

dTs α

A

As

----- T Pb – Sb – T

– ( = 

where AS is the storage zone cross-sectional area [L2] and ssext is a 

storage zone. The main channel concentrations T, Pb, and Sb are gi

trations are defined by:

dt

dPsb v2

d2

---- Ps P

( – = 

dt

dSsb v2

d2

---- Ss S

( – = 
where

fsb source/sink for dissolution from the storage zone immo

gsb source/sink for sorption/desorption from the storage zo

v2 storage zone settling velocity [LT−1]; and

d2 effective storage zone settling depth [L].

Wa

Pb 

Pw 

v1

Psb 

Psw 

v2

St

L( )

Immo

Imm

Figure 10. Conceptual surface-water system, including transien

total component concentration is as defined in figure 1. The total 

mobile precipitate (Psw), immobile precipitate (Psb), mobile sorbe

Conceptual Surface-Water 
24 One-Dimensional Transport with Equilibrium Chemistry (OTE

The transient storage equations are easily incorporated into th

tialization phase (Step 1) involves the additional task of estimating 

These concentrations are used as input to the equilibrium submode

centrations for the storage zone (Step 2). The results from the equili

for the storage zone in a manner analogous to that used for the ma

the main channel and the storage zone in each stream segment. In 

and the storage zone are solved for the dependent variables (eqs. 1

10 is then solved for the total component concentration in the mai

pling procedure described by Runkel and Chapra (1993, 1994). Fin

concentration, Ts.

2.3.6 Settling of Solid Phases

Equations 11 and 12 model the settling of precipitated and so

depth, d1 is determined from the depth-area power function:

d 

1 aA

=

where a is the coefficient of the power function and b is the expon

angular channel, b=1 and a=1/w, where w is channel width; and (2) 

b=0 (d1=aA0=a).

2.4 Numerical Solution

The governing equations described above (eqs. 10–12 and 45

must be approximated using numerical methods. The finite-differen

consistent with those presented for the OTIS solute transport mod

nomenclature is presented herein, however, to further define the c

2.4.1 The Conceptual System — Segmentation

To implement a numerical solution scheme, the physical syst

in which the stream is subdivided into N discrete segments. Each of 

is conserved. Equations 10–12 and 45–47 therefore apply to each se

the subscripts i, i-1, and i+1 denote concentrations and parameters 

(i-1,i) and (i,i+1) define values at the segment interfaces. The leng

Figure 11. Segmentation scheme used to implement the numerical s

crete segments (control volumes).

Segment ...

Segmentation 

Q

Interfa

i-1, i

Segment

1

Segment

2

Segment

i-1

x
2.4.2 Boundary Conditions

Equation 10 is solved using finite-difference approximatio

ute transport model (Runkel, 1998), a central differencing schem

i is dependent on the concentrations in the neighboring segment

and last segments of the modeled system (segment i-1 is undefi

N; fig. 11).

Upstream Boundary Condition (Segment 1). The upstre

at the upstream boundary ( , fig. 12). This boundary conditi

upstream end of the modeled system. As such, observed concen

tion, as discussed in Section 3.3.4 and Section 4.

Figure 12. Upstream boundary condition defined in terms of a fixe

Downstream Boundary Condition (Segment N). Direct 

eled system results in the introduction of a concentration term f

Because segment N+1 is outside the modeled system, concentra

task is accomplished by defining a dispersive flux at the downs

D 

C

 ̃ ∂

∂x

------ ⎝ ⎠ ⎛ ⎞

N N, + 1

= 

where DSBOUND is a user-supplied value for the dispersive flu

the concentration gradient (change in concentration with respec

Because non-zero concentration gradients are likely to exist in 

should be placed well downstream of the spatial locations of in

DSBOUND is minimized. User specification of DSBOUND is d

and the potential error is discussed in Section 4.1.

Figure 13. Downstream boundary condition defined in terms of a 

C

 ̃ bc

Segment 1 — Upstre

flow

C

 ̃

bc 

Segment N — Downstre

flow

C 

 ̃ N – 1
26 One-Dimensional Transport with Equilibrium Chemistry (OTE

2.4.3 Initial Conditions

In addition to the boundary conditions described in Section 2

tion of initial conditions that define the state of the modeled system 

segment of the stream network (fig. 11). This initialization is acco

( , fig. 12) and a steady-state solution of the governing transpor

assumed (Runkel, 1998; physical transport only — in the absence 

steady-state solution neglects the effects of reaction, the assigned i

upstream boundary condition on the modeled system. The model m

ary condition in place prior to reaching a quasi-steady state in whi

transport and reaction. This use of model “spin-up” (Thornton and 

Section 4.4.

C

 ̃ bc
3 User’s Guide

This section provides information on the use of the OTEQ 

described in Section 2.2 based on user-specified model parame

ified items are described in Sections 3.1–3.5.

3.1 Conceptual System, Revisited

Before giving a detailed description of the model’s input re

terms of the conceptual system. Figure 14 depicts the modeled 

as a continuous distance along which the physical model param

constant dispersion coefficient, lateral inflow concentration, an

system reflects both its inherent variability and the availability o

reach, whereas a stream with a well-characterized variation in c

number of reaches in the modeled system is specified by the N

Figure 14. Conceptual system that includes one or more reaches. 

segments.

Each reach is subdivided into a number of computational 

over which the governing mass-balance equations apply. For a 

that DELTAX is determined from the reach length, RCHLEN, 

DELTAX = 

Additional program variables are depicted in figure 15, wh

reach begins at the upstream boundary of the system, an incom

defined. USBC is the upstream boundary condition, which is di

table 16). The program variable denoting the starting stream dis

that these three variables apply only to the first reach of the mo

Reach Reach . . . Re

Segment Segment

. . .

Re

21 

12 

Conceptual System — Re

RCH

DELTAXj
28 One-Dimensional Transport with Equilibrium Chemistry (OTE

Figure 15. The first reach in the conceptual system and the required 

The remaining program variables shown in the figure, QLAT

reach in the modeled system (an exception is the case of unsteady 

resents flow entering the channel through surface inflows, overlan

flow carries a solute concentration, CLATIN. The final variable, Q

from the main channel due to ground-water recharge or surface-wa

that is equal to the simulated main channel concentration. Both QL

[L3 T−1 L−1].

3.2 Input/Output Structure

This section describes the input and output files associated wit

detailed descriptions are provided in Sections 3.3–3.5.

The input/output structure of OTEQ is depicted in figure 16. A

user-specified filenames for several input and output files. The par

conditions, and model parameters that remain constant throughout 

model parameters that can potentially vary during the simulation (

tional area). The control, parameter, and flow files are very similar 

model (Runkel, 1998). The MINTEQ input file contains the chem

tion 3.3.6). Finally, the MINTEQ database files are used to define 

tion 3.3.7).

The filenames for the control file (control.inp) and MINTEQ 

error.dbs) are set within the software and cannot be changed. The 

specified by the user in the control file.

4

Loss of water due to evaporation is not explicitly considered within OTEQ. Im

2004); contact the author of this report for additional details.

Conceptual System 

QSTART

12

x = XSTART

QL

USBC

CL

QL
Figure 16. OTEQ Input/Output files.

Also shown in figure 16 are the output files created by OT

an “echo” of the user-specified simulation options and model p

and any error messages generated during model execution. In a

each solute. If precipitation is being modeled, a precipitate outp

Similarly, when sorption is modeled, one sorption output file is 

described in Section 3.5.1, the solute, precipitate, and sorption o

specified print locations (Section 3.3.4).

Concentration-distance output files are optional output file

profile of the simulated solute concentrations at the completion 

output files are provided in Section 3.5.2.

The filenames for the solute output files are specified by t

output files are created by adding a “.pre” and “.sor” extension t

centration-distance output files are created by adding a “.d” ext

3.3 Input Format

As described in Section 3.2, several input files must be ass

input file is described in terms of a set of record types. Within ea

parameters. In general, record types refer to lines (rows) in the 

record. In most cases, the specific columns for a given field are 

tions to the column requirements are noted in the text, when a giv

columns is not required. In addition to column specifications, e

precision (D), integer (I), or character (C). Double precision fie

integer fields require the entry of a numeric value without a deci

Example input files are presented in Section 4.

Input/Outpu

O

control.inp

flow file

parameter file

MINTEQ

input file

MINTEQ

database

files
30 One-Dimensional Transport with Equilibrium Chemistry (OTE

3.3.1 Units

In Section 2.2 the governing equations and model parameters 

[T]. Within OTEQ, specific units are assigned to the various input 

subject to the following:

• Length Units [L]. Any unit of length may be used when 

the length unit must be consistent for all model parameter

denoted by an L in the input file descriptions that follow.

• Time Units. All model parameters that require the specifi

example) are defined in terms of seconds. Note, however, 

time (TSTART) are specified in terms of hours.

• Concentration Units. Concentration units of moles per lit

(boundary conditions, lateral inflow concentrations, etc.).

3.3.2 Internal Comments

Model users may document their work by placing comments 

pound sign (#) placed in column number 1 will be treated as comm

input files described in Section 4.

3.3.3 The Control File

The control file, control.inp, is used to specify the filenames 

file is described here; an example control file is presented in Sectio

control file consists of four record types. Record type 1 specifies t

name of the MINTEQ input file, record type 3 specifies the name 

solute output files. Although table 1 lists four record types, the con

repeated for each solute (output files are created for each solute).

Table 1. The OTEQ control file.

[C, Character]

Record type Input variable Format Colum

1 FILE C 1–20 

2 FILE C 1–20 

3 FILE C 1–20 

41 FILES C 1–20 

1Record type 4 is used once for each solute modeled (NSOLUTE tim
3.3.4 The Parameter File

The parameter file specifies print options, boundary condi

the run. The format of the parameter file is given in tables 2–16

below.

Record Type 1 — Simulation Title

The first record in the parameter file is a simulation title o

Record Type 2 — Print Step

The print step specifies the time interval at which results a

needed every 15 minutes, a value of 0.25 hour is entered for the 

differ if the requested step is not an even multiple of the integra

step equal to the nearest multiple of the integration time step.

Record Types 3–5 — Time Parameters

The next step in constructing the parameter file is to enter 

TSTEP input variable. TSTEP is the integration time step (Δt) u

performed to determine a time step that yields an accurate solu

The remaining time parameters, TSTART and TFINAL, a

denotes the simulation start time, and TFINAL specifies the simu

and TFINAL are specified in hours.

Record Type 6 — Distance at the Upstream Boundar

Record type 6 specifies XSTART, the distance at the upst

XSTART is used to determine the distances at various locations 

modeled area, XSTART is commonly set to 0.0.

Record Type 7 — Downstream Boundary Condition

The downstream boundary condition, DSBOUND, is set u

sented by the downstream boundary condition is set to zero. As 

the concentration gradient at the downstream boundary is equal 

should be such that the location of the downstream boundary is s

tional details on the downstream boundary condition are provid

Record Type 8 — Relative Error Tolerance

The relative error tolerance, TOL, is specified using recor

gence of the sequential iteration procedure (Section 2.2.3, eq. 1

Record Type 9 — Chemistry Option

As described in Section 2.2.3, the equilibrium submodel is 

chemical composition of water in the main channel. Section 2.3

position of the transient storage zone. OTEQ is configured such 

tations for the storage zone. This option is implemented using t

computations will be restricted to the main channel and chemic

If ICHEM is set to 2, equilibrium computations will be conduc

important to note that the physical process of transient storage i

0.0, record type 11).

Record Type 10 — The Number of Reaches

As discussed in Section 3.1, the modeled system is divide

record type 10. For each reach, spatially constant parameters ar
32 One-Dimensional Transport with Equilibrium Chemistry (OTE

Record Type 11 — Reach Specific Parameters

Record type 11 sets the parameters for each reach (table 3). U

input variable. The first field in the record indicates the number of 

specifies the length of the reach, RCHLEN. These two parameters 

Section 3.1. 

The dispersion coefficient (DISP), storage zone cross-section

(ALPHA), are specified in fields 3–5. The transient storage mecha

Application 1 (Section 4.1). To avoid division by zero, AREA2 m

tion results, however, provided ALPHA is 0.0. 

The coefficient in the depth-area power function (Section 2.3

tion, B1, complete record type 11. Because record type 11 is used 

for each reach in the stream network.

Table 2. The parameter file — record types 1–10.

[C, Character; D, Double precision; I, Integer; L, chosen length unit (Section 3.3.1

Record 

type Input variable Format Column 

1 TITLE C 1–80 — 

2 PSTEP D 1–13 hours 

3 TSTEP D 1–13 hours 

4 TSTART D 1–13 hour 

5 TFINAL D 1–13 hour 

6 XSTART D 1–13 L 

7 DSBOUND D 1–13 (L per sec

8 TOL D 1–13 — 

9 ICHEM I 1–5 — 

10 NREACH I 1–5 — 
Record Types 12 and 13 — Solute Information

Governing equations for the reactive transport model are d

chemical components within the modeled system. With the exc

all components within the equilibrium submodel are subject to 

these transported components are known as solutes.

The number of solutes, NSOLUTE, is specified in record 

LUTE is equal to the number of components in the MINTEQ in

Record type 13 defines each solute (table 5). The first field 

ciated with a given solute. The ID specified here must also be pre

between MINTEQ component numbers (ID) and specific chem

contain flags which indicate whether a given solute is involved 

and (or) reactions involving an external source/sink (EFLAG). A 

reactions, whereas a value of 0 indicates that the solute does not 

that make up the precipitate species defined using record type 15

sorption reactions (see Section 3.4.3). Because record type 13 d

Table 3. The parameter file — record type 11.

[D, Double precision; I, Integer; L, chosen length unit (Section 3.3.1

Input variable Format Column Units 

NSEG I 1–5 — 

RCHLEN D 6–15 L 

DISP D 16–25 L2

 per sec

AREA2 D 26–35 L2 

ALPHA D 36–45 per secon

A1 D 46–55 — 

B1 D 56–65 — 

Table 4. The parameter file — reco

[I, Integer]

Input variable Format C

NSOLUTE I 1–5 
34 One-Dimensional Transport with Equilibrium Chemistry (OTE

Record Type 14 — Number of Precipitates

OTEQ allows for the precipitation of one or more precipitate 

ified using record type 14 (table 6). When precipitation is modeled, 

within the equilibrium submodel (possible solids are defined in th

tation reactions are not modeled, NPRECIPS is set to 0 and record 

Record Type 15 — Precipitate Definition (optional)

Note: Record type 15 is omitted if precipitation is not modeled (N

Precipitate species are defined using record type 15 (table 7). 

fication number of the precipitate species, IDPRECIP. IDPRECIP 

a possible solid within the MINTEQ input file (Sections 3.3.6 and 3

and storage zone settling velocities, PSETTLE and PSETTLE2, re

and 47). Record type 15 is specified for each precipitate species (i

Table 5. The parameter file — record type 13.

[I, Integer; Record type 13 is used once for each solute mod

Input variable Format Column 

ID I 1–5 M

PFLAG I 6–10 Pr

SFLAG I 11–15 So

EFLAG I 16–20 Ex

Table 6. The parameter file — record type 

[I, Integer]

Input variable Format Column 

NPRECIPS I 1–5 

Table 7. The parameter file — record type 15.

[D, Double precision; I, Integer; L, chosen length unit (Section 3.3.1); Reco

Input variable Format Column Units 

IDPRECIP I 1–8 — 

PSETTLE D 9–18 L per second 

PSETTLE2 D 19–28 L per second 
Record Type 16 — Sorption Option

As described in Section 2.3.3, OTEQ may model sorption o

is controlled by the sorption option, ISORB (record type 16, tab

indicates sorption to a static surface, whereas a value of 2 indica

and dynamic surfaces is requested by setting ISORB to 3. Sorp

currently possible.

Record Types 17–20 — Sorption Definition: Static Su

Note: Record types 17–20 are omitted if sorption to a static sur

Record types 17–20 are specified if sorption to a static surf

be in chemical equilibrium or kinetically limited. The kinetic ra

kinetically limited sorption, a kinetic rate coefficient is defined 

LAMBS 

where Γ is the fraction of the equilibrium quantity allowed to s

and 36) and Δt is the integration time step in seconds (TSTEP × 

Record type 18 specifies the sorbent solid concentration as

24 and 25). Record type 18 is specified for each reach (it is use

In record type 13, SFLAG is used to indicate which solute

the static surface, record types 19 and 20 are specified for each 

IDSOLID, the MINTEQ component number. Record type 20 is 

IDSOLID in the main channel (SORBB) and the storage zone (

trations in reach one, record type 20 is repeated for the remaini

record types 19 and 20 is then specified for the remaining solut

specified in the same order as the MINTEQ component numbe

Table 8. The parameter file — reco

[I, Integer]

Input variable Format C

ISORB I 1–5 

Table 9. The parameter file — record types 17–18.

[D, Double precision]

Record 

type Input variable Format Column 

17 LAMBS D 1–13 p

181 SOLCON D 1–13 g

1

Record type 18 is used once for each reach (NREACH times).
36 One-Dimensional Transport with Equilibrium Chemistry (OTE

Record Type 21 — Sorption Definition: Dynamic Surfac

Note: Record type 21 is omitted if sorption to a dynamic surface i

Record type 21 is specified if sorption to a dynamic surface is 

11), IDSORB, specifies the MINTEQ component number of the so

specifies MWSORB, the molecular weight of the dynamic surface

Record Type 22 — Redox Option

Section 2.3.4 presents an iterative procedure whereby oxidati

setting the redox option, IREDOX, to 1 (record type 22, table 12) a

24. Oxidation/reduction will not be modeled if IREDOX is set to 0

Table 10. The parameter file — record types 19–20.

[D, Double precision; I, Integer]

Record 

type

Input 

variable Format Column Un

191

1

The block of record types 19–20 is used for each solute with SFLAG=1.

IDSOLID I 1–5 — 

201,2 SORBB D 1–13 moles p

201,2 SORB2B D 1–13 moles p

2Record type 20 is used once for each reach (NREACH times).

Table 11. The parameter file — record type 21.

[D, Double precision; I, Integer]

Input variable Format Column Units 

IDSORB I 1–7 — 

MWSORB D 8–19 gram per mole 

Table 12. The parameter file — record 

[I, Integer]

Input variable Format Colu

IREDOX I 1–5 
Record Types 23 and 24 — Redox Definition (option

Note: Record types 23 and 24 are omitted if oxidation/reductio

Record type 23 (table 13) is used to specify which solutes 

specify the MINTEQ component numbers of the two solutes in

24 specifies the target fraction used in the iterative procedure, T

variable; record type 24 is therefore specified for each reach (it is 

procedure distributes mass between IDRED and IDRED2 such t

IDRED is equal to THETA (where 0.0 < THETA < 1.0).

Record Types 25 and 26 — Output Specifications

Record types 25 and 26 control the format and type of outp

output the time-series of solute concentration at a number of fix

3.2) and the spatial distribution of solute concentration at a fixe

Within OTEQ, the fixed locations used for time-series out

ifies the number of print locations, NPRINT. Fields 2–4 set the 

print option (PRTOPT). Record type 26 specifies the distance o

If the interpolation option (IOPT) is set to 1, the concentrati

the centers of the two segments closest to the print location. If I

to the concentration of the nearest upstream segment. For many 

ulation results. Selection of the interpolation option (IOPT=1) m

used. Interpolation is not recommended when the print location

ized by physical parameters that differ markedly from those of 

The distance option, DOPT, is used to request the creation 

distance files will be created that contain the solute concentratio

type 5). If DOPT is set to 0, concentration-distance files are not c

a steady-state solution is of interest (see Applications 4 and 5, S

The print option, PRTOPT, determines the format of the o

output for the main channel only. Solute concentrations in both 

is set to 2. 

Table 13. The parameter file — record types 23 and 24.

[D, Double precision; I, Integer]

Record 

type Input variable Format Column 

23 IDRED I 1–5 

23 IDRED2 I 6–10 

241 THETA D 1–10 

1

Record type 24 is used once for each reach (NREACH times).
38 One-Dimensional Transport with Equilibrium Chemistry (OTE

Record Types 27 and 28 — Upstream Boundary Condit

The final record types in the parameter file specify the time-v

number of boundary conditions (NBOUND) and the boundary con

NBOUND upstream boundary conditions are characterized by a s

type 28. Each field in record type 28 has a free format such that pl

at which a boundary condition goes into effect (USTIME) is speci

denotes the boundary values corresponding to USTIME. The USB

boundary value for solute one is placed to the right of USTIME, th

for solute one, etc.). After specifying USTIME and USBC for the 

subsequent change in the boundary condition (it is used NBOUND 

Table 14. The parameter file — record types 25–26.

[D, Double precision; I, Integer; L, chosen length unit (Section 3.3.1)]

Record 

type Input variable Format Column 

25 NPRINT I 1–5 — 

25 IOPT I 6–10 — 

25 DOPT I 11–15 — 

25 PRTOPT I 16–20 — 

261 PRTLOC D 1–13 L 

1Record type 26 is used for each print location (NPRINT times).

Table 15. The parameter file — record type 27.

[I, Integer]

Input variable Format Column Units 

NBOUND I 1–5 — 

IBOUND I 6–10 — 
Interpretation of USBC is dependent on the type of upstre

tions may be specified using the boundary condition option, IB

imposed at the upstream boundary for IBOUND equal to 1. Un

centration, Cbc, described in Section 2.4.2. With a step concentra

and subsequently updated (to USBCj+1, USBCj+2,...) at the app

Figure 17. Upstream boundary condition options. The upstream b

step flux profile, or a continuous concentration profile.

Table 16. The parameter file — record type 28, upstrea

[D, Double precision; L, chosen length unit (Section 3.3.1); Record 

Input variable Format Units1

1

Units for USBC are dependent on the value of IBOUND in reco

USTIME D hour 

USBC2 D moles per liter; 

(moles per liter)(L3

 per 

2

USBC repeats horizontally (the value for solute one is placed to 

USBC for solute one, etc.). 

C @ x=XSTART

QC @ x=XSTART

C @ x=XSTART

Upstream Boundary 

Concentration — Step

IBOUND = 1

Flux — Step

IBOUND = 2

Concentration — Continuous

IBOUND = 3
40 One-Dimensional Transport with Equilibrium Chemistry (OTE

A step flux profile is imposed at the upstream boundary for I

QCbc, where Q is the volumetric flow rate at the upstream bounda

the first boundary value (USBCj) is divided by Q to obtain an initi

sequently updated at the appropriate times (USTIMEj+1, USTIME

sented by Runkel (1998).

A continuous concentration profile is imposed at the upstream 

updated at each model time step. The value assigned to Cbc is dete

To allow for interpolation, the time of the last boundary condition 

lation end time, TFINAL (record type 5). An example of a continu

(1998).

3.3.5 The Flow File

The flow file defines the model parameters that can potential

rate (Q, QSTART), lateral flow rates (QLATIN, QLATOUT), main 

concentration (CLATIN). Collectively, these parameters are know

The format of the flow file depends on the nature of the param

with respect to time, the flow is steady, whereas if any of the flow v

For steady flow conditions, the flow variables are specified on 

flow variables are specified using flow locations. The formats of t

The Flow File — Steady Flow

When the flow variables are constant with respect to time, the 

tions, the specification of a steady flow file is fairly straightforwar

Record Type 1 — Change-in-Flow Indicator. The first record ty

record type is used to define the flow regime. As shown in table 17

ables. Here we are concerned with a steady flow regime in which th

Record Types 2 and 3 — Flow Variables. The remaining record 

Record type 2 sets the volumetric flow rate at the upstream bounda

type 3, where each field has a free format such that placement of i

the first three fields of this record specify the lateral flow rates (QLA

a particular reach. The fourth field, used to indicate the lateral inflo

solute modeled. Because the flow variables vary from reach to reach

Table 17. Steady flow file — record type 1.

[D, Double precision]

Input variable Format Column Units 

QSTEP D 1–13 hour 

Table 18. Steady flow file — record type 2.

[D, Double precision; L, chosen length unit (Section 3.3.1)]

Input variable Format Column Units 

QSTART D 1–13 L3

 per seco
The Flow File — Unsteady Flow

With the allowance for time-variable flow parameters com

unsteady flow file differs markedly from that presented for the 

its steady counterpart, as several record types must be repeated 

The format used here is compatible with standard routing m

son, 1989). As discussed below, record types 2–8 allow for the s

the stream channel. These record types may be created by refor

Record Type 1 — Change-in-Flow Indicator. For an unstead

changes in the flow values occur. QSTEP therefore defines the t

change every 15 minutes, for example, QSTEP is set to 0.25 ho

Record Types 2 and 3 — Flow Locations. When detailed field 

varying flow variables at various locations along the stream. The

routing models. The format of the flow file allows the user to d

specified. These flow locations often correspond to the reaches 

required, however.) Using the data provided at the flow location

for each segment in the stream network. Lateral inflow rates an

locations.

Record type 2 indicates the number of flow locations, NFL

This record is used NFLOW times, once for each flow location

locations. The requirements, listed below, are checked internall

• The flow locations must be entered in ascending (down

first, the third must be downstream from the second, et

• The first flow location must be placed at the upstream 

XSTART, where XSTART is the starting location spec

location is analogous to the QSTART parameter in the 

• The last flow location must be placed at or below the d

Table 19. Steady flow file — record type 3.

[D, Double precision; L, chosen length unit (Section 3.3.1); 

(NREACH times)]

Input variable Format Units 

QLATIN D L3 second−1

 L−1 

QLATOUT D L3 second−1

 L−1 

AREA D L2 

CLATIN1

1

CLATIN repeats horizontally (the lateral inflow concen

centration for solute two is placed to the right of CLATIN fo

D moles per liter 

Table 20. Unsteady flow file — record type 1.

[D, Double precision]

Input variable Format Column Units 

QSTEP D 1–13 hour 
42 One-Dimensional Transport with Equilibrium Chemistry (OTE

Record Types 4–7 — Lateral Flows, Flows, and Areas. In contr

(QLATIN, CLATIN) do not necessarily correspond to specific reac

each flow location. The value specified is used for all segments in 

diately upstream (from location j-1 to location j). Note that this sch

tions are placed at the end of each reach. As stated earlier, the flow 

interpolate values for the segments within the network.

Values for lateral inflow, flow rate, area, and lateral inflow co

7, respectively. Record types 4–7 all use a free format such that pl

in table 22, the input fields in these record types (QLATIN, Q, AR

values for flow location 1 appear first, while those for flow locatio

7 is used once for each solute (record type 7 is used NSOLUTE tim

Record types 4–7 repeat for each change in the flow variable

(QSTEP = 0.25 hour), these record types will repeat four times for 

unsteady flow file may contain a large number of records.

Table 21. Unsteady flow file — record types 2 and 3.

[D, Double precision; I, Integer; L, chosen length unit (Section 3.3.1)]

Record 

type Input variable Format Column 

2 NFLOW I 1–5 — 

31 FLOWLOC D 1–13 L 

1

Record type 3 is used once for each flow location (NFLOW times).

Table 22. Unsteady flow file — record types 4–7.

[D, Double precision; L, chosen length unit (Section 3.3.1)]

Record 

type Input variable Format Units 

41,2

1

The block of record types 4–7 repeats for each change in the flow va

2

QLATIN, Q, AREA, and CLATIN repeat horizontally for each flow 

ues for location 2 appear to the right of those for flow location 1, etc.).

QLATIN D L3 second−1 L−

51,2 Q DL3

 per second 

61,2 AREA D L2 

71,2,3

3

Record type 7 is used once for each solute (NSOLUTE times).

CLATIN D moles per liter 
3.3.6 The MINTEQ Input File

The MINTEQ input file defines the chemical components

rium submodel. Unlike the other input files, the MINTEQ input 

file is created using PROTEQ, an adaptation of the problem defi

others, 1991). Because the MINTEQ input file is generated by 

presented here. Additional information on the use of PROTEQ i

4.2.2, 4.3.3, and 4.5.1), and by Allison and others (1991).

Record Types 1–4 — Temperature and Ionic Strength 

Equilibrium constants used within the equilibrium submod

Database values are referenced to 25° C and zero ionic strength 

types 1–4 (table 23) are used to specify temperature, ionic stren

The first record type is used to specify the water temperatu

corrected equilibrium constants based on a power function or th

temperature is spatially and temporally constant, and TEMPC s

Record type 2 is used to specify the ionic strength option, I

librium submodel; if ISOPT is set to 1, fixed ionic strength is a

(FIONS) is specified using record type 3. Fixed ionic strength is 

putational demands of the equilibrium submodel. In addition, re

ute to ionic strength need not be included as transported solutes

poral variation in ionic strength is not considered; as with TEM

The activity coefficient option, KKDAV, is specified usin

to develop activity coefficient corrections for the equilibrium co

As noted by Allison and others (1991), Debye-Huckel paramete

parameters are unavailable, the Davies equation is used.

Record Type 5 — Maximum Iterations

Record type 5 is used to specify ITMAX, the maximum n

model.

Record Type 6 — Sorption Option

The sorption option, IADS, is specified using record type 6

rithms (Allison and others, 1991), only the “diffuse layer mode

ymous with the “generalized two layer model” described in Sec

file record type 16), IADS is set to 0. When sorption is modeled 

invoked.

Record Types 7 and 8 — Surface Definition for Sorpt

Note: Record types 7 and 8 are omitted if sorption is not mode

As described in Section 2.3.3, sorption may involve static 

NADS, is specified using record type 7. When only a single sur

NADS is set to 1. When both static and dynamic surfaces are co

specific surface area of each surface, SSA, and is used NADS t

Record Type 9 — Component Definitions

Record type 9 is used to specify the component identificat

guess (GX) for each component. Record type 9 is used once for 

13) and additional times for the sorption components (Section 2

Additional Record Types

In addition to the record types presented above, several ot

record types include those needed to specify possible solids, ga

MINTEQ database files. As with the other record types, these a

ware, PROTEQ (see Section 3.4.3). Examples of these addition
44 One-Dimensional Transport with Equilibrium Chemistry (OTE

Table 23. The MINTEQ input file — record types 1–6.

[D, Double precision; I, Integer]

Record 

type Input variable Format Column 

1 TEMPC D 1–5 degrees 

2 ISOPT I 1–5 — 

3 FIONS D 1–9 moles p

4 KKDAV I 1–5 — 

5 ITMAX I 1–5 — 

6 IADS I 1–5 — 

Table 24. The MINTEQ input file — record types 7 and 8.

[D, Double precision; I, Integer]

Record 

type Input variable Format Column 

7 NADS I 1–5 — 

81 SSA D 1–13 mete

1Record type 8 is used once for each surface (NADS times).

Table 25. The MINTEQ input file — record type 9.

[D, Double precision; I, Integer; Record type 9 is used once for each com

Input variable Format Column Units1 

IDX I 1–7 — 

T D 9–18 moles per liter 

GX D 20–26 log(moles per 

1

Units for T and GX are in terms of moles of sites per gram of sorbent 
3.3.7 The MINTEQ Database Files

The database files distributed with OTEQ are based on the 

others, 1991). The version 3 files have been updated to provide 

EQC (Parkhurst and Appelo, 1999), as described in Appendix 

information is provided by Allison and others (1991).

Required Database Files

Several of the MINTEQ database files are required for OT

comp.dbs, error.dbs, and analyt.dbs. Three of these files (therm

thermodynamic data (enthalpy values and equilibrium constant

the equilibrium submodel are contained in error.dbs. A final da

used to develop temperature-corrected equilibrium constants.

Auxiliary Database Files

Several auxiliary database files are used by the problem de

reactions. Databases gases.dbs and redox.dbs are used to defin

Two additional database files are used to define sorption re

includes 40 reactions between aqueous components and a solid 

1990). These reactions may be added to the MINTEQ input file 

dlm2.dbs, is identical to feo-dlm.dbs, except that the reactions 

surface are added to the MINTEQ input file when concurrent s

ISORB=3, parameter file record type 16). Use of the sorption d

3.4 Input File Preparation and Model Execution

Instructions for routine execution of the OTEQ solute tran

vided in this subsection. These instructions assume that the user 

5.3). Prior to running OTEQ, the control, parameter, flow, and 

below.

3.4.1 Preparation of the Control File

The control file is created using a text editor (vi or emacs, 

(table 1).

3.4.2 Preparation of the Parameter and Flow Fil

As with the control file, the parameter and flow files are pr

files are described in Sections 3.3.4 and 3.3.5 (tables 2–22). Wh

tion of the various record types is generally straightforward. Ad

or when using alkalinity to determine dissolved inorganic carbo

flow file record type 3 (steady flow), and flow file record type 

MINTEQ runs. Use of MINTEQ to develop these record types i

described below.

Prior to executing MINTEQ in stand-alone mode, an inpu

distributed with MINTEQ (Allison and others, 1991). The PRO

5.3.3) by entering the prodef command at the shell prompt:

prodef

This command will start PRODEF; the user then creates an inpu

file is then used to conduct a stand-alone MINTEQ run. This ru

command at the shell prompt:

minteq

When MINTEQ execution is complete, the output file is inspect

4.2–4.5.
46 One-Dimensional Transport with Equilibrium Chemistry (OTE

3.4.3 Preparation of the MINTEQ input file — Use 

The MINTEQ input file used with OTEQ is typically generat

definition software (PRODEF) distributed with MINTEQ (Allison 

guided through a series of text-based prompts and menus. At the c

input file is created.

The PROTEQ session is initiated from the OTEQ work area 

prompt:

proteq

This command will start PROTEQ and bring up a prompt that ask

filename specified should be the same as that provided in the contr

acters. After specifying the MINTEQ input filename, the user may 

template or “seed” file is used if the user specifies a filename at th

from scratch, a simple carriage return is entered at the seed file pro

menu known as Edit Level I.

Edit Level I

Edit Level I is used to set record types 1–5 in the MINTEQ inp

options and the corresponding default values. For most OTEQ appl

Menu item 1 is used to specify temperature (TEMPC, MINTEQ in

strength (ISOPT and FIONS, MINTEQ input file record types 2 a

values for most applications.

Edit Level II, Menu Item 1 — Specify Aqueous Compon

After leaving Edit Level I, the user is presented with a main m

1 of Edit Level II allows the user to select aqueous components an

MINTEQ input file is developed for each of the modeled solutes.

Aqueous components are selected by responding to a series o

prompted for the component’s total concentration. Although this t

PROTEQ to develop an initial log activity guess. The total concentr

at the upstream boundary of the modeled system (Section 3.3.4, rec

solutes defined within the parameter file (Section 3.3.4, record typ

The remaining PROTEQ menus are not required for the most 

when precipitation and(or) adsorption reactions are considered, as 

Edit Level II, Menu Item 8 — Specify Possible Solids (op

If precipitation is being modeled, a possible solid is specified 

(Section 3.3.4, record type 15). Possible solids are specified by se

Edit Level II, Menu Item 3 — Specify Adsorption Model 

If sorption is being modeled (ISORB>0, parameter file recor

using menu item 3 of Edit Level II. The first step is to select the diff

After selecting the diffuse layer model, sorption surfaces and reac

Defining Sorption Surfaces. Depending on the application, 

dynamic surface is considered (ISORB=1 or 2), only a single surf

considered (ISORB=3), two surfaces are defined. The following c

• For a given sorption surface, the first site is the high-affini

• If concurrent sorption onto static and dynamic surfaces is 

and surface number two is the dynamic surface.

A sorption surface is added using menu item 1 from the sorption op

record type 8) of the surface is entered first. The user is then asked 

affinity site, site type 1), where the concentration of site type 1 is e

the sorbent (M).5 A second site type (low-affinity site, site type 2) i

The process of defining a surface with two site types is repeated if a 

a surface composed of hydrous ferric oxide (HFO; Dzombak and M
Defining Sorption Reactions. The final step in setting up 

between the aqueous components and the sorption surface. OT

others, 1999; Runkel and Kimball, 2002) have assumed that the 

Dzombak and Morel (1990). The instructions which follow are 

rium submodel. As discussed by Allison and others (1991), oth

Sorption reactions for the HFO surface are added using m

HFO database (feo-dlm.dbs) to be added to the MINTEQ inpu

submodel. If two surfaces are modeled (ISORB=3, static and d

database for the second surface (feo-dlm2.dbs).

3.4.4 Execution of OTEQ

Before running OTEQ, the required input files (control, pa

OTEQ work area (Section 5.3.3). To run OTEQ, enter the follo

oteq &

Use of the ampersand (&) after the command name causes the 

to the current window during model execution. Status of the ot

ation (Sections 2.2.3 and 2.3.4) and time step information is pe

tion. Users may inspect Part V of echo.out during program exe

or more echo.out, for example).

As discussed in Section 5.5.3, OTEQ error messages are w

the model run, users should inspect echo.out to verify model in

inspect the echo file for the following items:

• Verification of User Input. Users should inspect echo

the input files have been read correctly by OTEQ. User 

rectly by OTEQ, resulting in erroneous simulation resu

• Verification of the Flow Distribution. OTEQ sets the 

supplied information from the flow file (QSTART, QLA

Volumetric flow rates at the user-specified print locatio

used to verify that the spatial flow distribution has been 

• Verification of θ. The iterative procedure used to imple

yield the target fraction specified by the user (θtarget; TH

of echo.out to see how closely the simulated fraction (

locations.

• Warnings and Error Messages. Additional warnings 

5

This value (NS/M) is used to set T in record type 9 of the MINTEQ input fi

total concentration of the sorption component as given by equation 25 (Section 

Table 26. Sorption parameters — hydrous ferric oxide (

Parameter 

Specific surface area, SSA [meter2

 per gram HFO] 

Molecular weight of HFO, M [gram HFO per mole HFO] 

Site density, NS, low-affinity site [moles sites per mole HF

Site density, NS, high-affinity site [moles sites per mole HF
48 One-Dimensional Transport with Equilibrium Chemistry (OTE

3.5 Output Analysis

This section describes the mechanics of interpreting the solut

and sorption), and concentration-distance output files. These files 

and 42, Section 2) can be calculated from the available output. The 

mats (Sections 3.5.1 and 3.5.2), utility programs to reformat the ou

ulation results (Section 3.5.4).

3.5.1 The Solute and Solid Output Files

As discussed in Section 3.2, the OTEQ solute transport model 

files (precipitate and(or) sorption output files) are also created for s

(NPRECIPS>0 and(or) ISORB>0, Section 3.3.4). The following pa

In the discussion that follows, each line is described in terms of a 

 The format of the solute output file is illustrated in figure 18. E

(mobile) solute concentrations at the user-specified print locations

T 

w T Pb – Sb – 

= = 

in the main channel, and

T 

ws Ts Psb – Ssb – 

= = 

in the storage zone (figs. 1 and 10, Section 2). The first field of eac

(fig. 18). The next NPRINT fields give the total waterborne solute 

locations (PRTLOC, Section 3.3.4). If the print option (PRTOPT, 

waterborne solute concentrations in the storage zone (Tws) at the N

Figure 18. Example solute output file.

The solid output files include a time series of total solid conce

concentrations (immobile precipitated or immobile sorbed). Water

sorbed) may be obtained by difference. The format of the precipit

line is the simulation time, in hours (fig. 19). The next NPRINT fie

(P) at the NPRINT print locations; a second set of NPRINT fields g

(Pb) at the NPRINT print locations. Two additional sets of NPRINT 

for the total precipitated concentration in the storage zone (Ps) and 

tions in the storage zone (Psb) (these additional sets of fields are n

ogous to the precipitate output file shown in figure 19; simulated c

sorbed concentration in the main channel (S), immobile sorbed con

in the storage zone (Ss), and immobile sorbed concentration in the 

6

When modeling pH, one of the solutes represents the component associated wit

in the main channel and storage zone are output rather than Tw and Tws.

Time [hour]

Solute Out

Main Channel, Tw

@PRTLOC1 @PRTLOC2... @PRTLO

Total Waterbo

1.3000000000E+01 2.5514970966E-06 1.9435631478E-06 2.753164

1.3010000000E+01 2.5520255698E-06 1.9423423860E-06 2.753301

1.3020000000E+01 2.5525415204E-06 1.9411901809E-06 2.753437

1.3030000000E+01 2.5530448754E-06 1.9400188841E-06 2.753572

1.3040000000E+01 2.5535358426E-06 1.9387703773E-06 2.753707

1.3050000000E+01 2.5540149403E-06 1.9374143319E-06 2.753841

1.3060000000E+01 2.5544830222E-06 1.9359409690E-06 2.753974

1.3070000000E+01 2.5549413879E-06 1.9343511318E-06 2.754107
Figure 19. Example precipitate output file, for the case of PRTOPT

3.5.2 Concentration-Distance Output Files

Concentration-distance output files are generated when th

contain the spatial profiles of the concentrations included in the 

ulation. The format of each concentration-distance output file is 

and 19), with the exception of column one where distance is ou

tain one line for each segment in the modeled system.

3.5.3 The Post-Processor, POSTEQ

The solute, solid, and concentration-distance output files de

a simple post-processor known as POSTEQ. Using the solute a

files for each solute at each print location. Separate distance ver

tion-distance files are available (DOPT=1, Section 3.3.4). Thes

data using plotting utilities such as Xgraph and Grace (Section 

To see a list of available post-processing options, type the 

posteq -h

Post-processing is initiated by typing:

posteq

Storage zone solute concentrations are included in the post-proc

vided PRTOPT=2, Section 3.3.4).

3.5.4 Plotting Alternatives

Model users have several alternatives for plotting results fr

post-processor described in Section 3.5.3 and the Xgraph plotti

described here. For the example, assume the name of the solute 

at a distance of 100. In this case, the post-processor creates an o

for the remaining print locations). Simulation results and obser

xgraph -bb -tk -m so4.out.100.xgr obs10

where obs100.dat contains the observed sulfate data. Other plo

mann.ac.il/Grace) and various spreadsheet programs (the outpu

Time [hour]

Precipit

Total Precipitate Main Channel

@PRTLOC1 @PRTLOC2... @PRT

1.3000000000E+01 4.3909141799E-05 5.6662887202E-05 0.00000

1.3010000000E+01 4.4017580688E-05 5.6981759332E-05 0.00000

1.3020000000E+01 4.4126020648E-05 5.7298998783E-05 0.00000

1.3030000000E+01 4.4234462044E-05 5.7615024805E-05 0.00000

1.3040000000E+01 4.4342905272E-05 5.7930179528E-05 0.00000

1.3050000000E+01 4.4451350471E-05 5.8244723377E-05 0.00000

1.3060000000E+01 4.4559797437E-05 5.8558839074E-05 0.00000

1.3070000000E+01 4.4668245812E-05 5.8872648588E-05 0.00000
50 One-Dimensional Transport with Equilibrium Chemistry (OTE

4 Model Applications

In this section, several applications of OTEQ are presented. E

tures of the model. The reader is encouraged to review all of the app

capabilities. Additional applications of the model are presented els

Runkel and others, 1999; Runkel and Kimball, 2002; Runkel and 

Each application has four basic components: (1) a description 

a brief problem statement, (3) example input files, and (4) simulat

of the input requirements presented in Section 3.3. These example 

used within the various files. Due to space limitations, not all of th

input files is available as part of the software distribution (Section 

meters as the unit of length (L; Section 3.3.1).

4.1 Application 1: Time-Variable Simulation of a 

The primary purpose of this application is to illustrate the for

some of the numerical issues associated with the use of OTEQ.

Chapman (1982) presents a hypothetical example in which a 

fied form of Chapman’s example is considered here, in which a 200

uniform flow and a cross-sectional area of 0.4 meter2. The volume

ficient is 2.0 meter2 per second. Advection and dispersion are the 

gible (α = 0). 

Prior to the injection period, the stream is in equilibrium with 

no solid phases are present on the bed. At time equal to 0.01 hour, 

upstream boundary (x=0) is 0.3 mole per liter. At time equal to 0.0

at 0.09 hour, a second injection begins in which Na2SO4 is injecte

liter. This second injection ends at 0.11 hour.

The double pulse injection is depicted in figure 20. Due to the 

arated in the upstream regions of the stream. As the pulses travel d

of the two pulses occurs. As a result of this mixing, the solution be

(CaSO4⋅2H2O). A precipitate therefore forms and is subject to settl

streambed, and gypsum that has settled to the streambed redissolv

ification of four solutes (Ca, Cl, Na, and SO4) and a solid phase (g

 

Figure 20. Upstream boundary condition for double pulse injection. U

at background (0.0002 mole per liter) during nonpulse time periods (no

0.01 0.03 0.05 

0.3

0.0001

Time in 

Concentration at 

x=0, in moles per lite

r

CaCl2 

Double Pulse 
4.1.1 The Control File — Application 1

Example input files corresponding to the double pulse inje

the names of the parameter, MINTEQ, flow, and solute output f

params.inp, parameters used within the equilibrium submodel a

(record types 1–3). Simulated solute concentrations are written 

ute). Control files for the remaining applications (Sections 4.2–4

discussed further.

Figure 21. Control file for Application 1.

4.1.2 The Parameter File — Application 1

Figures 22 and 23 depict the parameter file for the double pu

every 0.002 hour (record type 2). The integration time step is 0

at 0.7 hours (record types 3, 4, and 5). The upstream boundary o

stream boundary condition (record types 6 and 7). The relative e

and storage zone chemistry is not modeled (record types 8 and 

The stream is divided into two reaches (record type 10). R

specific parameters (dispersion coefficient, storage zone area, e

length of reach one is specified such that the reach ends at the s

an additional 100 meters downstream to reduce any error introd

and 3.3.4). The storage zone exchange coefficient is set to zero 

storage zone cross-sectional area (AREA2) is set to a non-zero 

affect simulation results, provided ALPHA equals zero). For th

gular cross section (B1=1.0) that is 1.33 meters wide (A1=1/wi

7

For this and the subsequent examples, meters are used as the unit of length 

Control File for 

Record Type

1

2

3

4

###################################

#

# OTEQ control file

#

#

# line name of the:

# ---- --------------

# 1 parameter file

# 2 MINTEQ input file

# 3 flow file

# 4 to 3+NSOLUTE solute output files

#

###################################

params.inp

minteq.inp

q.inp

cl.out

na.out

ca.out

so4.out

4

4

4
52 One-Dimensional Transport with Equilibrium Chemistry (OTE

Figure 22. Parameter file for Application 1, record types 1–15.

Parameter File for Ap

Record Type

1

2

3

4

5

6

7

8

9

11

11

14

15

12

13

#######################################

# OTEQ parameter file

#######################################

#

# Parameters: TITLE followed by miscellaneo

# right of field.

#

# Status : Required (these fields are always 

#

#######################################

Double pulse of calcium chloride and sodium s

0.002 | PSTEP [hour]

0.001 | TSTEP [hour]

0.00 | TSTART [hour]

0.70 | TFINAL [hour]

0.0 | XSTART [L]

0.0 | DSBOUND [(moles per liter)(L per 

0.001 | TOL

1 | ICHEM [1=w/o S.Zone chem, 2=w/

2 | NREACH

#######################################

#

# Parameters: Reach definition & physical pa

# Status : Required (these fields are alway

# Note : NS = NSEG

#

# Parameters for each reach (J=1,NREACH):

#

#NS RCHLEN DISP AREA2 ALPHA 

# | | | | 

#######################################

200 200. 2.0 0.05 0.00 

100 100. 2.0 0.05 0.00 

#######################################

#

# Parameters: Solute definition. NSOLUTE f

# for each solute (J=1,NSOLUTE)

# Status : Required (these fields are alway

#

#######################################

 4 NSOLUTE

180 0 0 0 (ID, PFLAG, SFLAG, EFLA

500 0 0 0

150 1 0 0

732 1 0 0

#######################################

#

# Parameters: Precipitate definition. NPREC

# parameters for each precipitate 

# Status : NPRECIPS required. Additiona

# precipitated species required fo

# 

#NPRECIPS

#IDPREC PSETTLE PSETTLE2

# | |

#######################################

1

6015001 2.0e-4 0.0

10

13

13

13
Figure 23. Parameter file for Application 1, record types 16–28.

Four solutes are modeled (record type 12); record type 13 i

various flags for each solute. Specification of record type 13 for 

files in record type 4 of the control file (Cl, Na, Ca, and SO4, re

SO4 to indicate their role in the formation of the gypsum (CaSO

to zero (sorption reactions and external source/sinks are not con

(NPRECIPS=1) and 15. Sorption and oxidation/reduction react

22); record types 17–21 and 23–24 are therefore omitted.

Simulation results are printed at two locations, 100 and 20

boundary concentration is in terms of a step concentration prof

NBOUND=5, record type 27). The first, third, and fifth bounda

absence of an injection; the second and fourth boundary condit

(record type 28; fig. 20). Upstream boundary concentrations (U

(there are two moles of Cl in each mole of CaCl2, such that the 

Parameter File for Appl

Record Type

16

25

26

26

27

#####################################

#

# Parameters: Sorption definition.

# Status : ISORB required. Additional l

#

# ISORB [0=no sorption,1=static,2=dynam

#####################################

0

#####################################

#

# Parameters: Redox definition

# Status : IREDOX required. Additiona

#

# IREDOX [0=no redox,1=based on % dis

#

#####################################

0

#####################################

#

# Parameters: Print definition.

# Status : Required (these fields are alw

# Note : NPR = NPRINT

#

#NPR IOPT DOPT PRTOPT

# | | |

#####################################

2 0 0 1

100.0 (PRTLOC for I = 1, NPRIN

200.0

#####################################

#

# Parameters: Upstream Boundary Condit

# Status : Required (these fields are alw

# Note : NBN = NBOUND

#

#NBN IBOUND

# |

#####################################

5 1 

0.00 0.0002 0.0002 0.0001 0

0.01 0.6000 0.0002 0.3000 0

0.03 0.0002 0.0002 0.0001 0

0.09 0.0002 0.6000 0.0001 0

0.11 0.0002 0.0002 0.0001 0

#

# Cl Na Ca S

22

28

28

28

28

28
54 One-Dimensional Transport with Equilibrium Chemistry (OTE

4.1.3 The Steady Flow File — Application 1

The flow file for the double pulse injection is shown in figure 

3.3.5), and flow at the upstream boundary is 0.1 meter3

 per second 

the flow variables for each reach. Flow is uniform (QLATIN=0), a

reaches. Lateral inflow concentrations (CLATIN) for each solute a

consequence when there is no inflow (QLATIN=0).

Figure 24. Steady flow file for Application 1.

4.1.4 The MINTEQ Input File — Application 1

The MINTEQ input file created using PROTEQ (Section 3.4

using PROTEQ for the double pulse injection generally follows th

set record types 1–5. Record type 2 sets ISOPT equal to zero so th

(examples of fixed, rather than computed, ionic strength are prese

aqueous components (record type 9) corresponding to the four trans

equal to the initial concentrations at the upstream boundary (recor

(IDX, record type 9) for the aqueous components must match the 

13, parameter file). After specifying the aqueous components, men

2H2O) as a possible solid. The identification number of the possib

the precipitate definition (IDPRECIP, record type 15, parameter fi

8

In addition to the four transported solutes, a fifth component, H, is automaticall

independent of pH, the addition of H is not of consequence; H is therefore not defi

simulations (where H is a transported component) are provided in Sections 4.2–4.5

Steady Flow File for 

Record Type

1

3

2

3

#######################################

#

# OTEQ steady flow file

#

#######################################

0.0 QSTEP [hour]

0.1 QSTART [L^3 per second]

#######################################

# for I = 1, NREACH

#

#QLATIN QLATOUT AREA (CLATIN J=1

# 

#######################################

0.0 0.0 0.4 0.0 0.0 0.0 0.0

0.0 0.0 0.4 0.0 0.0 0.0 0.0
Figure 25. MINTEQ input file for Application 1.

4.1.5 Simulation Results — Application 1

Results from the double pulse injection are shown are sho

print location (200 meters) are plotted versus time. Because of t

to the sulfate pulse (fig. 26). Shortly after the peak of the calciu

tion is oversaturated with respect to gypsum due to elevated conc

precipitates and mobile precipitate is present in the water colum

in a gradual increase in immobile precipitate (Pb, fig. 26). As th

urated with respect to gypsum, and the calcium and sulfate that 

hours).

MINTEQ Input File 

Record Type

1

2

6

#####################################

#

# MINTEQ input file for use with OTE

#

#####################################

25.00 | TEMPC

 0 | ISOPT

0.00 | FIONS

 1 | KKDAV

 40 | ITMAX

#####################################

#

# Sorption info - Surface Definition

#

#####################################

 0 | IADS

#####################################

#

# Components

#

#####################################

 330 0.000E+00 -7.00 /H+1

 180 0.600E+00 -0.22 /Cl-1

 500 0.600E+00 -0.22 /Na+1

 150 0.300E+00 -0.52 /Ca+2

 732 0.300E+00 -0.52 /SO4-2

#####################################

#

# Possible Solids

#

#####################################

 

 5 1

6015001 4.8480 -0.2610 /GYPS

3

4

5

9

9

9

9

9
56 One-Dimensional Transport with Equilibrium Chemistry (OTE

Figure 26. Simulated concentrations of calcium (Ca) and sulfate (SO

4.1.6 Numerical Issues — Application 1

The accuracy of the simulation results presented above is hig

niques. As discussed in Sections 2.2 and 2.4, finite-difference techn

ing transport. Numerical aspects of finite-difference techniques (s

and are not discussed herein. Model users who are unfamiliar with 

provided by Thomann and Mueller (1987), Chapra and Canale (19

cific to OTEQ are discussed here. First, a relatively small time ste

centration fronts associated with the double pulse injection (fig. 20

tial iteration convergence failures (Section 2.2.3) that affect the acc

V of the echo output file when convergence failures occur. In gene

dent. Steady-state applications such as those presented in Sections 

variable simulations (Sections 4.1–4.3), for example.

A second numerical issue concerns the treatment of the down

reach 2 is a dummy reach 100 meters in length that is included suc

of interest (the last print location is at 200 meters; the downstream 

meters). Figure 27 illustrates the effect of the downstream bounda

downstream boundary too close to last print location. As shown in t

in a peak total waterborne (dissolved plus mobile precipitate) sulf

reach 2 is reduced to 50 meters, visually identical results are obtain

result in errors in the peak concentration (fig. 27). As with the inte

dummy reach is application dependent. Because of the presence of 

cations with large dispersion coefficients will tend to be sensitive to 

model runs with differing dummy reach lengths should be conduct

condition, as shown here.

Time, in 

Concentration, in moles per liter

0.01

0.02

0.03

0.0

0.04

0.1 0.2 0.3 

0.00

Double Pulse Injectio
Figure 27. Peak total waterborne (dissolved plus mobile precipita

for reach 2. The peak concentration is overestimated when reach l

4.2 Application 2: Time-Variable Simulation of 

In the foregoing example, a relatively simple geochemical 

not affected by pH. A more complex problem is considered here

purpose of the example is to illustrate how OTEQ can be used to 

solute (Section 2.3.1). As shown below, simulation of pH requi

upstream boundary, a quantity that is obtained from a stand-alo

also used to determine the total inorganic carbon concentration 

this example is the modeling of oxidation/reduction (Section 2.

Broshears and others (1996) describe a pH-modification e

Leadville, Colorado. During the experiment, a concentrated solu

5.6 hours. The experiment was designed to raise the acidic water

injection rate resulted in a downstream pH of approximately 4.2

pH (> 5.5). During this second stage of the injection, hydrous o

dissolved concentrations decreased. Decreases in dissolved iron 

oxidation of ferrous iron at the elevated pH. The pH of St. Kev

minated, and precipitated mass residing on the streambed rediss

Kevin Gulch are provided elsewhere (Kimball and others, 1992

and others, 2001).

OTEQ simulations of geochemical systems such as St. Kev

to the complete suite of observed cations and anions. Inclusion o

species, and solution of the coupled transport/equilibrium prob

to reduce simulation run time, OTEQ simulations often utilize 

interest. Prior to finalizing results, a final OTEQ simulation that 

use of the reduced set. In the case considered here, the set of so

chemical reactions. Stand-alone MINTEQ runs, for example, in

tion (concentrations of the uncomplexed species, Na+ and Cl−, 

tion, reactive solutes such as copper do not appreciably affect t

set of solutes considered in this problem is therefore limited to 

iron (Fe(II)), ferric iron (Fe(III)), aluminum (Al), and sulfate (S

Effect of Downstream 

Time, 

Concentration, in moles per liter

0.035

0.030

0.26

0.040

0.28 0.30 
58 One-Dimensional Transport with Equilibrium Chemistry (OTE

4.2.1 The Parameter and Flow Files — Application 

Record types 1–12, 16–21, and 25–26 of the parameter file are 

to Application 1 and the User’s Guide (Section 3.3.4) for further in

precipitation caused by the pH modification is simulated by consid

line gibbsite [Al(OH)3]. The precipitation flag in record type 13 is 

the formation of precipitates.9

 Record types 14 and 15 are used to 

As shown in figure 28, record type 22 of the parameter file is 

2.3.4) so that the oxidation of ferrous iron is considered. Record ty

reduction, Fe(II) and Fe(III) (MINTEQ component numbers 280 an

tion of total dissolved iron concentration that is Fe(II) for each reac

that ferrous iron comprises 70–90 percent of total dissolved iron d

1992; Broshears and others, 1996).

Figure 28. Partial listing of the parameter file for Application 2.

Record types 27 and 28 of the parameter file are similar in for

are of note. First, four boundary conditions are specified, with the f

post-injection periods. Boundary conditions two and three corresp

pH. Boundary concentrations for Fe(II), Fe(III), Al, and SO4 are se

from samples collected above the injection. Because total excess h

Morgan, 1972), the boundary concentration for H must be determin

total inorganic carbon (CO3) is not directly measured and the boun

stand-alone MINTEQ run, fixed pH and a fixed partial pressure of 

lieu of total component concentrations for H and CO3. These fixed 

and total inorganic carbon (CO3) by mass balance, as described be

Execution of the stand-alone MINTEQ run proceeds as follow

a MINTEQ input file as described by Allison and others (1991). Ed

and to fix pH at 3.6, the observed pH above the injection site. Edit L

components with total concentrations equal to the boundary conce

to specify the fixed partial pressure of CO2 gas, so that total inorg

9

The precipitation flag for total excess hydrogen (H) is set to 1 because of the pr

Fe(OH)3 causes an increase in the excess hydrogen concentration (eq. 15, Section 

10The ionic strength specified should correspond to the ionic strength of the wat

from the complete suite of cations and anions. Ionic strength is usually determined 

observed solutes, wherein ionic strength is “computed” rather than “fixed.”

Record Types 22–24 from Applicat

Record Type

22

23

24

24

#######################################

#

# Parameters: Redox definition.

#

# Status: IREDOX required. Additional lines 

# reaction are required for IREDOX = 

#

#IREDOX [0=no redox,1=based on fraction di

#IDRD IDRD2 [for IREDOX>0]

#THETA [for IREDOX>0, J=1,NREACH]

#

#######################################

1

280 281

0.80

0.80
spheric carbon dioxide.11 The partial pressure of CO2 is equal t

ating an input file with PRODEF, MINTEQ is executed and the 

are set using the total component concentrations from the equilib

H, the value obtained from the MINTEQ output file is used as 

tions. For CO3, the MINTEQ-derived value is used to set the p

one and four). CO3 boundary concentrations for boundary cond

plus the CO3 added by the Na2CO3 injection. Input and output 

are available as part of the software distribution (Section 5.2).

The steady flow file for this application is similar to that p

Application 1 and the User’s Guide (Sections 3.3.5) for further 

4.2.2 The MINTEQ Input File — Application 2

Creation of the MINTEQ input file using PROTEQ gener

and Section 3.4.3. Edit Level I is first used to set record types 1

4.2.1, ionic strength is fixed at a level corresponding to the com

aqueous components (record type 9) corresponding to the six s

(Al(OH)3)

13 are defined as possible solids (Edit Level II), with 

precipitate definition (record type 15, parameter file).

4.2.3 Simulation Results — Application 2

Simulation results for the St. Kevin pH modification are co

iron (Fe(II)+Fe(III)) and dissolved aluminum concentrations de

above 5.0 (results for pH are similar to those presented by Bros

dissolved iron, an initial simulation that does not consider the o

results. As shown by the dashed line in figure 29a, the simulation 

This overestimation occurs because the bulk of the remaining d

Fe(OH)3 solid phase. Final simulation results (solid line, fig. 29

Fe(OH)3 precipitates, improving the correspondence between s

(fig. 29b), simulated and observed concentrations are in close a

agreement is immediately following the termination of the Na2C

~3.7). When the pH drops, precipitated mass residing on the stre

aluminum concentrations. Within the equilibrium-based model

observed data. This overestimation and the overestimation of th

on the redissolution of precipitated mass that is not considered 

11In many applications, total inorganic carbon is determined from alkalinity 

and CO3 must be set based on equilibrium with atmospheric CO2. An example 

12At sea level, air is 0.03 percent CO2, such that CO2 partial pressure equals 

spheric pressure is approximately 69 percent of the atmospheric pressure at sea 

atm. × 0.69).

13Another potential solid phase for the precipitation of Al(OH)3 is amorphou

fit to observed data.
60 One-Dimensional Transport with Equilibrium Chemistry (OTE

Figure 29. Simulated and observed concentrations of (a) total [Fe(II)

Time, 

Dissolved Concentration, in micromoles per liter

8 10 12 

St. Kevin Gulch p

x = 24 

b. Aluminum

a. Iron [Fe(II)+Fe(III)]

30

15

10

5

20 

100

200

 0
4.3 Application 3: Time-Variable Simulation of 

The previous two applications illustrate the use of OTEQ to 

example, OTEQ is used to simulate the sorption of solutes onto 

passage of a solute pulse. Other concepts introduced in this exa

flow file, and the determination of inorganic carbon from alkal

that may arise when modeling advection-dominated systems.

Copper sulfate (CuSO4) is often used to control algal bloom

1981, for example). In this hypothetical example, a 1,000-meter 

joins the main stem. The tributary is fed by a small lake that has 

concentrations of Cu and SO4 at the tributary outlet during a 15

main stem by sorption onto streambed sediments. The main ste

tact with the streambed and sorption/desorption reactions are in 

gible due to a lack of waterborne solid phases. Copper and sulf

Iron oxides are present in the streambed sediments, such that th

(Sections 2.3.3, 3.3.7, and 3.4.3).

As in Application 2, the number of simulated solutes is re

problem at hand. Main stem and tributary waters are circumneu

concentrations (with the exception of copper) are negligible. Si

carbon (CO3), copper (Cu), sulfate (SO4), calcium (Ca), and m

100.0 hours, with elevated concentrations of Cu and SO4 enteri

4.3.1 The Parameter File — Application 3

Record types 1–9 are used to specify various simulation pa

The tributary input to the main stem is modeled using the segme

three main reaches plus an additional dummy reach for consider

type 10). Tributary input is represented by reach 2. Each reach i

record type 11). Six solutes are defined using record types 12 a

solutes with defined reactions in the Dzombak and Morel (1990

(NPRECIPS=0, record type 12).

Sorption to the streambed is modeled using the static surfa

face is specified using record types 16–20 (fig. 31). Record typ

1) and sorption/desorption reactions are in a state of chemical eq

centration of the static surface (SOLCON) is specified for each 

and 20 are used to specify the initial sorbed concentrations for e

to specify the MINTEQ component number for SO4, the first s

type 13). Initial sorbed concentrations (SORBB) for SO4 are the

types (record type 19 and record type 20 used NREACH times) 

Cu, and Ca). The initial sorbed concentrations may be determin

Specification of the remaining record types in the paramet

2. Two differences are of note. First, due to the larger segment 

at segment boundaries (fig. 30), the interpolation option is invok

upstream boundary concentrations (USBC, record type 28) for 

on a stand-alone MINTEQ run (Sections 3.4.2 and 4.2.1). Exec

cation 2, except that CO3 is determined from alkalinity rather t

concentration is entered when using PRODEF (menu item 6, E

14For the case considered here, a prior OTEQ simulation was conducted with 

els. Sorbed concentrations at the completion of the simulation were then used to 

OTEQ simulation are available as part of the software distribution (Section 5.2
62 One-Dimensional Transport with Equilibrium Chemistry (OTE

Figure 30. Segmentation scheme for the hypothetical stream with tr

2

1

3

4

Reaches and Segments

Lake-fed Tributary

pH=7.0

Alkalinity=70 mg/liter CaCO3

Q=0.1 meter3 per second

Upstream Boundary

pH=8.2

Alkalinity=90 mg/liter 

Q=0.5 meter3 per seco

1

Reach Boundary

Segment Boundary

Reach Number
Figure 31. Partial listing of the parameter file for Application 3.

Record Types 16–20 from App

Record Type

16

17

18

18

##################################

#

# Parameters: Sorption definition.

#

# Status: ISORB required. Additional l

#

##################################

1 ISORB [0=no sorption, 1=stat

999.0 LAMBS

0.05 SOLCON (J=1,NREACH)

0.05

0.05

0.05

#

# SORBB, initial sorbed concentration on 

#

732

1.0265e-7 0.0

2.8759e-7 0.0

3.7581e-7 0.0

3.7581e-7 0.0

#

# SORBB, initial sorbed concentration on 

#

330

-2.0803e-5 0.0

-1.0875e-5 0.0

-9.8279e-6 0.0

-9.8279e-6 0.0

#

# SORBB, initial sorbed concentration on 

#

231

4.7210e-7 0.0

2.0026e-6 0.0

2.5169e-6 0.0

2.5169e-6 0.0

#

# SORBB, initial sorbed concentration on 

#

150

1.0188e-5 0.0

4.1481e-6 0.0

3.2769e-6 0.0

3.2769e-6 0.0

18

18

19

20

20

20

20

19

20

20

20

20

19

20

20

20

20

19

20

20

20

20
64 One-Dimensional Transport with Equilibrium Chemistry (OTE

4.3.2 The Unsteady Flow File — Application 3

The primary purpose of the unsteady flow file is to allow for 

unsteady. In this example, a steady-flow regime is considered (flo

constant), but the unsteady flow file is needed to accommodate the 

flow file for this purpose is described below; consideration of uns

and others (1998).

A partial listing of the unsteady flow file is shown in figure 3

the flow variables change. In this example, QSTEP is set to 5 hou

represented by three QSTEPs. Four flow locations are defined usin

ing (downstream) order, with the first and last locations at the upstr

second and third flow locations bracket the tributary, so that the ch

utary outlet are considered.

After defining the flow locations, record types 4–6 are used to 

nel cross-sectional area for each flow location. Record type 7 speci

location. The block of record types (record types 4–6 and record t

hours of simulation time. The first block represents hours 0–5, wh

levels. Blocks 2–4 represent hours 5–20 when the tributary concent

ing blocks (20–100 hours) are identical to the first block and repre

levels.

Lateral inflow concentrations (record type 7) at the third flow 

input. Inflow concentrations for Cu, SO4, Ca, and Mg are set to ob

collected at the tributary outlet. Inflow concentrations for H and CO

identical to that for the upstream boundary condition (Section 4.3.

4.3.3 The MINTEQ Input File — Application 3

Creation of the MINTEQ input file using PROTEQ generally 

Section 3.4.3. Edit Level I is first used to set record types 1–5, and 

sponding to the six solutes. The relevant sorption information is the

3, the diffuse layer model is selected from the available sorption al

face is added using menu item 1 of the sorption option menu. The s

best estimate of Dzombak and Morel (1990) (see Section 3.4.3 an

concentration of 5.618×10−5 moles per gram. The surface definitio

sorption option menu) with a concentration of 2.247×10−3 moles p

are based on the best estimates of site density and molecular weig

and table 26). The final step in setting up the sorption problem is t

ponents and the static surface. Sorption reactions for the static surfa

This option allows the sorption database (feo-dlm.dbs) to be adde

tions for the static surface.
Figure 32. Unsteady flow file for Application 3.

Unsteady Flow File 

Record Type

1

4

5

6

7

(Record types 4–7 repeat for each 

####################################

# OTEQ unsteady flow file

####################################

5.0 QSTEP [hour] 

####################################

# Flow Locations

####################################

 4 NFLOW

 0.0 (FLOWLOC for I = 1, NFL

 500.0

 530.0

1100.0

####################################

#

# repeating block for each QSTEP:

#

# line 1 - record type 4, QLATIN at each f

# line 2 - record type 5, Q at each flow loc

# line 3 - record type 6, AREA at each flow 

# lines 4-9 - record type 7, CLATIN for SO

#

#Flow Loc:

# 1 2 3 4

####################################

# 0-5 hours (background conditions)

 0.00 0.00 3.33e-3 0.00

 0.50 0.50 0.60 0.60

 2.00 2.00 2.40 2.40

 0.00 0.00 6.81e-5 0.00

 0.00 0.00 2.09e-3 0.00

 0.00 0.00 1.00e-7 0.00

 0.00 0.00 1.74e-3 0.00

 0.00 0.00 5.50e-4 0.00

 0.00 0.00 7.80e-5 0.00

# 5-10 hours (elevated Cu and SO4)

 0.00 0.00 3.33e-3 0.00

 0.50 0.50 0.60 0.60

 2.00 2.00 2.40 2.40

 0.00 0.00 7.00e-5 0.00

 0.00 0.00 2.08e-3 0.00

 0.00 0.00 2.00e-6 0.00

 0.00 0.00 1.74e-3 0.00

 0.00 0.00 5.50e-4 0.00

 0.00 0.00 7.80e-5 0.00

# 10-15 hours (elevated Cu and SO4)

 0.00 0.00 3.33e-3 0.00

 0.50 0.50 0.60 0.60

 2.00 2.00 2.40 2.40

 0.00 0.00 7.00e-5 0.00

 0.00 0.00 2.08e-3 0.00

 0.00 0.00 2.00e-6 0.00

 0.00 0.00 1.74e-3 0.00

 0.00 0.00 5.50e-4 0.00

 0.00 0.00 7.80e-5 0.00

2

3

3

3

3

7

7

7

7

7

4

5

6

7

7

7

7

7

7

4

5

6

7

7

7

7

7

7
66 One-Dimensional Transport with Equilibrium Chemistry (OTE

4.3.4 Simulation Results — Application 3

Simulation results for the hypothetical sorption problem are sh

70 meters downstream from the tributary outlet (600 meters) and a 

simulation (squares, fig. 33) shows the copper concentrations that 

sorption is modeled, arrival of the copper pulse at 600 meters is de

caused by sorption to streambed sediments located between the tri

meters). Copper concentrations at 600 meters begin to increase whe

the elevated tributary input, such that the upstream segments no lo

upstream segments is especially pronounced at 1,000 meters, wher

33). This extended delay results in a lower peak concentration and 

meters.

Figure 33. Copper concentrations at 600 and 1,000 meters, resulting 

4.3.5 Numerical Issues — Application 3

As described in Section 2.4, the governing transport equation

employed by the OTIS solute transport model (Runkel, 1998). The

which the finite-difference grid is fixed in space (Runkel and Chap

rian finite-difference methods are known to produce numerical os

advection-dominated systems are considered (Gray and Pinder, 19

low flow are not advection-dominated, due to the mixing induced b

of the numerical oscillations that arise when modeling advection-d

Time, i

Copper Concentration, in micromoles per liter

Simulated Copper C

0.1

0.2

20 40 

0

0

0.3

0.4

r

re

c
Figure 34a shows the spatial profile of copper concentratio

corresponds to the simulation described above (Section 4.3.4), w

and 10-meter segment lengths are used. Copper concentrations i

concentrations of the tributary relative to reach 1. Figure 34b il

increasingly dominated by advection (as D is decreased). A dis

the copper concentration at 495 meters that precedes the conce

persion coefficient to 0.5 meter2 per second causes oscillating lo

OTEQ may be inappropriate for some highly advective systems

grid. For the case considered here, the oscillations produced wh

eliminated by reducing the segment length from 10 to 2.5 mete

Figure 34. Spatial profiles of simulated copper concentration at 4 

simulations showing oscillations that arise when D is decreased.

Dis

Copper Concentration, in micromoles per liter

Effect of Advection

0.010

0.015

200 4

0.005

0

0.020

0.025

D = 0.5 meter2 p

D = 1.0 meter2 p

D = 2.0 meters2 p

a. Spatial Cu Profile

b. Spatial Cu Profile

0.006

0.009

0.012

0.015

450 460 

 (D = 2 meters2 
68 One-Dimensional Transport with Equilibrium Chemistry (OTE

4.4 Application 4: Steady-State Simulation of Exi

In addition to the time-variable simulations shown in the prev

state solute concentrations that result from a constant loading scen

steady-state solution of the governing transport equations (Section 

tions must therefore be based on time-variable runs of sufficient du

“spin-up,” Thornton and Rosenbloom, 2005; Section 2.4.3). An ad

calibrated model to simulate potential remedial actions.

Kimball and others (1991) describe a synoptic study conducte

the Rocky Mountains of Colorado. Saint Kevin Gulch receives acid

the toe of a large mine dump. Instream metal concentrations incre

downstream, a circumneutral tributary (pH 6.5) known as Shingle 

decrease downstream of the confluence with Shingle Mill Gulch as 

study, water samples were collected from the acidic springs (363–

inflows, and numerous instream sites. Flow rates for the instream 

transport parameters (D, A, AS, and α) were estimated using tracer 

Data from the study provide detailed spatial profiles of flow a

application, OTEQ is used in steady-state mode to evaluate the eff

is a two-step process. In the first step, the hydrologic and geochemi

conditions are quantified. A model of existing conditions is develo

(Broshears and others, 1993) and the inflow chemistry provided b

observed spatial profiles of pH, iron, and aluminum. In the second 

reflect the proposed remedial action. In the example considered he

through the addition of CaCO3. This remedial action is evaluated b

centrations associated with the springs are modified accordingly.

Simulated solutes in this application include total excess hyd

ferric iron (Fe(III)), aluminum (Al), and sulfate (SO4). Precipitatio

modeled using ferrihydrite (Fe(OH)3) and microcrystalline gibbsit

set of solutes is intentionally limited to those listed above for illustr

include a more complete set of solutes so that unforeseen solute in

4.4.1 Quasi-Steady-State Simulations

To model steady-state conditions, a single boundary condition 

upstream boundary concentrations (USBC, record type 28, parame

meters. In addition, a steady flow file is developed using the estim

types 2 and 3, steady flow file). Given a constant loading scenario 

inflow concentrations, and a steady-flow regime), simulated solute 

model time steps. Once steady state is achieved, solute concentrat

profiles corresponding to the steady-state solution. These spatial co

equal to 1 (DOPT, record type 25, parameter file) so that concentra

end of the model run.

The time to reach steady state is application dependent. Becau

a small number of time steps may be needed to reach steady state. 

ever, when settling of waterborne solids (precipitated or sorbed ma

ered. Model users must therefore inspect the solute output files (Se

steady state. If concentrations have not reached steady state, solute 

not correspond to steady-state conditions. For the case considered 

subject to settling (PSETTLE>0, record type 15, parameter file). S

results at the most downstream location of interest (1,804 meters)

1,804 meters quickly reach steady-state levels. In contrast, 50 hour

waterborne Fe(III). OTEQ simulations of steady-state conditions a

tion time (TSTART=0, TFINAL=80, record types 4 and 5, parame
Figure 35. Time required to reach quasi-steady-state conditions. 

4.4.2 Modeling Existing Conditions and Remed

 The parameter, flow, and MINTEQ input files are similar 

tions. The system is modeled using seven reaches. Reach 1 is a 

record type 3, steady flow file). Flow in reach 2 (26–363 meter

during the synoptic. Inflow concentrations (CLATIN, record ty

that the diffuse inflows are chemically similar to the stream. Rea

at the toe of the mine dump. Direct use of observed data from th

centrations at 500 meters. Lateral inflow concentrations are the

need for this upward adjustment suggests the presence of a gro

springs. Reach 4 (500–510) is a tributary reach (as in Applicati

inflow concentrations for reach 4 are set using observed data fro

liter as CaCO3. Reaches 5 and 6 (510–1,557 meters) introduce 

observed data for a sampled inflow at 1,281 meters. Reach 7 (1,

(QLATOUT>0, record type 3, steady flow file). To minimize er

2.4.2 and 3.3.4), reach 7 extends 100 meters downstream of the 

applications, lateral inflow concentrations for H and CO3 in each 

and 4.2.1).

Use of the reach chemistry described above results in a sim

and aluminum. This calibrated model of existing conditions can 

various remedial actions. A hypothetical remediation plan in w

model this remedial action, 5.6×10−3 moles per liter of CaCO3 

A stand-alone MINTEQ run indicates that addition of CaCO3 r

tem. Reach 3 inflow chemistry (CLATIN, record type 3, steady 

dissolved component concentrations from the equilibrated mas

used to simulate the effects of remediation. Input and output fil

are available as part of the software distribution (Section 5.2).

T

Fe(III) Concentration, in micromoles per liter

Time to Quas

60

0 20 

70

80

50

40

30
70 One-Dimensional Transport with Equilibrium Chemistry (OTE

4.4.3 Simulation Results — Application 4

Simulation results for existing conditions are shown in figure 

meters as the acidic springs enter the main stem (fig. 36a–c). Shing

a decrease in metal concentrations and an increase in pH. The proc

tions are quantified through model application. For iron, the pH in

precipitation of ferrihydrite (fig. 36b). The decrease in iron at 500 

total waterborne iron) and chemical precipitation (the difference be

decrease in concentration is due solely to dilution (fig. 36c; total w

ing an absence of precipitated gibbsite).

The potential effects of remediation on pH, iron, and aluminu

tions, the inflow entering reach 3 (363–500 meters) has a higher pH 

4

 versus 3.15×10−3 moles per liter) due to treatment of the acidic s

concentrations decrease (fig. 36e), relative to existing conditions, 

causes further increases in pH and decreases in iron concentration

iron concentration downstream of Shingle Mill Gulch are attribute

trations in the reach 3 inflow waters are unaffected by treatment (al

concentrations in the vicinity of the springs are therefore unaffecte

affect aluminum concentrations further downstream, as gibbsite pr

Gulch where pH exceeds 4.8. In summary, treatment results in the 

nitude of the removal differ for the two solutes. Iron reductions ar

system), whereas aluminum reductions are due to instream remov

4.5 Application 5: Steady-State Simulation of Sor

In Application 3, OTEQ is used to simulate the time-variable 

orption that occurs following the passage of a solute pulse. As in A

steady-state solute concentrations that result from a constant loadin

the streambed are saturated such that sorption to the streambed is n

lying water column). Reactions in the water column, however, may 

acts as a sorptive surface. The purpose of this application is theref

tates as detailed in Section 2.3.3. In addition, this application illust

gen and total inorganic carbon for the upstream boundary and late

Runkel and Kimball (2002) describe a synoptic study conduc

near Silverton, Colorado. Mineral Creek receives acidic and metal

ground-water discharge zones located along the water course, resu

from the synoptic study provide detailed spatial profiles of flow an

application, OTEQ is used in steady-state mode to model a subsec

Kimball (2002). The upstream boundary of the subsection is located 

in the formation of waterborne solid phases that are transported do

the high velocity of Mineral Creek and the small size of the solid m

boundary has important implications for solute transport, as precip

data at the upstream boundary suggests the presence of precipitate

precipitated phase. Additional sorption and desorption of metals co

ported into reaches with different pH and(or) metal sources (for ex

Runkel and Kimball, 2002). These sorption reactions are modeled 

2.3.3, 3.3.7, and 3.4.3) as in Application 3.

Simulated solutes in this application include total excess hyd

ferric iron (Fe(III)), aluminum (Al), sulfate (SO4), and arsenic (As

(Fe(OH)3) and microcrystalline gibbsite (Al(OH)3) as possible sol

limited to those listed above for illustrative purposes. Simulations 

Runkel and Kimball (2002). As in Application 4, simulations of st

cient duration to obtain a quasi-steady state (Section 4.4.1).

15Dissolved and sorbed arsenic species in Mineral Creek are assumed to be in th

(As(V)) is the dominant form of arsenic in oxygenated surface waters (Hem, 1985
Figure 36. Spatial profiles: (a) pH, existing conditions; (b) iron, exi

tion; (e) iron, remediation; (f) aluminum, remediation.

3

4

5

0

100

200

0 500 1,000 1,500 

25

75

125

St. Kevin Gulch, A

Dis

Existing Conditions

Iron [Fe(II)+Fe(III)] 

Concentration, in micromoles per liter

Aluminum (Al) Concentration, 

in micromoles per liter pH, in standard units

a. 

b. 

c. 

observe

simulat

observed dissolve

simulated 

dissolved (C) 

simulated

total waterborne 

observed dissolved

simulated 

dissolved (C) 

(=total waterborne, 
72 One-Dimensional Transport with Equilibrium Chemistry (OTE

4.5.1 The Parameter, Flow, and MINTEQ Input File

 The parameter, flow, and MINTEQ input files are similar in 

tions. Within the parameter file, several differences are of note. Fi

subsection of the Runkel and Kimball (2002) study reach that is co

solid-phase settling in Mineral Creek is modeled with settling velo

settling, solute concentrations reach steady-state values over a rela

−TSTART, record types 4 and 5) required to reach quasi-steady st

4.4.1). Third, sorption onto waterborne solid phases is modeled usin

flag in record type 13 is first set to 1 for the three solutes with defi

SO4, and As(V)). Record type 16 is then used to invoke the dynam

is used to specify the MINTEQ component number of the solute tha

weight of the sorbent (MWSORB). As shown in figure 37, IDSOR

and MWSORB is set to the value provided by Dzombak and More

to precipitated Fe(III) that is present in the water column (ferrihyd

are determined from a series of stand-alone MINTEQ runs as desc

Figure 37. Partial listing of the parameter file for Application 5.

The steady flow file for this application is generally similar t

of flow at the upstream boundary and the cross-sectional area of e

for Fe(II), Fe(III), Al, SO4, and As(V) are set to observed total wat

concentrations for H and CO3 are determined from a series of stan

Creation of the MINTEQ input file using PROTEQ generally 

Section 3.4.3. Edit Level I is used to set record types 1–5; Edit Le

ferrihydrite (Fe(OH)3) and microcrystalline gibbsite (Al(OH)3) as 

tion for the dynamic surface is identical to that presented for Appl

site type is defined with a concentration of 1.124×10−4 moles per g

of freshly precipitated iron oxides (Runkel and others, 1999). Follo

ually edited to change the default value of the surface complexatio

ation constant is made to improve the correspondence between sim

Kimball, 2002).

Record Types 16 and 21 from Applic

Record Type

16

####################################

#

# Parameters: Sorption definition.

#

# Status: ISORB required. Additional line

#

####################################

2 ISORB [0=no sorption, 1=static, 

#

# dynamic surface

#

# IDSORB MWSORB

# |

#------------------------------

21 281 89.0
4.5.2 Specification of H and CO3: The Case of W

In Applications 2–4, total excess hydrogen (H) and total in

and lateral inflows are set based on output from a single stand-al

is relatively straightforward due to the absence of waterborne s

equal to dissolved concentrations at the upstream boundary and 

significant waterborne solid concentrations present at the upstr

stances, a more detailed approach to setting H and CO3 is requi

reproduced by the model at the locations of interest (upstream b

trated below for the upstream boundary condition; the same ap

interest has waterborne solid phases.

Determination of H and CO3, Step 1

Determination of H and CO3 concentrations at the upstrea

MINTEQ run is used to determine H and CO3 concentrations fo

is used to create a MINTEQ input file. Edit Level I is used to sp

II is used to define Fe(II), Fe(III), Al, SO4, and As(V)16 as aqu

concentrations specified in record type 28 of the parameter file 

up to be calculated based on equilibrium with atmospheric CO2 

MINTEQ is executed and the concentrations of H and CO3 list

distribution) are recorded for use in Step 2.

Determination of H and CO3, Step 2

Step 2 involves a series of stand-alone MINTEQ runs in w

determine the total excess hydrogen concentration that yields th

As before, the MINTEQ input file is created using PRODEF, as 

strength. Edit Level II is used to define Fe(II), Fe(III), Al, SO4

define H and CO3 as aqueous components, with total concentra

completed by using Edit Level II to define the possible solids t

gibbsite). After creating an input file, MINTEQ is executed and 

MINTEQ runs are then conducted in which the total excess hydr

observed and estimated pH. Results from this trial and error pro

shown in table 27. As shown in the table, three stand-alone MI

upstream boundary (4.43). The H concentration used as input fo

are then used to set the upstream boundary concentrations with

16The MINTEQ component for As(V) is H3AsO4.

17For the case of the Mineral Creek upstream boundary, the acidic water has 

when alkalinity is present at the upstream boundary (or in the inflow of interest

observed alkalinity (Edit Level I, Section 4.3.1).

Table 27. Stand-alone MINTEQ runs 

MINTEQ run H concen

v

1 −1.019×1

2 −1.519×1

3 −2.019×1
74 One-Dimensional Transport with Equilibrium Chemistry (OTE

4.5.3 Simulation Results — Application 5

Simulation results for the Mineral Creek synoptic study are sh

boundary condition described in the previous section results in sim

boundary (888 meters). The simulated precipitation of ferrihydrite 

that are comparable to the observed data (fig. 38a and b). The pres

itation provides a dynamic surface onto which arsenic sorption occ

38c). Further downstream, simulated values are generally consiste

observed and simulated values for pH, dissolved iron, and dissolv

underestimate the observed data in reaches 3–7 (1194–1989 meter

high metal concentrations (Runkel and Kimball, 2002).

Figure 38. Spatial profiles of simulated and observed: (a) pH, (b) iron

7

300

1,000 1,500 2,00

4

Mineral Cr

Dista

Iron [Fe(II)+Fe(III)] 

 Concentration, in 

in micromoles per liter

pH, in standard units

a.

3

2

1

0

200

100

0

5

3

b.

c.

Arsenic Concentration, micromoles per liter
5 Software Guide

This section provides information on how to obtain and in

are described in Section 5.5.

5.1 Supported Platforms

Executable binary files are available for personal compute

the supported operating systems and hardware platforms is pro

compile the source code as described in Section 5.4.

Column 3 of table 28 contains ARC, the notation used thro

architecture. Personal computer users running Linux, for exam

the specific instructions given in Sections 5.2–5.4.

5.2 Software Distribution

OTEQ may be obtained over the Internet at http://water.us

hardware-specific executables, and example input files from the 

executable and the set of example input files. These files are sto

presented in table 29.

Table 28. Supported systems.

Hardware 

Personal computer (Intel x86) Linu

Unix workstation (SUN) Solar

Table 29. Files to download.

File contents 

Executable files 

Source code 

Database files, documentation, and example input files 
76 One-Dimensional Transport with Equilibrium Chemistry (OTE

5.3 Installation

The files described in Section 5.2 may be used to install the O

procedure consists of three steps: (1) creating the OTEQ directory s

areas. All three steps should be completed on the user’s target platfo

within Microsoft Windows is known to cause installation problem

5.3.1 Creating the OTEQ Directory Structure

The OTEQ directory structure is created using the downloade

the base-directory. The base-directory is a user-selected directory 

After moving the files, unzip the tar files and extract the directory 

(1) Unzip the tar files:

gunzip oteq.dde.tar.gz

gunzip oteq.src.tar.gz 

gunzip oteq.ARC.tar.gz 

(2) Extract the directory structure:

tar -xovf oteq.dde.tar

tar -xovf oteq.src.tar 

tar -xovf oteq.ARC.tar 

5.3.2 Updating the User’s Path

After completing the installation procedure, oteq, minteq, a

oteq1.40/bin (fig. 39). A path to these executables is created by 

will be running one of two types of shells: a C shell or a Bourne sh

the multitude of available shells and the complexity of the startup 

(for example, specific shells may use files other than .cshrc and .p

installation should contact their system administrator. Also note th

user opens a new window or the startup file is executed using the 

C Shells

Model users running C shells such as csh and tcsh can updat

~/.cshrc file:

setenv OTEQ_HOME base-directory

set path=($OTEQ_HOME/oteq1.40/bin $path)

where base-directory is the user-selected directory where the OTE

Bourne Shells

Model users running Bourne shells such as sh, ksh, and bash 

of their ~/.profile file:

OTEQ_HOME=base-directory

export OTEQ_HOME

PATH=$OTEQ_HOME/oteq1.40/bin:$PATH

export PATH

where base-directory is the user-selected directory where the OTE

18Where ARC equals LIN or SOL (see tables 28 and 29).
Figure 39. OTEQ directory structure.

5.3.3 Creating User Work Areas

User work areas are user-defined directories from which OT

will complete multiple model runs, using multiple directories to 

ated over time.

Creation of user work areas is illustrated using a simple ex

only; model users may name their directories as they see fit. To 

the OTEQ and MINTEQ work areas. For example, one might cr

of the Red River. As discussed in Section 3.4.2, one or more stan

A subdirectory and work area for the stand-alone MINTEQ run

cd redriver

mkdir minteq

cd minteq

mkminteqlinks

where the final command listed above (mkminteqlinks) create

to the actual database files stored in base-directory/dbs (fig. 39

of an OTEQ work area proceeds in a similar manner, by issuing 

base-d

OTEQ Directo

oteq

bin 

oteq

EXPLANATION

executable

directory

dbs d

proteq

posteq

minteq

prodef
78 One-Dimensional Transport with Equilibrium Chemistry (OTE

mkdir run

cd run

mkoteqlinks

where the final command listed above (mkoteqlinks) creates a se

execution.20 If the user wishes to conduct multiple OTEQ runs, the 

for other directories (run2, run3, etc., for example).

5.4 Compilation

As discussed in Section 5.1, executable binary files for the OT

platforms and operating systems. As a result, compilation of the so

following situations:

• Executable binaries are not available for your specific har

• The dimensions of the modeled system exceed the maxim

• Modifications have been made to the underlying compute

The model is compiled using the make utility. To compile OTEQ

make FFLAGS=" " -f makefile.oteq

on the command line while in base-directory/oteq1.40/src (inser

To create a specific executable, append ../bin/oteq, ../bin/prote

MINTEQ and PRODEF, type:

make -f makefile.minteq

on the command line while in base-directory/oteq1.40/src. To cr

prodef to the above command. Note that compilation flags for cod

and PRODEF, due to issues with the MINVAL subroutine. When c

base-directory/oteq1.40/bin.

5.5 Software Overview

5.5.1 Model Development

The OTEQ solute transport model is written in ANSI standar

hardware platforms (Section 5.1). Specific development environm

19Users may later remove these links using the rmminteqlinks command. Use of 

not take up disk space).

20Users may later remove these links using the rmoteqlinks command. Use of rm

take up disk space).

Table 30. Development environments.

Processor Operating system 

Dual Core Xeon 5050 Mandriva 2007 Portland Gr

SUN SPARC 10 Solaris 2.6 SPARComp
5.5.2 Include Files

The source code used to develop OTEQ consists of many s

files are used. Use of include files allows program information 

fied by editing the include files rather than each individual rout

Maximum Dimensions — fmodules.inc

Under Fortran-77, the dimension of each vector and array 

knowledge of the maximum dimension for each model paramet

tant task because excessively small values limit program applic

To address this problem, the maximum dimensions for the 

Increases or decreases in the maximum dimensions are made b

Section 5.4. Default values for the maximum dimensions are gi

a user-supplied input variable. This correspondence is shown p

When running the model, the input variables may not excee

3.3.4 — Record type 25), for example, may not exceed the max

its given maximum, program execution is terminated and an er

increase the appropriate maximum value (by editing fmodules

Logical Devices — lda.inc

 In Fortran-77, a unit number is assigned to each file used f

device assignments (ldas), must be specified for each read and 

are shared between the input and output subroutines using a For

file lda.inc.

Table 31. Maximum dimensions and default values from fm

Dimension Default maximum 

MAXADS 2 Ad

MAXBOUND 200 Up

MAXCOMP 26 Co

MAXFLOWLOC 30 Fl

MAXPRECIP 10 Pr

MAXPRINT 50 Pr

MAXREACH 70 St

MAXSEG 8000 St

MAXSOLUTE 16 So

MAXSPEC 200 Sp
80 One-Dimensional Transport with Equilibrium Chemistry (OTE

5.5.3 Error Checking

OTEQ’s input subroutines perform several tests to validate the 

ten to echo.out and program execution is terminated. The error-ch

• The number of reaches, NREACH, must not exceed the m

• The number of segments, IMAX, must not exceed the ma

• The number of print locations, NPRINT, must not exceed 

• The number of solutes, NSOLUTE, must not exceed the m

• The number of upstream boundary conditions, NBOUND

• The number of flow locations, NFLOW, must not exceed 

• The number of precipitates, NPRECIPS, must not exceed 

• Sorption information for solutes sorbing to the static surfa

solute definition.

• The number of precipitates specified in the parameter file, 

specified in the MINTEQ input file.

• A given print location, PRTLOCi

, must lie within the mod

• The flow locations, FLOWLOC, must be entered in ascen

• The first flow location, FLOWLOC1, must be placed at the 

the downstream boundary.

• The kinetic rate coefficient for sorption, LAMBS, must be 

of the integration time step (in seconds).

• The chemistry option, ICHEM, must equal 1 or 2.

• Solutes specified in the parameter file must appear as TYP

• Precipitates specified in the parameter file must be define

• The sorption option must be 0, 1, 2, or 3.

• The component Id specified for the dynamic surface, IDS

• The interpolation option, IOPT, must equal 0 or 1.

• The distance option, DOPT, must equal 0 or 1.

• The print option, PRTOPT, must equal 1 or 2.

• The boundary condition option, IBOUND, must equal 1, 

• The redox option, IREDOX, must equal 0 or 1.

• The component Ids specified for redox, IDRED and IDRE

• The sorption option in the MINTEQ input file, IADS, mu

• The storage zone cross-sectional area, AREA2, must be g

• The coefficient of the depth-area power function, A1, mu
References Cited

Alley, W.M., and Smith, P.E., 1982, Distributed routing rainfal

user’s manual: U.S. Geological Survey Open-File Report 82

Allison, J.D., Brown, D.S., and Novo-Gradac, K.J., 1991, MIN

ronmental systems —Version 3.0 user’s manual: U.S. Enviro

Ball, J.W., and Nordstrom, D.K., 1991, WATEQ4F Users manu

lating speciation of major, trace and redox elements in natura

Bencala, K.E., 1983, Simulation of solute transport in a mounta

sorption: Water Resources Research, v. 19, no. 3, p. 732–73

Bencala, K.E., and Walters, R.A., 1983, Simulation of solute tra

model: Water Resources Research, v. 19, no. 3, p. 718–724.

Bird, R.B., Stewart, W.E., and Lightfoot, W.E., 1960, Transpor

Broshears, R.E., Bencala, K.E., Kimball, B.A., and McKnight, 

ulations for a mountain stream, Saint Kevin Gulch, Colorado

92–4081, 18 p.

Broshears, R.E., Runkel, R.L., Kimball, B.A., McKnight, D.M

stream — Experimental pH increase and simulation of contr

nology, v. 30, no. 10, p. 3016–3024.

Brown, K.P., and Hosseinipour, E.Z., 1991, New methods for m

into surface waters: Journal of Environmental Science and H

Chapman, B.M., 1982, Numerical simulation of the transport an

Resources Research, v. 18, no. 1, p. 155–167.

Chapra, S.C., 1997, Surface water-quality modeling: New York

Chapra, S.C., and Canale, R.P., 1988, Numerical methods for e

Chen, C.W., Leva, D., and Olivieri, A., 1996, Modeling the fat

mental Engineering, v. 122, no. 10, p. 924–934.

Di Toro, D.M., 1976, Combining chemical equilibrium and phy

chemical Processes in Aquatic Ecosystems, Canale, R.P., ed

Dzombak, D.A., 1986, Toward a uniform model for the sorptio

Technology, unpublished Ph.D. dissertation, 521 p.

Dzombak, D.A., and Morel, F.M.M., 1990, Surface complexati

p.

Gray, W.G., and Pinder, G.F., 1976, An analysis of the numeri

v. 12, no. 3, p. 547–555.

Harrison, David, 1989, Xgraph Version 11.3.2: University of C

Hem, J.D., 1985, Study and interpretation of the chemical char

Supply Paper 2254, 264 p.

HydroGeoLogic, Inc., and Allison Geoscience Consultants, Inc

model for environmental systems — user manual supplemen

HydroGeoLogic, Inc., 1999, Diffuse-layer sorption reactions fo

Virginia, 15 p.

Jobson, H.E., 1989, Users manual for an open-channel streamfl

Water-Resources Investigations Report 89–4133, 73 p. 

Keefe, S. H., Barber, L.B., Runkel, R.L., Ryan, J.N., McKnigh

transport in constructed wetlands: Water Resources Research

Kimball, B.A., McKnight, D.M., Wetherbee, G.A., and Harnish

acidic stream (St. Kevin Gulch, Colorado, USA): Chemical G

Kimball, B.A., Broshears, R.E., Bencala, K.E., and McKnight, 

cesses in a stream affected by acid mine drainage, in U.S. Geo

ings of the technical meeting, Monterey, California, March 1

ical Survey Water-Resources Investigations Report 91–4034

Kuwabara, J.S., Leland, H.V, and Bencala, K.E., 1984, Copper 

Engineering, v. 110, no. 3, p. 646–655.
82 One-Dimensional Transport with Equilibrium Chemistry (OTE

McKnight, D.M., 1981, Chemical and biological processes control

A field study of the CuSO4 treatment of Mill Pond Reservoir, Bu

no. 3, p. 518–531.

McKnight, D.M., Kimball, B.A., and Bencala, K.E., 1988, Iron ph

ence, v. 240, p. 637–640.

McKnight, D.M., Kimball, B.A., and Runkel, R.L., 2001, pH depe

affected by acid mine drainage: Hydrological Processes, v. 15, p

Morel, F.M.M., and Hering, J.G., 1993, Principles and application

Morel, F., and Morgan, J., 1972, A numerical method for computi

ence and Technology, v. 6, no. 1, p. 58–67.

Nordstrom, D.K., Plummer, L.N., Langmuir, Donald, Busenberg, E

Revised chemical equilibrium data for major water-mineral reac

eds., Chemical modeling in aqueous systems II: Washington D.

31, p. 398–413.

Nordstrom, D.K., and May, H.M., 1996, Aqueous equilibrium dat

environmental chemistry of aluminum: Boca Raton Fla., CRC P

Nordstrom, D.K., and Archer, D.G., 2002, Arsenic thermodynamic 

lenwerk, K.G., eds., Arsenic in ground water: Boston, Kluwer A

Parkhurst, D.L., and Appelo, C.A.J., 1999, User’s guide to PHREE

reaction, one-dimensional transport, and inverse geochemical ca

tions Report 99–4259, 312 p.

Rubin, Jacob, 1983, Transport of reacting solutes in porous media 

tion and chemical nature of reactions: Water Resources Researc

Runkel, R.L., 1993, Development and application of an equilibrium

streams: Boulder, University of Colorado, unpublished Ph.D. di

Runkel, R.L., 1998, One-dimensional transport with inflow and sto

U.S. Geological Survey Water-Resources Investigations Report 

Runkel, R.L., 2000, Using OTIS to model solute transport in strea

Runkel, R.L., Bencala, K.E., Broshears, R.E., and Chapra, S.C., 19

equilibrium-based model: Water Resources Research, v. 32, no

Runkel, R.L., and Chapra, S.C., 1993, An efficient numerical soluti

streams: Water Resources Research, v. 29, no. 1, p. 211–215.

Runkel, R.L., and Chapra, S.C., 1994, Reply to comment on “An ef

solute transport in small streams”: Water Resources Research, v

Runkel, R.L., and Kimball, B.A., 2002, Evaluating remedial alterna

tive transport model: Environmental Science and Technology, v

Runkel, R.L., Kimball, B.A., McKnight, D.M., and Bencala, K.E.

plexation approach for trace metal sorption: Water Resources R

Runkel, R.L., Kimball, B.A., Walton-Day, K., and Verplanck, P.L

water quality—Red Mountain Creek, Colorado: Applied Geoch

Runkel, R.L., McKnight, D.M., and Andrews, E.D., 1998, Analysi

in an Antarctic stream: Journal of the North American Bentholo

Runkel, R.L., McKnight, D.M., Bencala, K.E., and Chapra, S.C., 

pH modification experiment: Water Resources Research, v. 32, 

Thomann, R.V., and Mueller, J.A., 1987, Principles of surface wate

p.

Thornton, P.E., and Rosenbloom, N.A., 2005, Ecosystem model sp

trial carbon and nitrogen cycle model: Ecological Modelling, v. 

Westall, J.C., Zachary, J.L., and Morel, F.M.M., 1976, MINEQL — 

rium composition in aqueous systems: Massachusetts Institute o

Yeh, G.T., and Tripathi, V.S., 1989, A critical evaluation of recent d

multichemical components: Water Resources Research, v. 25, n

Yeh, G.T., and Tripathi, V.S., 1991, A model for simulating transp

and demonstration: Water Resources Research, v. 27, no. 12, p. 
Glossary

Summary of variables used in Section 2. The fundamental units 

used herein.

A main channel cross-sectional area [L2]

AS storage zone cross-sectional area [L2]

C concentration of a component’s dissolved phase [mole

Cs storage zone concentration of a component’s dissolved 

Ctot sum of all dissolved species for the two components in

main channel concentration of an arbitrary phase [mol

fixed concentration at the upstream boundary [moles p

lateral inflow concentration of the arbitrary phase [mo

storage zone concentration of an arbitrary phase [mole

CC exp(−ΨF/RTa)

D dispersion coefficient [L2T–1]

F Faraday constant [96,485 coulomb mole−1]

H+ hydrogen ion

K equilibrium constant

Kint intrinsic surface complexation constant

L() transport operator

M molecular weight of the sorbent [gram sorbent per mol

M number of complexed species

M2+ divalent cation

N number of stream segments

NS site density [moles of sites per mole sorbent]

P concentration of a component’s precipitate phase (Pw+

Pb concentration of a component’s immobile precipitate p

Ps storage zone concentration of a component’s precipita

[moles per liter]

Psb storage zone concentration of a component’s immobile 

[moles per liter]

Psw storage zone concentration of a component’s mobile p

Pw concentration of a component’s mobile precipitate pha

PMe total precipitate concentration for the specified hydrou

Pb

Me immobile precipitate concentration of a specified hydr

Pw

Me mobile precipitate concentration of a specified hydrou

Q volumetric flow rate [L3T–1]

R molar gas constant [8.314 joules mole−1 K−1

]

S concentration of a component’s sorbed phase (Sw+Sb) 

Sb concentration of a component’s immobile sorbed phas

Ss storage zone concentration of a component’s sorbed p

Ssb storage zone concentration of a component’s immobile 

C

 ̃

C

 ̃ bc

C

 ̃ L

C

 ̃ S
84 One-Dimensional Transport with Equilibrium Chemistry (OTE

Ssw storage zone concentra

Sw concentration of a com

SA specific surface area [m

SC concentration of the so

S1 sorbed concentration as

S2 sorbed concentration as

S3 sorbed concentration as

SOH uncharged surface hydr

SOM+ positively charged surf

T total component concen

Ta absolute temperature [K

Ts total storage zone comp

TH total component concen

Z valence of a symmetric

a coefficient of the depth

ai stoichiometric coefficie

am stoichiometric coefficie

b exponent of the depth-a

c concentration of the un

ce molar electrolyte conce

d1 effective settling depth 

d2 effective settling depth 

fb source/sink term for di

fbm source/sink term for dis

fsb source/sink term for di

[moles per liter T−1]

fw source/sink term for pr

[moles per liter T–1]

gb source/sink term for so

[moles per liter T–1]

gsb source/sink term for so

[moles per liter T−1]

gw source/sink term for so

i subscript denoting strea

k superscript denoting th

m superscript denoting th

n superscript denoting tim

np number of solid precipi

pbm immobile precipitate co

pm total precipitate concen

qLIN lateral inflow rate [L3T

sext source/sink term repres

ssext source/sink term repres
[moles per liter T–1]

t time [T]

v1 main channel settling velocity [LT–1]

v1m settling velocity for precipitated species m [LT−1]

v2 storage zone settling velocity [LT–1]

w channel width [L]

x distance [L]

xi concentration of the ith complexed species [moles per 

* superscript denoting the equilibrium concentration in t

α storage zone exchange coefficient [T–1]

Δt integration time step [T]

Δx stream segment length [L]

ε dielectric constant of water

εo permittivity of free space [8.876×10−12 coulomb volt−

σ net surface charge density [coulomb meter−2]

σe relative error tolerance for sequential iteration procedu

Γ fraction of the equilibrium quantity that is allowed to so

step

Ψ surface potential [volts]

θtarget fraction of the total dissolved concentration that is asso

involved in oxidation/reduction
86 One-Dimensional Transport with Equilibrium Chemistry (O

Appendix 1. Modifications to the MINTE

The database files distributed with OTEQ are based on the 

others, 1991). The version 3 files have been updated to provide 

files are generally consistent with the databases of WATEQ (B

1999), and the published data of Nordstrom and others (1990), 

Details specific to each of the updated database files are discusse

the best available information, OTEQ users should note that:

“logK’s and enthalpies of reaction have been taken from v

made to determine the aqueous model that was used to dev

defined by the current data base files are consistent with th

the program should be considered to be preliminary. Caref

left to the users of the program” (Parkhurst and Appelo, 1

Given this potential uncertainty in the logK values and enthalpi

tivity analysis techniques to their model applications. For a disc

effects of reaction enthalpy, see Runkel and others (2007).

Aqueous Species - thermo.dbs

Enthalpy values and equilibrium constants (logKs) contain

tained in wateq4f.dat, a database file distributed with PHREEQ

were identical, as shown in table 32. Enthalpy and(or) logK val

wateq4f.dat for 33 aqueous species; OTEQ’s version of therm

(table 33). In addition to providing for consistency with PHREE

closer to those provided in the revised database distributed with 

ous species not listed in tables 32 and 33 have not been checke

original values from version 3 of MINTEQ.

Mineral Species - type6.dbs

Enthalpy and logK values contained in MINTEQ’s type6.

database file distributed with PHREEQC (Parkhurst and Appelo

in table 34. Enthalpy and(or) logK values in MINTEQ’s type6.

species (table 35). OTEQ’s version of type6.dbs has not been c

logK values remain at the original values from version 3 of MI

(to values from wateq4f.dat or other sources) when adding poss

Enthalpy and logK values for mineral species not listed in table

unlisted species remain at the original values from version 3 of 

Sorption to HFO - feo-dlm.dbs and feo-dlm2.dbs

The auxiliary database file for sorption of inorganic ions o

of MINTEQ (Allison and others, 1991) is based on the doctora

tributed with OTEQ (feo-dlm.dbs and feo-dlm2.dbs, Section 

Dzombak and Morel (1990). The updated files include 40 reacti

MINTEQ (HydroGeoLogic Inc., 1999) and PHREEQC (Parkhu
Table 32. Aqueous species for which MINT

identical to those in wateq4f.dat, a database 

MINTEQ ID Species 

201300 AgBr 

201301 AgBr2- 

201302 AgBr3-2 

201800 AgCL 

201801 AgCL2- 

201802 AgCL3-2 

201803 AgCL4-3 

202700 AgF 

203300 AgOH 

203301 Ag(OH)2- 

203800 AgI 

203801 AgI2- 

203802 AgI3-2 

203803 AgI4-3 

204910 Ag(NO2)2- 

204920 AgNO3 

207300 AgHS 

207301 Ag(HS)2- 

207302 Ag(S4)2-3 

207303 AgS4S5-3 

207304 Ag(HS)S4-2 

207320 AgSO4- 

902700 BF(OH)3- 

1505800 CaHPO4 

1505801 CaPO4- 

1505802 CaH2PO4+ 

1601300 CdBr+ 

1601301 CdBr2 
88 One-Dimensional Transport with Equilibrium Chemistry (OTE

1601800 CdCl+ 

1601801 CdCl2 

1601802 CdCl3- 

1601803 CdOHCl 

1602700 CdF+ 

1602701 CdF2 

1603300 CdOH+ 

1603301 Cd(OH)2 

1603302 Cd(OH)3- 

1603303 Cd(OH)4-2 

1603304 Cd2OH+3 

1604920 CdNO3+ 

1607320 CdSO4 

1607321 Cd(SO4)2-2 

2301800 CuCl2- 

2301801 CuCl3-2 

2311400 CuCO3 

2311401 Cu(CO3)2-2 

2311800 CuCl+ 

2311801 CuCl2 

2311802 CuCl3- 

2312700 CuF+ 

2313300 CuOH+ 

2313301 Cu(OH)2 

2313302 Cu(OH)3- 

2313303 Cu(OH)4-2 

2313304 Cu2(OH)2+2 

2317320 CuSO4 

2803300 FeOH+ 

Table 32. Aqueous species for which MINTEQ 

identical to those in wateq4f.dat, a database dis

MINTEQ ID Species 

[
2803301 FeOH3-1 

2803302 FeOH2 

2807300 Fe(HS)2 

2807301 Fe(HS)3- 

2807320 FeSO4 

2811800 FeCl+2 

2811801 FeCl2+ 

2811802 FeCl3 

2812700 FeF+2 

2812701 FeF2+ 

2812702 FeF3 

2813300 FeOH+2 

2813304 Fe2(OH)2+4 

2813305 Fe3(OH)4+5 

3300603 H4AsO3+ 

3300900 H2BO3-1 

3302702 H2F2 

3304900 NH3 

3305800 HPO4-2 

3305801 H2PO4- 

3307300 H2S 

3307301 S-2 

4107320 KSO4- 

4407320 LiSO4- 

4605800 MgPO4- 

4605801 MgH2PO4+ 

4605802 MgHPO4 

4703300 MnOH+ 

4703301 Mn(OH)3-1 

Table 32. Aqueous species for which MINT

identical to those in wateq4f.dat, a database 

MINTEQ ID Species 
90 One-Dimensional Transport with Equilibrium Chemistry (OTE

4704920 Mn(NO3)2 

4907320 NH4SO4- 

5001400 NaCO3- 

5007320 NaSO4- 

5401300 NiBr+ 

5401401 NiCO3 

5401402 Ni(CO3)2-2 

5401800 NiCl+ 

5401801 NiCl2 

5402700 NiF+ 

5403300 NiOH+ 

5403301 Ni(OH)2 

5403302 Ni(OH)3- 

5407320 NiSO4 

5407321 Ni(SO4)2-2 

6001300 PbBr+ 

6001301 PbBr2 

6001400 Pb(CO3)2-2 

6001401 PbCO3 

6001800 PbCl+ 

6001801 PbCl2 

6001802 PbCl3- 

6001803 PbCl4-2 

6002700 PbF+ 

6002701 PbF2 

6002702 PbF3- 

6002703 PbF4-2 

6003300 PbOH+ 

6003301 Pb(OH)2 

Table 32. Aqueous species for which MINTEQ 

identical to those in wateq4f.dat, a database dis

MINTEQ ID Species 

[
6003302 Pb(OH)3- 

6003303 Pb2OH+3 

6003304 Pb3(OH)4+2 

6003305 Pb(OH)4-2 

6004920 PbNO3+ 

6007320 PbSO4 

6007321 Pb(SO4)2-2 

7702700 SiF6-2 

8913301 U(OH)2+2 

8913302 U(OH)3+ 

8913303 U(OH)4 

8935804 UO2H2PO4)3 

9501300 ZnBr+ 

9501301 ZnBr2 

9501401 ZnCO3 

9501402 Zn(CO3)2-2 

9501800 ZnCl+ 

9501801 ZnCl2 

9501802 ZnCl3- 

9501803 ZnCl4-2 

9501804 ZnOHCl 

9502700 ZnF+ 

9503300 ZnOH+ 

9503301 Zn(OH)2 

9503302 Zn(OH)3- 

9503303 Zn(OH)4-2 

9507320 ZnSO4 

9507321 Zn(SO4)2-2 

Table 32. Aqueous species for which MINT

identical to those in wateq4f.dat, a database 

MINTEQ ID Species 
92 One-Dimensional Transport with Equilibrium Chemistry (OTE

Table 33. Aqueous species for which MINTEQ version 3 enthalpy 

distributed with PHREEQC. Enthalpy and logK values have been cha

[kcal/mole, kilocalories per mole]

MINTEQ ver. 3

(Allison and others, 1991)

& A

MINTEQ ID Species Enthalpy

[kcal/mole] LogK [

302700 AlF+2 0.000 7.010 1

302701 AlF2+ 20.000 12.750 1

302702 AlF3 2.500 17.020 2

302703 AlF4- 0.000 19.720 2

303300 AlOH+2 11.899 -4.990 11

303301 Al(OH)2+ 0.000 -10.100 0

303302 Al(OH)4- 44.060 -23.000 41

303303 Al(OH)3 0.000 -16.000 0

307320 AlSO4+ 2.150 3.020 6

307321 Al(SO4)2- 2.840 4.920 2

1507320 CaSO4 1.470 2.309 1

2311803 CuCl4-2 7.780 -4.590 7

2813301 FeOH2+ 0.000 -5.670 0

2813302 FeOH3 0.000 -13.600 24

2813303 FeOH4- 0.000 -21.600 0

2817320 FeSO4+ 3.910 3.920 5

2817321 Fe(SO4)2- 4.600 5.420 4

3300600 H2AsO3- 6.560 -9.228 6

3300601 HAsO3-2 14.199 -21.330 14

3300602 AsO3-3 20.250 -34.744 20

3300611 H2AsO4- -1.690 -2.243 -1

3300612 HAsO4-2 -0.920 -9.001 -0

3300613 AsO4-3 3.430 -20.597 3

3301400 HCO3- -3.617 10.330 -3

3301401 H2CO3 -2.247 16.681 -5

3302700 HF 3.460 3.169 3

3307620 HSeO4-1 4.200 1.906 5

3307700 H3SiO4- 8.935 -9.930 4

3307701 H2SiO4-2 29.714 -21.619 14

4603300 MgOH+ 15.935 -11.790 16

4701801 MnCl2 0.000 0.041 0

4701802 MnCl3- 0.000 -0.305 0

5002700 NaF 0.000 -0.790 2
Table 34. Mineral species for which MINTEQ 

to those in wateq4f.dat, a database distributed 

MINTEQ ID Species 

16001 GAMMA CD 

1002000 ACANTHITE 

1016000 GREENOCKITE 

1023000 CHALCOCITE 

1023001 DJURLEITE 

1023002 Anilite 

1023003 BLAUBLEI II 

1023100 BLAUBLEI I 

1023102 CHALCOPYRITE 

1028000 FES PPT 

1028001 GREIGITE 

1028002 MACKINAWITE 

1028003 PYRITE 

1047000 MNS GREEN 

1054000 MILLERITE 

1095000 ZNS (A) 

1095001 SPHALERITE 

1095002 WURTZITE 

2002000 AG2O 

2016000 CD(OH)2 (A) 

2016001 CD(OH)2 (C) 

2023000 CUPRITE 

2023100 CU(OH)2 

2023101 TENORITE 

2023102 DIOPTASE 

2028100 FERRIHYDRITE 

2028101 FE3(OH)8 
94 One-Dimensional Transport with Equilibrium Chemistry (OTE

2054000 NI(OH)2 

2054001 BUNSENITE 

2060000 MASSICOT 

2060001 LITHARGE 

2060002 PBO,.3H2O 

2060003 PLATTNERITE 

2060004 PB(OH)2 (C) 

2060005 PB2O(OH)2 

2077001 CRISTOBALITE 

2077003 SIO2(A,GL) 

2089300 UO3 (C) 

2089301 GUMMITE 

2089302 B-UO2(OH)2 

2089303 SCHOEPITE 

2095000 ZN(OH)2 (A) 

2095001 ZN(OH)2 (C) 

2095002 ZN(OH)2 (B) 

2095003 ZN(OH)2 (G) 

2095004 ZN(OH)2 (E) 

2095005 ZNO(ACTIVE) 

2095006 ZINCITE 

3006100 AS2O5 

3023000 CUPROUSFERIT 

3023100 CUPRICFERIT 

3028000 MAGNETITE 

3028100 HEMATITE 

3028101 MAGHEMITE 

3047100 BIXBYITE 

Table 34. Mineral species for which MINTEQ ve

to those in wateq4f.dat, a database distributed w

MINTEQ ID Species 
3050000 NATRON 

3060000 PB2O3 

3060001 MINIUM 

3089100 U4O9 (C) 

4002000 BROMYRITE 

4016000 CDBR2,4H2O 

4023000 CUBR 

4060000 PBBR2 

4060001 PBBRF 

4095000 ZNBR2, 2H2O 

4102000 CERARGYRITE 

4116000 CDCL2 

4116001 CDCL2,1H2O 

4116002 CDCL2,2.5H2O 

4116003 CDOHCL 

4123000 NANTOKITE 

4123100 MELANOTHALLI 

4123101 ATACAMITE 

4128100 FEOH)2.7CL.3 

4147000 MNCL2,4H2O 

4150000 HALITE 

4160000 COTUNNITE 

4160001 MATLOCKITE 

4160002 PHOSGENITE 

4160003 LAURIONITE 

4160004 PB2(OH)3CL 

4195000 ZNCL2 

4195001 ZN2(OH)3CL 

Table 34. Mineral species for which MINTEQ 

to those in wateq4f.dat, a database distributed 

MINTEQ ID Species 
96 One-Dimensional Transport with Equilibrium Chemistry (OTE

4195002 ZN5(OH)8CL2 

4202000 AGF.4H2O 

4210000 BAF2 

4216000 CDF2 

4223000 CUF 

4223100 CUF2 

4223101 CUF2,2H2O 

4260000 PBF2 

4280000 SRF2 

4289100 UF4 (C) 

4289101 UF4.2.5H2O 

4295000 ZNF2 

4302000 IODYRITE 

4306000 ASI3 

4316000 CDI2 

4323000 CUI 

4360000 PBI2 

4395000 ZNI2 

5002000 AG2CO3 

5015003 HUNTITE 

5023100 CUCO3 

5046000 ARTINITE 

5046002 MAGNESITE 

5046003 NESQUEHONITE 

5050001 THERMONATR 

5054000 NICO3 

5060000 CERRUSITE 

5060001 PB2OCO3 

Table 34. Mineral species for which MINTEQ ve

to those in wateq4f.dat, a database distributed w

MINTEQ ID Species 
5060002 PB3O2CO3 

5060003 HYDCERRUSITE 

5095000 SMITHSONITE 

5095001 ZNCO3,1H2O 

5123100 CU2(OH)3NO3 

5216000 CD(BO2)2 

5260000 PB(BO2)2 

5295000 ZN(BO2)2 

6002000 AG2SO4 

6016000 CD3(OH)4SO4 

6016001 CD3OH2(SO4)2 

6016002 CD4(OH)6SO4 

6016003 CDSO4 

6016004 CDSO4, 1H2O 

6016005 CDSO4,2.7H2O 

6023000 CU2SO4 

6023100 ANTLERITE 

6023101 BROCHANTITE 

6023102 LANGITE 

6023103 CUOCUSO4 

6023104 CUSO4 

6023105 CHALCANTHITE 

6041000 ALUM K 

6046000 EPSOMITE 

6047000 MNSO4 

6047100 MN2(SO4)3 

6050001 MIRABILITE 

6050002 THENARDITE 

Table 34. Mineral species for which MINTEQ 

to those in wateq4f.dat, a database distributed 

MINTEQ ID Species 
98 One-Dimensional Transport with Equilibrium Chemistry (OTE

6054000 NI4(OH)6SO4 

6054001 RETGERSITE 

6054002 MORENOSITE 

6060000 LARNAKITE 

6060001 PB3O2SO4 

6060002 PB4O3SO4 

6060003 ANGLESITE 

6060004 PB4(OH)6SO4 

6095000 ZN2(OH)2SO4 

6095001 ZN4(OH)6SO4 

6095002 ZN3O(SO4)2 

6095003 ZINCOSITE 

6095004 ZNSO4, 1H2O 

6095005 BIANCHITE 

6095006 GOSLARITE 

7002000 AG3PO4 

7010000 URANOCIRCITE 

7015000 NINGYOITE 

7015001 AUTUNITE 

7015002 FCO3APATITE 

7016000 CD3(PO4)2 

7023100 CU3(PO4)2 

7023101 CU3(PO4)2,3W 

7023102 TORBERNITE 

7028000 BASSETITE 

7028001 VIVIANITE 

7028100 STRENGITE 

7041000 K-AUTUNITE 

Table 34. Mineral species for which MINTEQ ve

to those in wateq4f.dat, a database distributed w

MINTEQ ID Species 
7046000 SALEEITE 

7047000 MN3(PO4)2 

7049000 URAMPHITE 

7050000 NA-AUTUNITE 

7054000 NI3(PO4)2 

7060000 PRZHEVALSKIT 

7060001 CLPYROMORPH 

7060002 HXYPYROMORPH 

7060003 PLUMBGUMMITE 

7060004 HINSDALITE 

7060005 TSUMEBITE 

7080000 SR-AUTUNITE 

7089301 H-AUTUNITE 

7095000 ZN3(PO4),4W 

8015000 URANOPHANE 

8054000 NI2SIO4 

8060000 PB2SIO4 

8095000 WILLEMITE 

8216000 CDSIO3 

8260000 PBSIO3 

8295000 ZNSIO3 

8450000 MAGADIITE 

8628000 GREENALITE 

1

Actual values in wateq4f.dat are equal to -1 times t

mineral reactions are written within PHREEQC.

Table 34. Mineral species for which MINTEQ 

to those in wateq4f.dat, a database distributed 

MINTEQ ID Species 
100 One-Dimensional Transport with Equilibrium Chemistry (OT

 

Table 35. Mineral species for which MINTEQ version 

wateq4f.dat, a database distributed with PHREEQC. Def

file are equal to the MINTEQ version 3 values.

[kcal/mole, kilocalories per mole]

PHR

(wateq4f.dat, Par

19

MINTEQ ID Species Enthalpy

[kcal/mole] 

2003000 ALOH3(A) 26.500 

6041001 ALUNITE 50.250 

6015000 ANHYDRITE 1.710 

5015000 ARAGONITE 2.589 

6010000 BARITE -6.350 

2003001 BOEHMITE 28.181 

2046000 BRUCITE 27.100 

5015001 CALCITE 2.297 

6080000 CELESTITE 1.037 

2077000 CHALCEDONY -4.720 

8646000 CHRYSOTILE 46.800 

8246000 CLINOENSTITE 20.049 

2003002 DIASPORE 24.681 

8215000 DIOPSIDE 32.348 

5015002 DOLOMITE 9.436 

4215000 FLUORITE -4.690 

8046000 FORSTERITE 48.578 

2003003 GIBBSITE(C) 22.800 

2028102 GOETHITE 14.480 

6015001 GYPSUM 0.109 

3047000 HAUSMANNITE 100.640 

5046001 HYDRMAGNESIT 52.244 

8603001 KAOLINITE 35.300 

6028000 MELANTERITE -4.910 
2015001 PORTLANDITE 31.000 

2047003 PYROCROITE 

2077002 QUARTZ -5.990 

1006001 REALGAR -30.894 

5089300 RUTHERFORDIN 1.440 

8646003 SEPIOLITE(C 10.700 

2077004 SIO2(A,PT) -3.340 

5080000 STRONTIANIT 0.400 

8215001 TREMOLITE 96.853 

3089101 U3O8(C) 116.000 

2089100 URANINITE 18.610 

5010000 WITHERITE -0.703 

1

Actual values in wateq4f.dat are equal to −1 times the valu

tions are written within PHREEQC.

Table 35. Mineral species for which MINTEQ vers

wateq4f.dat, a database distributed with PHREEQC. 

file are equal to the MINTEQ version 3 values.—Con

[kcal/mole, kilocalories per mole]

P

(wateq4f.dat, 

MINTEQ ID Species Enthalpy

[kcal/mole

Publishing support provided by:

Denver Publishing Service Center

For more information concerning this publication, contact:

Director, USGS Colorado Water Science Center

Box 25046, Mail Stop 415

Denver, CO 80225

(303)236-4882

Or visit the Colorado Water Science Center Web site at:

http://co.water.usgs.gov/
Runkel—One-Dimensional Transport with Equilibrium Chemistry (OTEQ) — A Reactive Transport Model for Streams and Rivers—Techniques and Methods 6–B6

Printed on recycled paper