NAME
phreeqc - A program for speciation, batch-reaction, one-
dimensional transport, and inverse geochemical calculations
SYNOPSIS
phreeqc [infile [outfile [database [screen_output]]]
OPTIONS
infile the name of the PHREEQC input file
outfile the name of the file to which PHREEQC output
will be written
database the name of the PHREEQC database
screen_output the name of the file to which screen output
will be directed
If no arguments are specified, the program prompts for the
input, output, and database file names.
If only infile is specified, then outfile defaults to
infile.out. If no database file is specified, the
phreeqc.dat database distributed with PHREEQC will be used.
ABSTRACT
PHREEQC is a computer program written in the C programming
language that is designed to perform a wide variety of low-
temperature aqueous geochemical calculations. PHREEQC is
based on an ion-association aqueous model and has
capabilities for (1) speciation and saturation-index
calculations; (2) batch-reaction and one-dimensional (1D)
transport calculations involving reversible reactions, which
include aqueous, mineral, gas, solid-solution, surface-
complexation, and ion-exchange equilibria, and irreversible
reactions, which include specified mole transfers of
reactants, kinetically controlled reactions, mixing of
solutions, and temperature changes; and (3) inverse
modeling, which finds sets of mineral and gas mole transfers
that account for differences in composition between waters,
within specified compositional uncertainty limits.
METHOD
For speciation and batch-reaction calculations, PHREEQC
solves sets of nonlinear mole-balance and mass-action
equations that define an ion-association model. A Newton-
Raphson formulation is used to iteratively arrive at a
solution to the equations. A robust numerical method is
provided by using an optimizing solver that allows both
equality and inequality equations. The solver is used to
obtain the intermediate estimates of changes in the unknowns
at each iteration.
For inverse modeling, a set of linear mole-balance equations
are solved. The equations contain additional unknowns that
account for uncertainty in the analytical data. The
optimizing solver is used to solve the linear equations
while maintaining the uncertainty terms within specified
limits.
For transport modeling, the partial differential equations
of transport are solved by an operator splitting scheme that
sequentially solves for advective and dispersive transport,
followed by chemical equilibration that is equivalent to
batch-reaction calculations for each cell. Time steps are
selected to maintain numerical accuracy. If kinetic
reactions are modeled, yet another splitting of operators is
implemented and a 5th order Runge-Kutta method is used to
integrate the ordinary differential equations of the kinetic
reactions.
HISTORY
Version 2.6 2002/04/22 - New batch release to correspond with
initial release of PhreeqcI, an graphical user interface
to PHREEQC that runs under Windows operating systems.
Changes to PHREEQC include: (1) All selected_output is
routed through a single routine. (2) Underscores "_" are
allowed inside square brackets, [A_bcd] for element
definitions. (3) Fixed bug match_elts_in_species, check
for "e-" was wrong. (4) Modified minteq.dat to put
CuS4S5-3, Cu(S4)2-3 in Cu(1) mole balance equations
instead of Cu(2). Before the change, the program would
not converge if Cu(2) were defined in an initial solution.
(5) Revisions to improve SOLID_SOLUTIONS convergence with
small numbers of moles of solids. (6) Made changes related
to dump file and PhreeqcI. (7) Printed number of iterations
now sums iterations in all kinetics calculations. (8) Fixed
bug with LA("H2O"), which was returning natural log of
activity of water.
Version 2.5 2001/10/01 - (1) In llnl.dat, fixed sign errors
in REE (rare earth elements) for some redox reactions and
removed some redundant species, generally ReeO2- was
retained and Ree(OH)4- was removed. (2) Added the
capability to use square brackets to define an "element"
name. The brackets act like quotation marks in that any
character string can be used within the brackets as an
element name. This was introduced to simplify expansion
of the model to isotopic species. [13C], [14C], and
[18O] are legal element names. (3) Added identifier
-activity_water for a species in SOLUTION_SPECIES data
block. This identifier has been added for future updates
that will allow isotopic calculations. It is intended to
be used only for isotopic variations of H2O, like D2O or
H2[O18]. It forces the activity coefficient for the
species to be activity(water)/55.5. This effectively
sets the activity of the species to the mole fraction in
solution. (4) Fixed bug in checking solid solutions for
presence or absence of elements in the system.
Programming error caused segmentation fault if an error
was detected under certain conditions. (5) Changed
return value of MOL to be molality of water if argument
is "H2O". Also changed return value of LA to be activity
of water if argument is "H2O". (6) Diffuse layer
calculation was incorrect if aqueous phase did not have 1
kilogram of water. Eq. 74 of manual has molality, but
code used moles. The code was corrected by adding the
mass of water to the formulation. (7) Stagnant zones
with first-order exchange approximation (1 stagnant cell,
exchange factor, and porosities defined) did not work
correctly if mobile and immobile cells did not have equal
volumes of water. The mixing factors were revised to
account for the masses of water in the stagnant and
mobile zones. (8) A fatal error was erroneously detected
if the database file had a DATABASE data block. DATABASE
data block is now ignored while reading the database
file. (9) Added identifier -bad_step_max to KINETICS
data block. An integer following -bad_step_max gives the
maximum number of times a rate integration may fail
before execution of the program is terminated. Default
is 500.
Version 2.4.2 2001/06/15 - Fixed spreadsheet bug. Program
was not ignoring columns that could not be identified as
either element names or allowed data (ph, pe, number,
description, etc.). Also, the program failed if a
spreadsheet solution number was negative.
Version 2.4.1 2001/06/04 - Fixed spreadsheet bugs with
isotopes.
Version 2.4 2001/06/01 - (1) Added structure for spreadsheet
for use by PhreeqcI. (2) Isotope value initialized
incorrectly if only an -uncertainty was defined in
SOLUTION_SPREAD. (3) Fixed segmentation violation when
primary and secondary master species were defined
improperly. (4) Corrected enthalpies of reaction in
llnl.dat. Previous release had erroneously had
enthalpies of formation in -delta_H parameter; the values
should be enthalpies of reaction. Enthalpies of reaction
were calculated from the enthalpies of formation and
these values are now included in the -delta_H parameter.
This change will have very little impact on calculations
because the analytical expression has precedence over
-delta_H in calculating temperature dependence of log K,
and nearly all species and minerals have an analytical
expression or lack both an analytical expression and an
enthalpy of reaction. (5) Corrected bugs in punch of
solid solution components that caused both selected
output and output file errors: moles were incorrect in
selected output, and total moles and mole fraction were
incorrect in output file. (6) Added surface complexation
constants for Fe+2; two complexes for weak sites and one
complex for strong sites. phreeqc.dat and wateq4f.dat
modified. (7) Comment for units of parameters for
calcite rate equation was wrong. Rate equation now uses
cm^2/L for area parameter. Previously the correct units
would have been 1/decimeter. phreeqc.dat and wateq4f.dat
modified. (8) Fixed a bug when rates were equal within
tolerance, but negative concentrations occurred because
of small initial concentrations. (9) Added -warnings to
PRINT keyword for specification of maximum number of
warnings to print. Negative number allows all warnings
to be printed. (10) Function CELL_NO in Basic now prints
a number equivalent to -solution in SELECTED_OUTPUT data
block. This does not change printing for ADVECTION or
TRANSPORT calculations. (11) Kinetics time is halved for
advective part of reaction in transport; time incorrectly
accounted for before. (12) -punch_ identifiers printed
-1 instead of the correct solution number for batch-
reaction calculations. (13) -high_precision is no longer
reset to false with every SELECTED_OUTPUT data block.
(14) SELECTED_OUTPUT file name stored for use by
PhreeqcI. (15) Alkalinity for NH3 corrected to 1.0 in
llnl.dat. (16) Fixed bug with USER_PRINT of kinetics.
Did not find correct kinetics information in some cases.
(17) Fixed bug in default values for SOLUTION_SPREAD.
Cannot use phase name and SI for pH or pe, and bug did
not allow PHREEQC to run. Now PHREEQC runs, but warns
that this is not allowed.
Version 2.3 2001/01/02 - (1) Added new keyword DATABASE. It
must be the first keyword in the input file. The
character string following the keyword is the pathname
for the database file to be used in the calculation. The
file that is specified takes precedence over any default
database name, including environmental variable
PHREEQC_DATABASE and command line arguments. (2) Fixed
bug in SOLUTION_SPREAD. If first heading in the
spreadsheet input was an identifier--pH, pe, units,
etc.--then the headings were interpreted as an identifier
and bad things happened. (3) Added new keyword to make
aqueous model similar to LLNL and Geochemists Workbench
when using llnl.dat as the database file. Values of
Debye-Huckel a and b and bdot (ionic strength
coefficient) are read at fixed temperatures. Linear
interpolation occurs between temperatures. (4) New
options for SOLUTION_SPECIES are -llnl_gamma a, where a
is the ion-size parameter. -co2_llnl_gamma, indicates
the temperature dependent function for the bdot term
given in -co2_coefs of LLNL_AQUEOUS_MODEL_PARAMETERS will
be used. Applies to uncharged species only. (5) Fixed
bug in basic interpreter. A number like "..524" would
cause an infinite loop. (6) Added function SURF to
Basic. SURF("element", "surface") gives the amount of
element in the diffuse layer for "surface". "surface"
should be the surface name, not the surface-site name
(that is, no underscore). (7) Fixed option to
"runge_kutta" from "runge-kutta" to match documentation
for KINETICS. (8) Fixed UO2+2 and Mn+2 reaction
stoichiometry for Hfo surface complexation in
wateq4f.dat. (9) Added option for an equilibrium-phase
to dissolve only. "dis" is added at the end of a line
defining an equilibrium-phase. No data fields may be
omitted. Should not be used when adding an alternative
reaction. (10) R-K integration failed when only the
final rate generated negative concentrations. (11) Allow
decimals in definition of secondary master species, for
example S(0.3). (12) Fixed bug if description was more
than about 85 characters; now allows about 400
characters. (13) Fixed bug for surface/exchange sites
related to phases. Was checking internal copies of
surfaces/exchange with negative numbers. (14) Fixed bug
in quick prep that did not set the correct pointer for
gas phases. (15) Fixed segmentation fault that occurred
if all elements for phase-boundary mineral were not in
the solution. Only applied to a phase used to define
concentration in an initial solution calculation. (16)
Added option to eliminate echo of input file in PRINT
data block. -echo_input T/F turns echoing on and off.
Default is on.
Version 2.2 2000/03/01 - (1) Fixed bug in MIX if no
solutions are defined. (2) Changed printout for surface;
only gives net surface charge when -diffuse_layer is
defined; prints correct value for the surface charge and
surface charge density for diffuse-layer calculation.
(3) Added new function EDL to Basic, which gives
information about surfaces. EDL("element", "surface")
gives the amount of element in the diffuse layer for
"surface" if -diffuse_layer is used, zero otherwise. The
"surface" should be the surface name, not the surface-
site name (that is, no underscore). Special values for
"element" include: "charge"--gives surface charge,
equivalents (all options), "sigma"--surface charge
density, C/m**2 (default or -diffuse_layer options),
"psi"--potential at the surface, Volts (default or
-diffuse_layer options), "water"--mass of water in the
diffuse layer, kg (-diffuse_layer option). Value of
function is zero for any undefined quantities. (4)
Changed distribution to be more consistent with other
USGS water-resources applications.
Version 2.1 2000/01/19 - Fixed problem in DOS version
preventing last digit of 3-digit exponent from printing
for PUNCH and PRINT commands included in USER_PUNCH and
USER_PRINT input data blocks, respectively.
Version 2.0 1999/12/15 - Enhancements added to extend
capabilities to simulate dispersion (or diffusion) and
stagnant zones in 1D-transport calculations, to model
kinetic reactions with user-defined rate expressions, to
model the formation or dissolution of ideal,
multicomponent or nonideal, binary solid solutions, to
model fixed-volume gas phases in addition to fixed-
pressure gas phases, to allow the number of surface or
exchange sites to vary with the dissolution or
precipitation of minerals or kinetic reactants, to
include isotope mole balances in inverse modeling
calculations, to automatically use multiple sets of
convergence parameters, to print user-defined quantities
to the primary output file and (or) to a file suitable
for importation into a spreadsheet, and to define
solution compositions in a format more compatible with
spreadsheet programs.
Version 1.6 1998/01/16 - Code fixes: (A) Scaled 1/u factors
by 1e-3 in inverse modeling, may help with precision
problems in solver. (B) Fixed bug in inverse modeling.
If carbon alkalinity was small relative to total
alkalinity, speciation in carbon-derivative calculation
failed. (C) Changed parameterization of Debye-Huckel A
and B to be consistent with Wateq4f. (D) Fixed pointer
fault if specified solution does not exist for initial
surface or initial exchange calculation. (E) Changed
databases: (1) Units for Delta_H for H2 and Goethite
incorrectly defaulted to kJ/mol. Now explicitly
kcal/mol. (2) Added analytical expression for H2S, NH3,
KSO4. (3) Added species CaHSO4+. (4) Added delta H for
Goethite. (5) changed logK for jarosite(ss). (F) Check
for misuse of "as CaCO3" for alkalinity. (G) Revise
guesses does not fail after 30 iterations. Now goes
through another 30 iterations to make calculated moles
less than total, but not necessarily within 5 orders of
magnitude of the total. (H) Added logic to use " " file
names. (I) Added logic to allow "Alkalinity" to be
printed to the selected output file through -totals
option of keyword SELECTED_OUTPUT. (J) Fixed bug in
SELECTED_OUTPUT where log activities were not printed
correctly for non-master species. (K) Fixed bug where
space for exchange and surface was always freed, but not
always allocated. (L) Fixed bug in selected_output file
where solution number was not printed correctly for
reaction and transport calculations. (M) Now stops
execution if error occurs in get_elements_in_species.
(N) In EQUILIBRIUM_PHASES, only accepted alternate-
formula phase name if it began with a capital letter.
Version 1.5 1997/01/31 - Code fixes: (A) after a reaction
calculation, saving a range of exchangers, surfaces, gas
phases, or equilibrium_phases now functions correctly;
(B) in reaction calculations, diffuse_layer calculations
is now reset properly as the value for surface, if
surface is used, or to FALSE if not used; (C) in
transport simulations, for cell n, mixing was only
correct using cells j < n. Cells j > n, had already been
updated with mixing and reactions. Now cell n+1, remains
as it was before mixing and reactions until after mixing
and reactions are calculated for cell n; (D) now
correctly handles situations where both mineral amounts
and aqueous concentrations for an element were very
small, in the range of 1e-14. Some parameters were
adjusted to try to improve the selection of minerals to
include in the Jacobian. As has been the case, but maybe
not stated, it is possible to lose mass at these
concentrations, because amounts of minerals less than
1e-14 are set to 0.0; (E) when an element was missing
from the system (not in solution, but mineral with zero
mass was specified), a warning message that used a NULL
pointer has been fixed.
Version 1.4 1996/10/18 - Fixed an error that occurred during
transport or reaction calculations. If a concentration
of an element was less than 1e-14 (Mn in equilibrium with
pyrolusite in the specific example), then under some
conditions a warning message would be printed ("WARNING:
Element in phase, Pyrolusite, is not in model.") and
erroneous results could be calculated. Error message was
changed to include the offending element. Fixed an error
in inverse modeling that occurred when dissolved oxygen
or O2(g) were part of an inverse model. The
stoichiometry of electron and alkalinity change per mole
of oxygen were wrong. This could have caused erroneous
mole-balance models when using versions 1.3 and lower.
Similar errors would occur with any element or valence
state that has a master species with stoichiometry 2 or
greater for the element or valence state (N2 for
example). Added analytical expressions to phreeqc.dat
and wateq4f.dat aqueous species NH3, KSO4-, HS-. Added
delta_H for goethite. Added species CaHSO4+ to
phreeqc.dat. All data from WATEQ4F manual.
Version 1.3 1996/07/30 - Included minteq.dat, a translated
MINTEQA2 database for all but DOM (dissolved organic
matter) species and some species that could not be
deciphered. Added check to make sure minimum number of
species are defined, e-, H+, H2O, H2, O2, and minimum
number of master species are defined, H(0), H(1), O(-2),
O(0), E. Removed limitation on number of strings. Also
upped initial dimension from 1500 strings to 3000. Fixed
various errors in the program including the following.
Fixed logic resulting in segmentation fault when parsing
equations. Fixed incomplete error message for invalid
output file name. Fixed divide by zero problem when
activity of a master species in an intermediate iteration
was very small. In inverse modeling, fixed bug with
-balances that did not interpret valence states with a
"+" sign correctly. N(+5) caused a segmentation fault,
although N(5) executed correctly.
Version 1.2 1996/03/22 - Fixed bug if number of aqueous
species for a solution > 250 added code to dynamically
allocate sufficient space. Fixed most warnings from
Borland C++ compiler. Fixed bug with opening files, if
nonexistent file is given. Closes all files before
exiting. Acrobat version of documentation included in
distribution.
Version 1.1 1995/11/22 - Fixed bug with "use" data.
Selection process for reactions was not working
correctly. Wrong solution could be used in reaction
calculation. Used purify to remove memory errors.
Version 1.0 1995/10/12 - First release.
DATA REQUIREMENTS
Proper use of the program requires adequate knowledge of
geochemistry and a proper formulation of the problem. Input
is arranged in keyword data blocks, which can appear in any
order. Data fields for a keyword are read in a free format,
thus they are not column dependent.
For speciation modeling, analytical data for a solution
composition (SOLUTION keyword) are needed.
For batch-reaction modeling, the initial solution
composition is required (SOLUTION or MIX data block). Other
equilibrium reactants may be defined with
EQUILIBRIUM_PHASES, EXCHANGE, SURFACE, GAS_PHASE, and
SOLID_SOLUTION data blocks. Nonequilibrium reactions may be
defined with KINETICS and RATES, REACTION, and
REACTION_TEMPERATURE data blocks.
For 1D transport modeling, the data for batch-reaction
modeling are needed for each cell in the modeled system. In
addition, physical information is needed about column
dimensions, time steps, boundary conditions, and
dispersivity.
For inverse modeling, the solution composition of the final
solution and one or more initial solutions are needed
(SOLUTION data block). Uncertainty limits must be defined
explicitly or by default for each element and element redox
state in the solutions. In addition, the identity and
composition of a set of plausible reactants and products are
needed.
Three default databases are included that contain the
definition of aqueous species, exchange species, surface
species, and mineral and phases for a set of elements. The
database phreeqc.dat contains information for Al, B, Ba, Br,
C, Ca, Cd, Cl, Cu, F, Fe, H, K, Li, Mg, Mn, N, Na, O, P, S,
Si, Sr, Zn. The database wateq4f.dat contains the
additional constituents Ag, As, Cs, Fulvate, Humate, I, Ni,
Rb, Se, and U. The database minteq.dat is derived from the
thermodynamic data of the program MINTEQA2. If additional
elements, species, or phases are needed, then chemical
reactions, log K, and data for the temperature dependence of
log K are needed for each additional species and phase.
SYSTEM REQUIREMENTS
PHREEQC is written in ANSI C. Generally, the program is
easily installed on most computer systems. The code has
been used on UNIX-based computers and on IBM-compatible
computers with processors running at 100 megahertz or
faster.
DOCUMENTATION
Parkhurst, D.L., and Appelo, C.A.J., 1999, User's guide to
PHREEQC (Version 2)--a computer program for speciation,
batch-reaction, one-dimensional transport, and inverse
geochemical calculations: U.S. Geological Survey Water-
Resources Investigations Report 99-4259, 312 p.
RELATED DOCUMENTATION
Charlton, S.R., Macklin, C.L. and Parkhurst, D.L., 1997,
PHREEQCI--a graphical user interface for the geochemical
computer program PHREEQC: U.S. Geological Survey Water-
Resources Investigations Report 97-4222, 9 p.
Charlton, S.R., and Parkhurst, D.L., 2002, PhreeqcI--A graphical user
interface to the geochemical model PHREEQC: U.S. Geological Survey
Fact Sheet FS-031-02, 2 p.
Parkhurst, D.L., Thorstenson, D.C., and Plummer, L.N., 1980,
PHREEQE--a computer program for geochemical calculations:
U.S. Geological Survey Water-Resources Investigations
Report 80-96, 195 p. (Revised and reprinted, 1990.)
Plummer, L.N., Parkhurst, D.L., Fleming, G.W., and Dunkle,
S.A., 1988, A computer program incorporating Pitzer's
equations for calculation of geochemical reactions in
brines: U.S. Geological Survey Water-Resources
Investigations Report 88-4153, 310 p.
Plummer, L.N., Prestemon, E.C., and Parkhurst, D.L., 1991,
An interactive code (NETPATH) for modeling NET
geochemical reactions along a flow PATH: U.S.
Geological Survey Water-Resources Investigations Report
91-4078, 227 p.
Plummer, L.N., Prestemon, E.C., and Parkhurst, D.L., 1994,
An interactive code (NETPATH) for modeling NET
geochemical reactions along a flow PATH--version 2.0:
U.S. Geological Survey Water-Resources Investigations
Report 94-4169, 130 p.
REFERENCES
Appelo, C.A.J., and Postma, D., 1993, Geochemistry,
groundwater and pollution: Rotterdam, Netherlands, and
Brookfield, Vermont, A.A. Balkema.
Appelo, C.A.J., and Willemsen, A., 1987, Geochemical
calculations and observations on salt water intrusions.
I: A combined geochemical/mixing cell model: Journal of
Hydrology, v. 94, p. 313-330.
Parkhurst, D.L., and Plummer, L.N., 1993, Geochemical
models, in Alley, W.M., ed., Regional ground-water
quality: New York, Van Nostrand Reinhold, chap. 9, p.
199-225.
Plummer, L.N., 1984, Geochemical modeling: A comparison of
forward and inverse methods, in Hitchon, B., and Wallick,
E.I., eds., Proceedings First Canadian/American
Conference on Hydrogeology--Practical Applications of
Ground Water Geochemistry, Banff, Alberta, Canada:
Worthington, Ohio, National Water Well Association, p.
149-177.
TRAINING
PHREEQC is taught as part of the courses Geochemistry for
Ground-Water Systems (GW3021TC) at the USGS National
Training Center.
CONTACTS
Operation:
U.S. Geological Survey
David Parkhurst
Denver Federal Center, MS 413
Lakewood, CO 80225
dlpark@usgs.gov
Distribution:
U.S. Geological Survey
Hydrologic Analysis Software Support Program
437 National Center
Reston, VA 20192
h2osoft@usgs.gov
Official versions of U.S. Geological Survey water-resources
analysis software are available for electronic retrieval via
the World Wide Web (WWW) at:
http://water.usgs.gov/software/
and via anonymous File Transfer Protocol (FTP) from:
water.usgs.gov (path: /pub/software).
The WWW page and anonymous FTP directory from which the
PHREEQC software can be retrieved are, respectively:
http://water.usgs.gov/software/phreeqc.html
--and--
/pub/software/geochemical/phreeqc
See
http://water.usgs.gov/software/ordering_documentation.html
for information on ordering printed copies of USGS
publications.
SEE ALSO
netpath(1) - Interactive program for calculating NET
geochemical reactions and radiocarbon dating
along a flow PATH
phreeqci(1) - Graphical user interface for PHREEQC
phrqpitz(1) - A program for geochemical calculations in
brines
wateq4f(1) - A program for calculating speciation of major,
trace, and redox elements in natural waters