max/iphreeqc
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/iphreeqc.git synced 2025-12-15 16:18:22 +01:00
Code Issues Packages Projects Releases Wiki Activity
iphreeqc/revisions
David L Parkhurst 60a1544019 Copying new classes (cxx) and cpp files to src 
Will remove cpp and header files and make phreeqc an external directory.




git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/trunk@785 1feff8c3-07ed-0310-ac33-dd36852eb9cd
2006-02-16 00:21:39 +00:00

1402 lines
61 KiB
Plaintext
Raw Blame History

------------------------------------------------------------
Features not documented in WRIR 99-4259.
------------------------------------------------------------
svn 675:
Added PRINT option to print the species that contribute
to alkalinity. Alkalinity distribution is printed in
the output file following the distribution of species.
Default at program startup is false.
PRINT
-alkalinity true
------------------------------------------------------------
Version 2.12:
------------------------------------------------------------
* Made aqueous activity coefficients the default activity
coefficients for exchange species when using the
Pitzer formulation. New option in EXCHANGE is
-pitzer_exchange_gammas T/F, default is true;
defining "false" sets exchange activity coefficients
to 1.0. Option has no effect for ion-association
model (non-Pitzer).
* Added multiplier format to REACTION steps and KINETICS steps,
which simplifies definition of multiple equal reaction increments.
This definition:
INCREMENTAL_REACTIONS true
REACTION
H2O 1
-36 3*-4 2*-.25 -0.19 4*-0.1 3*-0.05 moles
is equivalent to this definition:
INCREMENTAL_REACTIONS true
REACTION
H2O 1
-36 -4 -4 -4 -.25 -.25 -0.19 -0.1 -0.1 -0.1
-0.1 -0.05 -0.05 -0.05 moles
* Added Pitzer activity formulation. Use pitzer.dat database
to invoke the Pitzer model. Should have same capabilities
as ion-association model except explicit diffuse layer
calculation is not implemented with the Pitzer model.
New keyword is PITZER with following options:
PITZER
-MacInnes T/F # uses MacInnes assumption or unscaled for
# individual activities and activity coefficients
-B0
Na+ Cl- 0.0765 -777.03 -4.4706 0.008946 -3.3158E-6
-B1
Na+ Cl- 0.2664 0 0 6.1608E-5 1.0715E-6
-B2
Mg+2 SO4-2 -37.23 0 0 -0.253
-C0
Na+ Cl- 0.00127 33.317 0.09421 -4.655E-5
-THETA
K+ Na+ -0.012
-LAMDA
Na+ CO2 0.1
-ZETA
H+ Cl- B(OH)3 -0.0102
-PSI
Na+ K+ Cl- -0.0018
A five-term expression for temperature dependence is available
for all Pitzer parameter values:
P = A0 + A1*(1/TK - 1/TR) + A2log(TK/TR) + A3/(TK-TR) +
A4*(TK*TK - TR*TR),
where TK is Kelvin and TR is 298.15.
* Cl1mp is a new multiple precision version of routine cl1,
a simplex-based optimization routine. Cl1mp was develeped
by using the Gnu Multiple Precision package (gmp).
Calculations are carried out to about 30 significant
digits. Cl1mp may help in some situations where roundoff
errors are a problem, but it is still possible that roundoff
errors will cause cl1mp to fail to find a solution to an
optimization problem. The mp version has the following
options in INVERSE_MODELING:
-multiple_precision T/F--causes the mp version
to be used in inverse modeling calculations.
-mp_tolerance 1e-12--tolerance for mp version of
cl1. As in cl1, numbers less than the
tolerance are considered to be zero.
1e-12 is the default.
-censor_mp 1e-20--as calculations occur in the
linear equation array, elements less
than this value are set to zero. Default
is 1e-20. A value of 0.0 causes no
censoring to occur.
------------------------------------------------------------
Version 2.11:
------------------------------------------------------------
* A new database, minteq.v4.dat, has been translated from
version 4.02 of MINTEQA2 and is included in all
distributions. The database minteq.dat from earlier
version of MINTEQA2 has been slightly revised and is
also included.
------------------------------------------------------------
Version 2.10:
------------------------------------------------------------
No new features.
------------------------------------------------------------
Version 2.9:
------------------------------------------------------------
* Added new keyword COPY that allows a data entity
to be copied from one index to a new index
or to a range of indicies. Format is
COPY keyword index index_start[-index_end]
where keyword may be SOLUTION
EQUILIBRIUM_PHASES
EXCHANGE
GAS_PHASE
KINETICS
MIX
REACTION
REACTION_TEMPERATURE
SOLID_SOLUTION
SURFACE
* Added new Basic functions
b$ = PAD(a$, 20) pads a$ to a total of 20 characters
with spaces and stores result in b$. PAD returns
a copy of a$ if a$ is more than 20 characters.
i = INSTR(a$, b$) sets i to the character position of
string b$ in a$, 0 in not found.
b$ = LTRIM(a$) trims white space from beginning of
string a$ and stores result in b$.
b$ = RTRIM(a$) trims white space from end of string
a$ and stores result in b$.
b$ = TRIM(a$) trims white space from beginning and
end of string a$ and stores result in b$.
* Added new Basic function SYS that calculates the to
total amount of an element in all phases (solution,
equilibrium_phases, surfaces, exchangers, solid solutions,
and gas phase). KINETIC reactions are not included.
The function has two forms: (1) one element name as an
argument (variable names are user specified)
10 t = SYS("As")
the function will return the total arsenic in the system.
(2) 5 argumens
10 t = SYS("As", count_species, names$, types$, moles)
will return the total arsenic in the system to tot; count_species--
the number of species that contain arsenic, including
solution, equilibrium_phases, surfaces, exchangers, solid solutions,
and gas phase species; names$--a character array that has the
name of each species; type$--a character array that specifies the
type of phase for the species, aq, equi, surf, ex, s_s, gas, diff.
Diff refers to the amount of the element in the diffuse layer of
a surface when the explicit diffuse layer calculation is used;
moles--an array containing the number of moles of the element in
the species. The sum of moles(i) is equal to tot.
SYS has several special arguments for the form
SYS("arg", count, names$, types$, values)
arg is one of the options listed below.
count is a single numeric value and is the number of elements
in the following arrays.
name$ is an array of string values.
type$ is an array of string values.
values is an array of numeric values.
Values of arg:
elt_name returns total number of moles of element in system.
count is the number of species for the element in
the system, including aqueous, exchange, surface,
equilibrium_phase, solid solution component, and
gas phase "species".
Arrays are filled for each "species"; values are moles.
"elements" returns total number of moles of dissolved elements other
than H and O.
count is number of elements, valence states,
exchangers, and surfaces.
Arrays are filled for each element and valence state,
type is "dis"; exchanger, type is "ex",
and surface, type is "surf". Values are moles.
"phases" returns maximum saturation index of all phases.
count is number of phases in system.
Arrays are filled for each phase; values are
saturation indexes.
"aq" returns sum of moles of all aqueous species.
count is number of aqueous species in system.
Arrays are filled with each aqueous species;
values are moles.
"ex" returns sum of moles of all exchange species.
count is number of exchange species in system.
Arrays are filled with each exchange species;
values are moles.
"surf" returns sum of moles of all surface species.
count is number of surface species in system.
Arrays are filled with each surface species;
values are moles.
"s_s" returns sum of moles of all solid solution components.
count is number of solid solution components in system.
Arrays are filled with each solid solution component;
values are moles.
"gas" returns sum of moles of all gas components.
count is number of gas components in system.
Arrays are filled with each gas component;
values are moles.
* Added new Basic function, DESCRIPTION, that has the value
defined for the description field of the SOLUTION keyword line.
* Added alternative ordinary differential equation solver
called CVODE, a set of C routines from the Lawrence
Livermore National Labs. CVODE is part of the SUNDIALS
package. CVODE is used in place of the Runge Kutta method
when "-cvode true" is used within a KINETICS data block.
KINETICS
-cvode true
------------------------------------------------------------
Version 2.8:
------------------------------------------------------------
No new features.
------------------------------------------------------------
Version 2.7:
------------------------------------------------------------
Changed format of selected output file:
Removed quotations surrounding strings in headings.
Removed quotations surrounding strings in state variable.
All fields are 12 or 20 places depending on
-high_precision.
Headings are not truncated even if longer than
specified precision.
For isotopes, missing value is -9999.9
Selected output is updated each simulation.
If a species or phase is defined
subsequent to the simulation where SELECTED_OUTPUT
was defined it will appear in the selected output
file in the simulation in which it is defined and
in subsequent simulations.
Added strings for each file, which can be extracted from the
executable file with the "ident" command.
Fixed null pointer for isotope_ratios if Basic routine
was undefined.
Fixed problem in C++ if structure name is same as member name.
logk member of logk structure was renamed to log_k.
Added identifier -add_constant to PHASES, EXCHANGE_SPECIES,
SOLUTION_SPECIES, and SURFACE_SPECIES.
-add_constant -0.301
log K is augmented by the specified constant.
Theory and implementation of isotopes in PHREEQC is documented in:
Thorstenson, D.C., and Parkhurst, D.L., 2002, Calculation of
individual isotope equilibrium constants for implementation in
geochemical models: U.S. Geological Survey Water-Resources
Investigations Report 02-4172, 129 p.
Added KEYWORDS:
ISOTOPES
Element
-isotope isotope_name units standard_ratio
-total_is_major T/F (OPTION IS DISABLED!!)
CALCULATE_VALUES
Name
-start
Basic statements, must have SAVE
-end
ISOTOPE_RATIOS (for printing)
Name=Calculate_values_name Isotope_name
ISOTOPE_ALPHAS (for printing)
Name=Calculate_values_name Named_logk=named_expression_name
Basic functions:
calc_value("calc_value_name") evaluates a definition of CALCULATE_VALUES
lk_named("name") log10(K) of definition in NAMED_EXPRESSIONS
lk_phase("name") log10(K) of definition in PHASES
lk_species("name") log10(K) of definition in (SOLUTION, EXCHANGE, SURFACE)_SPECIES
sum_gas("template","element") Sum of element in gases with specified template
template="{C,[13C],[14C]}{O,[18O]}2" includes all CO2 gases
sum_species("template","element") Sum of element in aqueous, exchange, and surface species with
specified template
sum_s_s("s_s_name","element") Sum of element in a specified solid solution
PRINT keyword:
-initial_isotopes T/F
-isotope_ratios T/F
-isotope_alphas T/F
-censor_species 1e-8 # omit species from Distribution of Species if less than
# relative minimum of an element or element redox state
# total concentration
SELECTED_OUTPUT keyword:
-calculate_values name1 name2 ...
-isotopes minor_isotope1 minor_isotope2 ....
Added functions LK_SPECIES, LK_NAMED, LK_PHASE for Basic
interpreter. LK_SPECIES("CaHCO3+") returns the
log k for the association reaction for the ion pair
CaHCO3+ at the current temperature. The log K is
for the reaction as defined in the database or
input file. Similarly,
LK_NAMED("Log_alpha_18O_CO2(aq)/CO2(g)") returns the
value for the log K at the current temperature using
expressions defined in NAMED_LOG_K data block;
LK_PHASE("Calcite") returns the value of log K
for calcite at the current temperature for the
dissociation reaction defined in the database or
input file. Values are "log10" values.
Example for Basic program:
10 PRINT "Log10 KCalcite: ", LK_PHASE("Calcite")
20 PRINT "Log10 KCaHCO3+: ", LK_SPECIES("CaHCO3+")
30 PRINT " 1000ln(alpha): ", LK_NAMED("Log_alpha_18O_CO2(aq)/CO2(g)")*LOG(10)*1000
NAMED_EXPRESSION--New keyword data block.
This data block was implemented to facilitate isotopic
calculations. It allows analytical expressions that
are functions of temperature to be defined. The purpose
is to separate the fractionation factors from the log
K, so that the fractionation factor or its temperature
dependence can be easily modified. The named expression
can be added to a log K for a species or phase by the
-add_logk identifier in SOLUTION_SPECIES
EXCHANGE_SPECIES, SURFACE_SPECIES, or PHASES data
block. Log K, Delta H, and analytical expressions for a
log K can be defined with identifiers -log_k, -delta_h,
and -analytical_expression as described in SOLUTION_SPECIES
in WRIR 99-4259. Fractionation factors are often defined
as 1000*ln(alpha). The identifier -ln_alpha1000 can be used
to enter data in this form. The analytical expression is the
same as defined in SOLUTION_SPECIES, but the result of the
expression is converted to log10(alpha) by dividing by
1000*ln(10) before it is summed into log K values.
NAMED_EXPRESSIONS
Log_K_calcite # CaCO3 + 2H3O+ = Ca+2 + 3H2O + CO2
log_k 8.201
delta_h -8.035 kcal
-analytic 292.29 0.015455 -24146.841 -94.16451 2248628.9
Log_alpha_18O_CO2(aq)/Calcite
-ln_alpha1000 3.8498 0.0 10.611e3 0.0 -1.8034e6
Log_alpha_13C_CO2(aq)/Calcite
-ln_alpha1000 2.72 0.0 0.0 0.0 -1.1877e6
------------------------------------------------------------
Added identifier -add_logk to SOLUTION_SPECIES
EXCHANGE_SPECIES, SURFACE_SPECIES, and PHASES data
block.
Allows a named expression to be added to the definition
of the log K for a species or phase. In the following
example, the log K for the phase Ca[14C][18O]3 is summed from
four parts, one defined with the log_k identifier and the
other three parts from expressions defined in NAMED_EXPRESSIONS.
The named expression is multiplied by the coefficient at the
end of the line before it is summed into the log K. A missing
coefficient is 1.0 by default.
PHASES
Ca[13C][18O]3
Ca[13C][18O]3 + 3CO2 + 2H3O+ = Ca+2 + 3H2O + 3CO[18O] + [13C]O2
log_k 0.903089986991 # 3*log10(2)
-add_logk Log_K_calcite 1.0
-add_logk Log_alpha_13C_CO2(aq)/Calcite 1.0
-add_logk Log_alpha_18O_CO2(aq)/Calcite 3.0
SOLUTION keyword:
At present, can only define isotopes in the units defined in ISOTOPES.
------------------------------------------------------------
Version 2.6:
------------------------------------------------------------
No new features.
------------------------------------------------------------
Version 2.5:
------------------------------------------------------------
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. For example, [Fe3], [13C],
and [N5] are now legal "element" names. All element
names without brackets must begin with a capital letter,
followed by zero or more lower case letters and underscores.
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.
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:
------------------------------------------------------------
------------------------------------------------------------
Added identifier -warnings to PRINT keyword.
An integer following -warnings gives the maximum number
of warnings to print into the output file. A negative
number allows all warnings to be printed.
Example: -warnings 20
------------------------------------------------------------
Changed the results of the function CELL_NO in Basic programs.
Function cell_no in Basic now prints a number equivalent
to -solution in SELECTED_OUTPUT data block. It gives the
solution number for initial solution calculations and the
solution being used in batch reaction calculations.
Result is the same as previous versions for ADVECTION or
TRANSPORT calculations.
------------------------------------------------------------
Version 2.3:
------------------------------------------------------------
DATABASE--New keyword data block
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.
LLNL_AQUEOUS_MODEL_PARAMETERS--New keyword data block
Added new keyword to make aqueous model similar to
EQ3/6 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.
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.
LLNL_AQUEOUS_MODEL_PARAMETERS
-temperatures
0.0100 25.0000 60.0000 100.0000
150.0000 200.0000 250.0000 300.0000
#debye huckel a (adh)
-dh_a
0.4939 0.5114 0.5465 0.5995
0.6855 0.7994 0.9593 1.2180
#debye huckel b (bdh)
-dh_b
0.3253 0.3288 0.3346 0.3421
0.3525 0.3639 0.3766 0.3925
-bdot
0.0394 0.0410 0.0438 0.0460
0.0470 0.0470 0.0340 0.0000
#cco2 (coefficients for the Drummond (1981) polynomial)
-co2_coefs
-1.0312 0.0012806
255.9 0.4445
-0.001606
------------------------------------------------------------
Added function SURF to Basic.
SURF("element", "surface") gives the amount of element
sorbed on "surface". "surface" should be the surface
name, not the surface-site name (that is, no underscore).
------------------------------------------------------------
Allow decimals in definition of secondary master species.
Some redox states do not average to integers,
for convenience in identifying them, decimal numbers
may be used within the parentheses that define the
redox state, example S(0.3) could be used in the
MASTER_SPECIES data block for the valence state of
aqueous species S6-2.
------------------------------------------------------------
Eliminate echo of input file in PRINT data block.
-echo_input T/F turns echoing on and off.
Default is true, initial value is true.
------------------------------------------------------------
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.
Example:
EQUILIBRIUM_PHASES
Dolomite 0.0 0.001 dis
------------------------------------------------------------
Version 2.2:
------------------------------------------------------------
Added function EDL to Basic.
EDL("element", "surface") gives the amount of
element in the diffuse layer for "surface", not
including sorbed species. "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.
"sigma" - surface charge density, C/m**2.
"psi" - potential at the surface, Volts.
"water" - mass of water in the diffuse layer, kg.
------------------------------------------------------------
End of Features not documented in WRIR 99-4259.
------------------------------------------------------------
------------------------------------------------------------
Revisions and bug fixes
------------------------------------------------------------
svn 675:
Added PRINT option to print the species that contribute
to alkalinity. Alkalinity distribution is printed in
the output file following the distribution of species.
Default at program startup is false.
PRINT
-alkalinity true
svn 655:
IAP and log K printed in Phase assemblage data
block were calculated from reactions rewritten to
new master species. Now the original data base
reaction is used to calculate IAP and log K.
Also fixed check that ensured all elements of
phase in are in solution before SI is calculated.
svn 631:
Fixed bug with alternate formula for equilibrium phase,
nothing happened if all other equations were satisfied
at beginning of reaction calculation.
svn 603:
Link gmp library statically.
svn 601:
Fixed statement on multiprecision.
svn 581:
Fixed bug in PhreeqcI that did not reinitialize
Chebyschev parameters leading to incorrect results
with Pitzer activity coefficients. Results were
correct on first run, but erroneous on subsequent
runs.
Added statement to identify multiprecision or
standard solver for inverse modeling.
svn 578:
Distribution changes. Fixed names in README file.
Modified Makefile to use specified version. Split
Linux and source distribution procedure.
Version 2.12 Date: Wed September 28, 2005
Executables and output files for Sun operating systems
are no longer provided.
Limited log activities of master species to be greater
than the smallest machine precision exponential number.
Avoids a matherr exception and allows trial of additional
parameter sets to attempt to solve the system of
equations.
Made aqueous activity coefficients the default activity
coefficients for exchange species when using the
Pitzer formulation. New option in EXCHANGE is
-pitzer_exchange_gammas T/F, default is true;
defining "false" sets exchange activity coefficients
to 1.0. Option has no effect for ion-association
model (non-Pitzer).
Edited phreeqc.dat to add -gamma expression for
CdX2 and PbX2.
Replaced O2(g) log K in phreeqc.dat and wateq4f.dat
with data from llnl.dat, which appears to be better.
Added multiplier format to REACTION increments, which
simplifies definition of multiple equal reaction increments.
REACTION
H2O 1
-36 3*-4 2*-.25 -0.19 4*-0.1 3*-0.05 moles
Added Pitzer activity formulation. Use pizer.dat database
to invoke the Pitzer model. Should have same capabilities
as ion-association model except explicit diffuse layer
calculation is not implemented with the Pitzer model.
Fixed bug in surface sites related to mineral and exchange
sites related to mineral. Did not have complete logic to
handle redox valence states in formula for species.
Modified to remove non standard usage of va_list.
Removed exchange master species from SYS("ex",..)
list of species. This species is fictive and should
not be included in the list.
Changed -redox in SOLUTION so that if one of the
redox states of a couple is not defined, then
redox defaults to pe.
Fixed buffer overrun in PhreeqcI with SOLUTION_SPREAD,
caused segmenatation fault with lines greater than
500 characters.
Added bigger string for some error messages to avoid
access violation in cvode.
Changed output_msg to warning_msg for combinations
of convergence parameters so that message would
be controlled by -warnings identifier.
Carriage returns are now stripped from Basic program
statements. Switching files from Windows to Unix sometimes
leaves extra carriage returns at the ends of lines, which
caused a syntax error for some Basic commands.
Two bugs were fixed in inverse modeling. (1) Potential
models are now checked for all equality and inequality
constraints. Previously some constraints were not checked.
(2) One loop of cl1 did not include the last row when
checking for the pivot element. Also printing of headers
was slightly modified for inverse modeling.
A new multiple precision version of cl1 was develeped by
using the Gnu Multiple Precision package (gmp). Calculations
are carried out to about 30 significant digits. The mp
version may help in some situations where roundoff errors are
a problem, but it is still possible that roundoff errors will
cause cl1mp to fail to find a solution to the optimization
problem. The mp version has the following options in
INVERSE_MODELING:
-multiple_precision T/F--causes the mp version
to be used in inverse modeling calculations.
-mp_tolerance 1e-12--tolerance for mp version of
cl1. As in cl1, numbers less than the
tolerance are considered to be zero.
1e-12 is the default.
-censor_mp 1e-20--as calculations occur in the
linear equation array, elements less
than this value are set to zero. Default
is 1e-20. A value of 0.0 causes no
censoring to occur.
Version 2.11 Date: Mon February 7, 2005
Fixed error in selected output file with mixing reaction.
MIX number was written to two columns, should be one.
Fixed memory leak with PAD function.
New database minteq.v4.dat has been translated from version
4.02 of MINTEQA2. An older version of the MINTEQA2 database is
retained in file minteq.dat.
Switched version control to subversion. Simplified,
reorganized makefiles.
Fixed bug with PRINT; -warnings n. Use of this option
generally eliminated all warning messages instead of
all messages after the nth. Default number of warning
messages printed in now 100 per simulation.
Fixed memory leaks in cvode that caused phreeqci to crash.
Now uses PHRQ_malloc in case of other memory leaks. Also
fixed potential memory error with PAD Basic function.
Saturation index phases that included water had wrong
value if distribution of species, exchange, or surface
not written also.
Fixed error message in cvode, if max iterations exceeded
the error caused a segmentation fault.
Made printing of parameter combination message a warning
message so that it could be turned off.
Version 2.10 Date: Tue November 2, 2004
Rearranged i/o for PHREEQC and reorganized driver
subroutine. The object of these changes is to make
the program more functional as a module for other
programs (PHAST) and eventually to produce a callable
C and Fortran module.
Fixed a problem with surface related to a phase, when
phase was not part of system (for example, Fe(OH)3a when
there is no iron in system.
Added convergence parameter set that requires mineral
transfers to produce positive concentrations in the
event that negative concentrations have been produced in
the prior Newton-Raphson iteration. (Fluorite example).
Fixed bug with kinetics formulas; did not account for
stoichiometric coefficient correctly when using
phase names. Generalized to allow multiple phase
names in the -formula definition.
Version 2.9 Date: Wed September 15, 2004
In inverse modeling, program terminates if sum of initial
solutions and phases is > 32.
Fixed bug with isotopes. Log activity estimate after initial
solution calculation was inf under some conditions. An initial
surface calculation failed when using D.
Changed saturation index print out to use reaction and log K
defined in PHASES definition. Previously, reaction could be
rewritten to predominant redox species.
Fixed incorrect print of elapsed time for kinetics in advection.
Added phrqproto.h prototype file and phrqtype.h for
switching compilation to long double.
Fixed incorrect printout of kinetics delta moles with
advection.
Added convergence parameter set that skips mineral
equations for first 5 iterations.
Added entity_exists for module.
Tried to fix bug with mix index incorrect (-2) for
mixing with kinetics.
Added new keyword COPY that allows a data entity
to be copied from one index to a new index
or to a range of indicies. Format is
COPY keyword index index_start[-index_end]
where keyword is one of the following:
SOLUTION
EQUILIBRIUM_PHASES
EXCHANGE
GAS_PHASE
KINETICS
MIX
REACTION
REACTION_TEMPERATURE
SOLID_SOLUTION
SURFACE
Numerical method had a bug with ionic strength, if
mass of water was not approximately 1. Routine
revise_guesses did not divide by the mass of
water. Also changed check in routine molalities
to check by molality, not moles.
Fixed a null pointer when surface was related to a
mineral and mineral was not in database.
Added new Basic functions
i = INSTR(a$, b$) sets i to the character position of
string b$ in a$, 0 in not found.
b$ = LTRIM(a$) trims white space from beginning of
string a$ and stores result in b$.
b$ = RTRIM(a$) trims white space from end of string
a$ and stores result in b$.
b$ = LTRIM(a$) trims white space from beginning and
end of string a$ and stores result in b$.
Added new Basic function SYS that calculates the to
total amount of an element in all phases (solution,
equilibrium_phases, surfaces, exchangers, solid solutions,
and gas phase). KINETIC reactions are not included.
The function has two forms: (1) one element name as an
argument (variable names are user specified)
10 t = SYS("As")
the function will return the total arsenic in the system.
(2) 5 argumens
10 t = SYS("As", count_species, names$, types$, moles)
will return the total arsenic in the system to tot; count_species--
the number of species that contain arsenic, including
solution, equilibrium_phases, surfaces, exchangers, solid solutions,
and gas phase species; names$--a character array that has the
name of each species; type$--a character array that specifies the
type of phase for the species, aq, equi, surf, ex, s_s, gas, diff.
Diff refers to the amount of the element in the diffuse layer of
a surface when the explicit diffuse layer calculation is used;
moles--an array containing the number of moles of the element in
the species. The sum of moles(i) is equal to tot.
SYS has several special arguments for the form
SYS("arg", count, names$, types$, values)
arg is one of the options listed below.
count is a single numeric value and is the number of elements
in the following arrays.
name$ is an array of string values.
type$ is an array of string values.
values is an array of numeric values.
Values of arg:
elt_name returns total number of moles of element in system.
count is the number of species for the element in
the system, including aqueous, exchange, surface,
equilibrium_phase, solid solution component, and
gas phase "species".
Arrays are filled for each "species"; values are moles.
"elements" returns total number of moles of all elements,
valence states, exchangers, and surfaces.
count is number of elements, valence states,
exchangers, and surfaces.
Arrays are filled for each element, valence state,
exchanger, and surface; values are moles.
"phases" returns maximum saturation index of all phases.
count is number of phases in system.
Arrays are filled for each phase; values are
saturation indexes.
"aq" returns sum of moles of all aqueous species.
count is number of aqueous species in system.
Arrays are filled with each aqueous species;
values are moles.
"ex" returns sum of moles of all exchange species.
count is number of exchange species in system.
Arrays are filled with each exchange species;
values are moles.
"surf" returns sum of moles of all surface species.
count is number of surface species in system.
Arrays are filled with each surface species;
values are moles.
"surf" returns sum of moles of all solid solution components.
count is number of solid solution components in system.
Arrays are filled with each solid solution component;
values are moles.
"gas" returns sum of moles of all gas components.
count is number of gas components in system.
Arrays are filled with each gas component;
values are moles.
Added new Basic function, DESCRIPTION, that has the value
defined for the description field of the SOLUTION keyword line.
Added alternative ordinary differential equation solver
called CVODE, a set of C routines from the Lawrence
Livermore National Labs. CVODE is part of the SUNDIALS
package. CVODE is used in place of the Runge Kutta method
when "-cvode true" is used within a KINETICS data block.
KINETICS
-cvode true
Fixed error in SOLUTION_SPREAD, defining -redox
did not set the default redox for the solutions
that were defined; pe was always used as default.
Modified code to allocate space differently for
pp_assemblage, exchange, surface, gas_phase,
kinetics, and s_s_assemblage. Enough space is allocated
at beginning of distribute_initial_conditions.
May speed up phast initialization and make better
use of available memory.
Changed gfw of water to 18 if isotopes of water are
included. Solvent is [1H]2[16O].
Fixed a bug in surface integration where order of ions
in the list of g's was incorrect.
Pyrite rate was not 0 if supersaturated in phreeqc.dat
and wateq4f.dat
Segmentation error if a surface species was not
defined with an equation that contained another surface
species. In this case, the surface master species
had been redefined to be an aqueous species (SOLUTION_SPECIES).
Version 2.8 Date: Tue April 15, 2003
Updated arsenic data in wateq4f.dat to be consistent with
Archer and Nordstrom (2002).
Revised Basic interpreter to allow lines of any length
and character strings of any length.
Renumbering basic statement that included the function
SURF in PhreeqcI caused SURF to be omitted and generated
a syntax error. SURF and other functions had not been
implemented in PhreeqcI.
Fixed a bug in the Basic Interpreter. If elements of
a dimensioned variable (character or number) were used on
both sides of an equation, the result was erroneously
stored in the last element of the variable used on the
right-hand side instead of the element specified on the
left-hand side.
Using comma in some fields caused an infinite loop.
Fixed bug with SOLUTION_SPREAD, Phreeqc was not
calculating solution numbers for solution_spread
solutions without solution numbers.
Fixed bug with stagnant zone calculations. If solutions
not defined for stagnant cells, PhreeqcI crashed.
Added new state for calculations, PHAST. Previously
phast used the state TRANSPORT, which caused some
erroneous results with temperature when TRANSPORT was
used in the PHREEQC part of the calculation.
Trying to define dump file in TRANSPORT caused a file
opening error. Fixed logic so now can open a file
and the name can include blanks.
Version 2.7 Date: Fri February 14, 2003
Initialized gfw in elements structure.
Fixed bug where "time" would be wrong for initial
solution calculation. Needed to initialize
rate variables for PhreeqcI.
Added print of simulation number to error file for
phreeqci
Limited printing of cell numbers to output file in advection
calculations. Cell numbers only printed if results
for cell are going to be printed.
PhreeqcI captured status messages for kinetics, which
made a very large error file in some cases and
a long wait to view the output file in PhreeqcI.
Now PhreeqcI does not capture these intermediate
status messages.
Removed old code related to redirecting error file
Corrected error in transport where wrong time step was used
for integration.
Changes to speed up transport algorithm.
Allow file names with spaces in selected_output file name and
dump_file name.
Modifications to work with RC1 phast log file.
Allow any characters in square brackets for element name.
- and + and perhaps others caused problems before.
Fixed log molality of water in species printout, was
equal to log activity of water. Also fixed
basic function for LM.
Changed solid solution prints to print 0 if solid solution
is not present.
Fixed bug if no rate name was defined before options
in RATES.
Fixed warning on Mac compilation in fpunchf.
Fixed bug if isotopes were used but H and O isotopes
were not defined.
Fixed bug where special initial solution calculations
were done at later calculation stages.
Needed to set initial_solution_isotopes = FALSE;
Fixed problem in C++ if structure name is same as member name.
logk member of logk structure was renamed to log_k.
Added identifier -add_constant to PHASES, EXCHANGE_SPECIES,
SOLUTION_SPECIES, and SURFACE_SPECIES.
-add_constant -0.301
log K is augmented by the specified constant.
Added punch_isotopes and punch_calculate_values to allow
printing isotope ratios and any CALCULATE_VALUES result.
Added KEYWORDS:
ISOTOPES
Element
-isotope isotope_name units standard_ratio
-total_is_major T/F (OPTION IS DISABLED!!)
CALCULATE_VALUES
Name
-start
Basic statements, must have SAVE
-end
ISOTOPE_RATIOS (for printing)
Name=Calculate_values_name Isotope_name
ISOTOPE_ALPHAS (for printing)
Name=Calculate_values_name Named_logk=named_expression_name
Basic functions:
calc_value("calc_value_name") evaluates a definition of CALCULATE_VALUES
lk_named("name") log10(K) of definition in NAMED_EXPRESSIONS
lk_phase("name") log10(K) of definition in PHASES
lk_species("name") log10(K) of definition in (SOLUTION, EXCHANGE, SURFACE)_SPECIES
sum_gas("template","element") Sum of element in gases with specified template
template="{C,[13C],[14C]}{O,[18O]}2" includes all CO2 gases
sum_species("template","element") Sum of element in aqueous, exchange, and surface species with
specified template
sum_s_s("s_s_name","element") Sum of element in a specified solid solution
PRINT keyword:
-initial_isotopes T/F
-isotope_ratios T/F
-isotope_alphas T/F
-censor_species 1e-8 # Omits print of species if less than relative criterion
SELECTED_OUTPUT keyword:
-calculate_values name1 name2 ...
-isotopes minor_isotope1 minor_isotope2 ....
Added functions LK_SPECIES, LK_NAMED, LK_PHASE for Basic
interpreter. LK_SPECIES("CaHCO3+") returns the
log k for the association reaction for the ion pair
CaHCO3+ at the current temperature. The log K is
for the reaction as defined in the database or
input file. Similarly,
LK_NAMED("Log_alpha_18O_CO2(aq)/CO2(g)") returns the
value for the log K at the current temperature using
expressions defined in NAMED_LOG_K data block;
LK_PHASE("Calcite") returns the value of log K
for calcite at the current temperature for the
dissociation reaction defined in the database or
input file. Values are "log10" values.
Example for Basic program:
10 PRINT "Log10 KCalcite: ", LK_PHASE("Calcite")
20 PRINT "Log10 KCaHCO3+: ", LK_SPECIES("CaHCO3+")
30 PRINT " 1000ln(alpha): ", LK_NAMED("Log_alpha_18O_CO2(aq)/CO2(g)")*LOG(10)*1000
Added NAMED_EXPRESSIONS data block. This data block was
implemented to facilitate isotopic calculations.
It allows analytical expressions that are functions
of temperature to be defined. The purpose is to
separate the fractionation factors from the log K,
so that the fractionation factor or its temperature
dependence can be easily modified. The named
expression can be added to a log K for a species
or phase by the -add_logk identifier in SOLUTION_SPECIES
EXCHANGE_SPECIES, SURFACE_SPECIES, or PHASES data
block.
Version 2.6 Date: Mon April 22, 2002
PhreeqcI released.
All selected_output is routed through a single routine.
Allow "_" inside square brackets, [A_bcd].
Fixed bug match_elts_in_species, check for "e-" was wrong.
Modified minteq.dat to put CuS4S5-3, Cu(S4)2-3 in
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.
Made revisions hopefully to improve SOLID_SOLUTIONS
convergence with small numbers of moles of
solids.
Made changes related to dump file and PhreeqcI.
Iterations now sums iterations in all kinetics calculations
Fixed bug with LA("H2O"), which was returning natural log
of activity of water.
Version 2.5 Date: Mon October 1, 2001
In llnl.dat, fixed sign errors in RRE (rare earth elements)
for some redox reactions and removed some redundant
species, generally ReeO2- was retained and Ree(OH)4-
was removed.
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.
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.
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.
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".
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.
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.
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.
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: Date: Fri June 15, 2001
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: Date: Mon June 4, 2001
Fixed spreadsheet bugs with isotopes.
Version 2.4: Date: Fri June 1, 2001
Added structure for spreadsheet for use by PhreeqcI.
Isotope value initialized incorrectly if only an -uncertainty was
defined in SOLUTION_SPREAD.
Fixed segmentation violation when primary and secondary master
species were defined improperly.
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.
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.
Added surface complexation constants for Fe+2; two complexes for
weak sites and one complex for strong sites. phreeqc.dat
and wateq4f.dat modified.
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.
Fixed a bug when rates were equal within tolerance, but negative
concentrations occurred because of small initial
concentrations.
Added -warnings to PRINT keyword for specification of maximum
number of warnings to print. Negative number allows
all warnings to be printed.
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.
Kinetics time is halved for advective part of reaction in
transport; time incorrectly accounted for before.
-punch_ identifiers printed -1 instead of the correct solution
number for batch-reaction calculations.
-high_precision is no longer reset to false with every
SELECTED_OUTPUT data block.
SELECTED_OUTPUT file name stored for use by PhreeqcI.
Alkalinity for NH3 corrected to 1.0 in llnl.dat.
Fixed bug with USER_PRINT of kinetics. Did not find correct
kinetics information in some cases.
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: Date: Tue January 2, 2001
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.
Fixed bug in SOLUTION_SPREAD. If first heading in
the spread-sheet input was an identifier--pH,
pe, units, etc--then the headings were interpreted
as an identifier and bad things happened.
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.
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.
LLNL_AQUEOUS_MODEL_PARAMETERS
-temperatures
0.0100 25.0000 60.0000 100.0000
150.0000 200.0000 250.0000 300.0000
#debye huckel a (adh)
-dh_a
0.4939 0.5114 0.5465 0.5995
0.6855 0.7994 0.9593 1.2180
#debye huckel b (bdh)
-dh_b
0.3253 0.3288 0.3346 0.3421
0.3525 0.3639 0.3766 0.3925
-bdot
0.0394 0.0410 0.0438 0.0460
0.0470 0.0470 0.0340 0.0000
#cco2 (coefficients for the Drummond (1981) polynomial)
-co2_coefs
-1.0312 0.0012806
255.9 0.4445
-0.001606
Fixed bug in basic interpreter. A number like "..524" would
cause an infinite loop.
Added function SURF to Basic.
SURF("element", "surface") gives the amount of
element sorbed on "surface". "surface"
should be the surface name, not the
surface-site name (that is, no underscore).
Fixed option to "runge_kutta" from "runge-kutta" to match
documentation for KINETICS.
Fixed UO2+2 and Mn+2 reaction stoichiometry for Hfo surface complexation
in wateq4f.dat.
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.
Example:
EQUILIBRIUM_PHASES
Dolomite 0.0 0.001 dis
R-K integration failed when only the final rate generated
negative concentrations.
Allow decimals in definition of secondary master species, for
example S(0.3).
Fixed bug if description was more than about 85 characters;
now allows about 400 characters.
Fixed bug for surface/exchange sites related to phases. Was
checking internal copies of surfaces/exchange with negative
numbers.
Fixed bug in quick prep that did not set the correct pointer
for gas phases.
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.
Added option to eliminate echo of input file in PRINT
data block. -echo_input T/F turns echoing on
and off. Default is on.
Release 2.2: Date: Wed March 1, 2000
Fixed bug in MIX if no solutions are defined.
Changed printout for surface.
Only gives net surface charge for diffuse layer
calculation.
Prints correct value for the surface charge and
surface charge density for diffuse-layer
calculation.
Added function EDL to Basic.
EDL("element", "surface") gives the amount of
element in the diffuse layer for "surface".
not including sorbed species. "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.
"sigma" - surface charge density, C/m**2.
"psi" - potential at the surface, Volts.
"water" - mass of water in the diffuse layer, kg.
Changed distribution to be more consistent with other USGS
water-resources applications.
Release 2.1: Date: Wed January 19, 2000
Added additional #ifdef's for PhreeqcI.
Fixed problem with formats for USER_PUNCH and
others with Microsoft C++ 3 digit
exponents.
Initial Release 2.0: Date: Wed December 15, 1999
Version: C_54 = Version 2.0
Reference in New Issue View Git Blame Copy Permalink