************************************************************ ************************************************************ * Features not documented in WRIR 99-4259. * ************************************************************ ************************************************************ ------------------------------------------------------------ Version 2.14.3: November 17, 2007 ------------------------------------------------------------ -------- svn 2312 -------- Added new option to PITZER datablock, use_etheta t/f. If true, the nonsymmetric mixing terms--cation/cation and anion/anion of different charge--are included; if false these terms are excluded from all equations. Default is true. PITZER -use_etheta true -------- svn 2270 -------- Added additional parameters Pitzer activity formulation for neutral species, MU and ETA. MU applies to nnn, nnn', nn'n'', nna, nn'a, nnc, nn'c interactions, where n, n', and n'' are neutral species, a is an anion and c is a cation. ETA applies to ncc' and naa' interactions. Also modified LAMDA for the special case of nn interactions (coefficients in osmotic and ln equations are different than other interaction types). Source of equations is Clegg and Whitfield, 1991, Activity coefficients in natural waters, Chapter 6, in Pitzer, K.S. (Ed.) Activity Coefficients in Electrolyte Solutions, 2nd Ed. CRC Press, Boca Raton. Removal of the 6 coefficient in last two terms of eq 35 and 36 (p. 2404) per Cleg and Whitfield, 1995, Geochimica et Cosmochemica Acta, v. 59, no. 12, pp 2403-2421. Order of species in definitions should not matter. PITZER -MU CO2 CO2 CO2 ? # nnn CO2 CO2 NH3 ? # nnn' CO2 B(OH)3 NH3 ? # nn'n'' CO2 CO2 Ca+2 ? # nnc CO2 CO2 Cl- ? # nna CO2 NH3 Ca+2 ? # nn'c CO2 NH3 Cl- ? # nn'a -ETA CO2 Ca+2 Mg+2 ? # ncc' CO2 Cl- SO4-2 ? # naa' As with all other Pitzer parameters, a five-term expression for temperature dependence is available: 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. A0 through A4 are defined in order. Any undefined values are assumed to be zero. -MU CO2 CO2 CO2 ? ? ? ? ? # nnn ------------------------------------------------------------ Version 2.14.2: September 17, 2007 ------------------------------------------------------------ Fixed logic of memory checking for PhreeqcI. This serious bug makes versions 2.14.0 and 2.14.1 unusable. ------------------------------------------------------------ Version 2.14.1: September 5, 2007 ------------------------------------------------------------ No new features. ------------------------------------------------------------ Version 2.14.0: August 30, 2007 ------------------------------------------------------------ No new features. ------------------------------------------------------------ Version 2.13.3: February 15, 2007 ------------------------------------------------------------ No new features. ------------------------------------------------------------ Version 2.13.2: February 1, 2007 ------------------------------------------------------------ No new features. ------------------------------------------------------------ Version 2.13.1: January 16, 2007 ------------------------------------------------------------ No new features. ------------------------------------------------------------ Version 2.13.0: November 3, 2006 ------------------------------------------------------------ -------- svn 1368 -------- (1) Added multicomponent diffusion (MCD) to transport capabilities. MCD allows different tracer diffusion coefficients for species, but calculates charge balanced transport. In the example, MCD is specified to be true, default tracer diffusion coefficient for species (Dw) is 1e-9, porosity is set to 0.3, porosity limit is set to 0.05, and an exponent of porosity (n) is set to 1.0. Effective diffusion coefficient is defined by the equation: De = Dw * porosity^n. Diffusion stops when the porosity falls below the porosity limit. TRANSPORT -multi_d true 1e-9 0.3 0.05 1.0 Added tracer diffusion coefficient to SOLUTION_SPECIES definitions, -dw identifier. SOLUTION_SPECIES H+ = H+ log_k 0.0 -gamma 9.0 0.0 -dw 9.31e-9 (2) Added phreeqd.dat database with diffusion coefficients (-dw) defined for aqueous species in database directory. (3) Added BASIC functions to obtain and modify the porosity in a cell. The functions can be used in BASIC programs defined with keyword RATES, USER_PRINT, USER_PUNCH and USER_GRAPH. get_por(cell_no) # returns the porosity in cell # 'cell_no' change_por(0.21, cell_no) # porosity of cell 'cell_no' # becomes 0.21 (4) Mobile surface and Donnan option in SURFACE. Mobile surfaces are meant for modeling transport of colloids. Only surfaces with a diffuse double layer can be transported (the ensemble must be electrically neutral). Surfaces related to equilibrium-phases and kinetics cannot be transported. Example 1: Use donnan assumption to calculate the explicit composition of the diffuse layer of surfaces. Thickness of the diffuse layer is defined to be 1e-7 meters. (Default thickness is 1e-8 meters.) Hfo (both sites Hfo_w and Hfo_s) is a surface that is transported with advection and dispersion. The diffusion coefficient of 1e-13 m^2/s is used with option -multi_d true in TRANSPORT. Sfo is an immobile surface (Dw = 0). SURFACE -donnan 1e-7 Hfo_w 97.5e-5 600 88e-3 Dw 1e-13 Hfo_s 2.5e-5 Sfo_w 97.5e-5 600 88e-3 Dw 0 Sfo_s 2.5e-5 Example 2: Define Donnan calculation information. Thickness of the diffuse layer is 1e-8 meters, and relative viscosity is 0.5. Relative viscosity only applies to multicomponent diffusion of solutes in the diffuse layer. (Default viscosity is 1.0.) SURFACE -donnan 1e-8 viscosity 0.5 Example 3: Define Donnan calculation information. Thickness of the diffuse layer is 1.5 Debye lengths, maximum fraction of water that can be in the diffuse layer is 0.9. (Default thickness in Debye lengths is 1, default limit is 0.8.) SURFACE -donnan debye_lengths 1.5 limit_ddl 0.9 When option '-only_counter_ions' is used in conjunction with with '-donnan', all the co-ions (with the same sign of charge as the surface) will be excluded from the DDL, and will be given a concentration of (near) zero in the DDL. (5) Added special BASIC function to change the diffusion coefficient of (part of) a SURFACE, and hence to change the status from mobile to immobile or immobile to mobile. Example 1: take a fraction 0.2 of 'Hfo', rename it to 'Sorbedhfo', with a diffusion coefficient of 0, in cell 'cell_no' USER_PRINT 10 print 'Changing surface in cell ', cell_no 20 change_surf("Hfo", 0.2, "Sorbedhfo", 0, cell_no) Example 2: change the diffusion coefficient of 'Hfo' to 1e-12 m2/s in cell 'cell_no' 10 change_surf("Hfo", 1, "Hfo", 1e-12, cell_no) This function can be used in BASIC programs defined with keywords RATES, USER_PRINT, USER_PUNCH and USER_GRAPH. For correct operation of 'change_surf', the surface components must have been defined with the same surface species (the association constants may differ) and the same diffuse layer thickness or Debye length. The surfaces will be adapted at the end of a calculation step. The result of change_surf does not show up in print or punch results of that step, but the reformatting is effective in the next timestep calculations. -------- svn 1337 -------- Added -add_logk to NAMED_EXPRESSIONS keyword. NAMED_EXPRESSIONS Log_alpha_14C_CO3-2/CO2(aq) -add_logk Log_alpha_14C_CO3-2/CO2(g) 1 -add_logk Log_alpha_14C_CO2(aq)/CO2(g) -1 -------- svn 1281 -------- Added new option to PITZER data block to allow definition of alpha1 and alpha2 for specific electrolytes. Entries are following -ALPHAS are Ion1, Ion2, alpha1, alpha2. Both alpha1 and alpha2 should be defined. Default is 0.0 for undefined values following Ion1 and Ion2. Example: PITZER -ALPHAS # # Defaults for ion valences in salts # # 1-N (only B1): alpha1 = 2 # 2-2: alpha1 = 1.4 alpha2 = 12.0 # 3-2, 4-2: alpha1 = 2 alpha2 = 50 # #Ion1 Ion2 Alpha1 Alpha2 Fe+2 Cl- 2 1 Fe+2 SO4-2 1.559 5.268 -------- svn 1279 -------- Added new Basic function OSMOTIC that returns the osmotic coefficient if the Pitzer model (PITZER keyword data block) is used or 0.0 if the ion-association model is used. Example: USER_PRINT 10 PRINT "Osmotic coefficient: ", OSMOTIC -------- svn 1245 -------- Enabled redox in Pitzer model with option in PITZER keyword. Typically, the option will be included in the pitzer database file. Example: PITZER -redox TRUE The default database for the Pitzer model does not contain any redox definitions and the default value for the option is FALSE. At a minimum, species O2 and H2 must be defined in the database or input file to allow redox calculations. -------- svn 1207 -------- Added option to force an equilibrium phase to be included in the equality constraints. Normally, the SIs of equilibrium phases are optimized to be negative and the sum of SIs is minimized. If -force_equality is used, then the phase must reach its target SI or the calculation fails with an error. Example: EQUILIBRIUM_PHASES Fix_pH -7 NaOH -force_equality Calcite 0 Dolomite 0 CO2 -3.5 One example of using the new option would be to ensure that a target pH is attained, as in the example above. -------- svn 1179 -------- New option (-sites_units density) allows alternative units (sites/nm^2) for definition of number of sites for a surface. This approach requires better consistency among the parameters as both the number of sites and the surface area are based on the mass. It makes more sense than the default, which requires the number of sites (first numeric item in a line) to be defined in units of moles, independently of the number of grams of sorbent. Example: SURFACE 1 -sites_units DENSITY SurfOH 2.6 600. 1.0 SurfaOH 2.6 30. 2.0 Explanation: In this example, Surf has a site density of 2.6 sites per nanometer squared, a specific area of 600 meters squared per gram, and a mass of 1 gram. Surfa has a site density of 2.6 sites per nanometer squared, a specific area of 30 meters squared per gram, and mass of 2 grams. -------- svn 1096 -------- Allows solids and gases in the equations for PHASES. This capability simplifies the definitions for gas and solid isotopic components. Solids must be identified with "(s)" and gases with "(g)". The first entity on the left- hand-side of the equation must be the stoichiometric formula of the solid or gas component being defined, optionally with (g) or (s). In turn gases and solids included in the equation must be defined with reactions that ultimately allow the defined species (C[18O]2(g) in this case) in terms of aqueous species. Example: PHASES C[18O]2(g) C[18O]2(g) + CO2(g) = 2CO[18O](g) log_k 0.602059991327962396 # log10(4) -------- svn 1092 -------- CD_MUSIC sorption model has been implemented. Still missing logic for surfaces related to equilibrium- phases and kinetics. Has explicit calculation of diffuse layer composition with Donnan assumption. Old diffuse-layer calculation will not be implemented. Example: SURFACE Goe_uniOH .000552 96.387 1 -capacitance 1.1 5 Goe_triO .000432 -cd_music -donnan Explanation: 1.1 5 are capacitances for the cd-music model for 0-1 and 1-2 planes, respectively. -cd_music specifies that the surface is a cd-music surface. -donnan optionally calculates the diffuse layer composition with the Donnan model. Example: SURFACE_SPECIES Goe_uniOH-0.5 + H+ + AsO4-3 = Goe_uniOAsO3-2.5 + H2O log_k 20.1 # eq 7 K1, Kin1 -cd_music -1 -6 0 0.25 5 Explanation: -cd_music gives the charge contribution of the surface species to the three planes. Plane 0 is -1 + 0.25*5; Plane 1 is -6 + (1-0.25)*5; Plane 2 (or d) is 0. -------- 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. Example: 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, moles. Example: 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 (moles) sum_s_s("s_s_name","element") Sum of element in a specified solid solution (moles) 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 * ************************************************************ ************************************************************ ------------------------------------------------------------ Version 2.14.3: November 17, 2007 ------------------------------------------------------------ -------- svn 2386 -------- Fixed bug in routine find_Jstag. Incorrect index (cell_no) caused segmentation violation in rare instances. -------- svn 2312 -------- Added new option to PITZER datablock, use_etheta t/f. -------- svn 2279 -------- Error lost 200 moles of mineral. Should only be a problem in some cases where moles of mineral is greater than 200. -------- svn 2270 -------- Added additional parameters Pitzer activity formulation for neutral species, MU and ETA. -------- svn 2269 -------- Fixed buffer overrun in SOLUTION_SPREAD when pasting. -------- svn 2268 -------- Fixed error in prep.c where realloc was called instead of PHRQ_realloc, which eliminated a memory leak. ------------------------------------------------------------ Version 2.14.2: September 17, 2007 ------------------------------------------------------------ -------- svn 2267 -------- Fixed logic of memory checking for PhreeqcI. This serious bug makes versions 2.14.0 and 2.14.1 unusable. ------------------------------------------------------------ Version 2.14.1: September 5, 2007 ------------------------------------------------------------ -------- svn 2219 -------- Updated transport.c to adjust transport in diffuse layer to be charge balanced for MCD calculation. ------------------------------------------------------------ Version 2.14.0: August 30, 2007 ------------------------------------------------------------ ------------- svn 2203-2204 ------------- Revised logic for using phqalloc memory checker. Compiler option USE_PHRQ_ALLOC now turns on memory checker. If USE_PHRQ_ALLOC is defined and NDEBUG is not defined, file name and line number are also used in memory checking. Model.c now uses same compile options as all other files. -------- svn 2199 -------- Initialized variables that caused problems when rerunning simulations in PfW and PhreeqcI. -------- svn 2138 -------- Fixed bugs in MCD calculation related to saving solutions after initialization. -------- svn 2055 -------- Profiled and optimized of code. Automatic string in Basic factor saves many mallocs. Reordered to minimize call to strcmp_nocase in basicsubs.c and xsolution_save. Minimized mallocs for solver. -------- svn 2051 -------- Fixed bugs in MCD calculation. -------- svn 2040 -------- Fixed warnings for type-punned with new gcc. Reverted to 2.12 for infilling solutions for transport. Only solutions are used, not additional reactants for solution 0 and n+1. Added digits to printout of REACTION stoichiometry. -------- svn 1852 -------- Fixed error in CDMUSIC surface related to a phase. Stoichiometry of H was calculated incorrectly. -------- svn 1837 -------- Initialize flag for MCD calculation. PhreeqcI would do MCD calculation after TRANSPORT was redefined not to do MCD calculation. ------------------------------------------------------------ Version 2.13.2: February 1, 2007 ------------------------------------------------------------ -------- svn 1700 -------- Fixed bug with redox elements in multi_diffusion stagnant zones. -------- svn 1683 -------- Worked on convergence problems when optimizer fails to find a solution. Censored values greater than 1e8. -------- svn 1637 -------- Fixed bug with dissolve only in ADVECTION calculations. -------- svn 1629 -------- Fixed bug with redox elements in multi_diffusions. Added Phreeqc For Windows changes from 2.13.1. ------------------------------------------------------------ Version 2.13.1: January 16, 2007 ------------------------------------------------------------ -------- svn 1600 -------- Fixed logic error that required rebuilding aqueous model when not necessary. Now runs faster (sometimes 3X) than version 2.13.0. -------- svn 1590 -------- Removed porosity from one statement to eliminate oscillations in multicomponent diffusion calculation. -------- svn 1558 -------- Dissolve-only option did not work correctly for stagnant cells in TRANSPORT calculations. The moles of equilibrium phases, kinetics, gas_phase, and solid solutions were not initialized at the beginning of each transport step. Thus, the printed values for delta moles for the step in the output and punch file were incorrect for stagnant-zone cells. "Dissolve only" was always tested relative to the number of moles initially in the cells, not the amount remaining at a given time step. -------- svn 1485 -------- Pitzer version with gas_phase did not work. Added gas_phase and cd_music to numerical derivative routine. svn 1368: (1) Added multicomponent diffusion in transport and SOLUTION_SPECIES. (2) Added BASIC functions to obtain and modify the porosity in a cell. (3) Aded mobile surface and Donnan option in SURFACE. (4) Added special BASIC function to change the diffusion coefficient of a SURFACE, and hence to change the status from mobile to immobile or immobile to mobile. svn 1337: Added -add_logk to NAMED_EXPRESSIONS keyword. svn 1306: Revised printing of distribution of species, pure phase assemblages, and solid solutions to use longer fields for names. More revisions to logic for using gases and solids in equations for phases. Revised logic for solid solutions with small (1e-25) amounts of component. svn 1282: Fixed bug when gas phase had no gas components. Looked the same as not having a gas phase at all. svn 1245: Enabled redox in Pitzer model with option in PITZER keyword. Typically, the option will be included in the pitzer database file. PITZER -redox TRUE The default database for the Pitzer model does not contain any redox definitions and the default value for the option is FALSE. At a minimum, species O2 and H2 must be defined in the database or input file to allow redox calculations. svn 1179: New option (-sites_units density) allows alternative units (sites/nm^2) for definition of number of sites for a surface. This approach requires better consistency among the parameters as both the number of sites and the surface area are based on the mass. It makes more sense than the default, which requires the number of sites (first numeric item in a line) to be defined in units of moles, independently of the number of grams of sorbent. SURFACE 1 -sites DENSITY SurfOH 2.6 600. 1.0 SurfaOH 2.6 30. 2.0 In this example, Surf has a site density of 2.6 sites per nanometer squared, a specific area of 600 meters squared per gram, and a mass of 1 gram. Surfa has a site density of 2.6 sites per nanometer squared, a specific area of 30 meters squared per gram, and mass of 2 grams. svn 1128: Fixed bug with value of time printed to selected output file when using cvode. Value printed was an intermediate integration time step, not time at end of integration. svn 1096: Allows solids and gases in the equations for PHASES. This capability simplifies the definitions for gas and solid isotopic components. Solids must be identified with "(s)" and gases with "(g)". The first entity on the left- hand-side of the equation must be the stoichiometric formula of the solid of gas component being defined, optionally with (g) or (s). In turn gases and solids included in the equation must be defined with reactions that ultimately allow the defined species (C[18O]2(g) in this case) in terms of aqueous species. C[18O]2(g) C[18O]2(g) + CO2(g) = 2CO[18O](g) log_k 0.602059991327962396 # log10(4) svn 1092: CD_MUSIC sorption model has been implemented. Still missing logic for surfaces related to equilibrium- phases and kinetics. Has explicit calculation of diffuse layer composition with Donnan assumption. Old diffuse-layer calculation will not be implemented. Example: SURFACE Goe_uniOH .000552 96.387 1 -capacitance 1.1 5 Goe_triO .000432 -cd_music -donnan 1.1 5 are capacitances for the cd-music model for 0-1 and 1-2 planes, respectively. -cd_music specifies that the surface is a cd-music surface. -donnan optionally calculates the diffuse layer composition with the Donnan model. SURFACE_SPECIES Goe_uniOH-0.5 + H+ + AsO4-3 = Goe_uniOAsO3-2.5 + H2O log_k 20.1 # eq 7 K1, Kin1 -cd_music -1 -6 0 0.25 5 -cd_music gives the charge contribution of the surface species to the three planes. Plane 0 is -1 + 0.25*5; Plane 1 is -6 + (1-0.25)*5; Plane 2 (or d) is 0. svn 1030: Fixed bug in tranport. Mixing was not printed when using -cvode in kinetics. svn 984: Fixed bug in transport when cell without a surface followed a cell with a diffuse-layer surface. Fixed bug in TOTAL function; code ran of the end of the list of master species; changed logic to recognize the end of the list. svn 874: Fixed bug in check_same_model. Thought surface calculation was the same even though -edl switch was different, which gave irratic results and possible crash. Now checks more carefully to make sure calculation for surfaces is really the same and reinitializes if not. svn 847: Fixed bug with DESCRIPTION function. Did not give correct solution description for reactions. svn 826: Update tally.c to avoid conflicts in C++ version of phast. svn 801: Wrote around underflow in fabs in subroutine reset. svn 794: Errors in minteq.v4.dat database. Several redox reactions had delta H listed as kcal instead of kJ. kcal is correct only for the following species H2, NO2-, and NH4+. 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