Frequently Asked Questions for PHREEQC and PhreeqcI

STABILITY CONSTANTS FOR SEAWATER: We focus our investigations in marine environments and we are concerned about the database file PHREEQC.DAT with respect to seawater instead of fresh- or groundwater. There is the dependence of log K data from ionic strength (or salinity) and it is not exactly clear which circumstances the listed dataset represents, but surely no marine environment.

The program seems not to correct log K data by itself, if the ionic strength increases by adding, i.e., Na and Cl. Do such corrections have to be added manually to the database file by changing the log K for the corresponding reaction? Are there any database files with log K coherent data for seawater available?

The log K's in phreeqc.dat are for infinite dilution, consistent with thermodynamics and the free energies of formation of the species. The log K is assumed to apply at all ionic strengths because the mass-action equations are written in terms of activities. All ionic strength corrections go into the activity coefficients of the aqueous species. Either a Debye-Huckel expression with an ion-specific b-dot ionic-strength term (b*I) or the Davies Debye-Huckel expression (which depends only on the charge of the species) is used to calculate activity coefficients.

If you use conditional constants for seawater, the mass action equations are written with molalities rather than activities and all ionic strength corrections go into the log K. Effectively, PHREEQC automatically uses
```
log K(conditional) = log K(thermodynmaic) * gamma(reactants)/gamma(products) 
```
where gamma(reactants) is a product of activity coefficients calculated from the ionic strength of the solution for the reactant side of the balanced reaction, each activtiy coefficient raised to the power denoted by the stoichiometric coefficient. Gamma(products) is similar for the product side.

The activity coefficient parameters of the major ions are derived predominantly from sodium and chloride mean-activity-coefficient data. They should work best in sodium chloride solutions like seawater, but "best" is probably not as good as some of the conditional constants for seawater. The Davies equation is used for activity coefficients of ion complexes.

If you want to use conditional constants, it would be necessary to redefine the log K's to the conditional constants and force activity coefficients to 1.0 (by defining a very large "a" term and a 0 "b" term in the Debye-Huckel expression for each ion. For example:
```
Na+ + SO4-2 = NaSO4-
        log_k           0.700 #redefine with conditional constant
        -gamma    1e6    0.0
```
I do not have a file that is set up with conditional constants for seawater. If you know the constants from other sources, it should not be hard to set up the database file.
TOTAL CO2 and TOTAL CARBON: I'm wondering why, in the "Description of Solution" section of my output, Total Carbon and Total CO2 are both given. For all of my simulations, their concentrations are the same. How are these parameters calculated, and are their any situations where they would be different?

Total CO2 can also be called "dissolved inorganic carbon" and it is the sum of all the CO2, HCO3-, and CO3-2 species plus their ion pairs and other complexes. Total carbon includes DIC plus methane, so Total CO2 and Total Carbon differ only if you either enter methane in a solution composition or generate methane in a reaction calculation,
HIGH IONIC STRENGTH and REDOX REACTIONS: I am currently looking at the geochemistry of two groups of saline perennial spring discharge in the Canadian High Arctic. The springs discharge at a constant temperature around +6.5 C and precipitate calcite along with smaller amounts of undetermined iron- and sulfur-based minerals. Sulfate-reducing bacteria are believed to exist in the springs, associated with a notable presence of H2S at the outlets. The pH of the waters rise with increasing distance from the spring outlets, leading me to believe that CO2 degassing is the dominant mechanism for calcite precipitation.

One of the major goals of this research is to examine the re-equilibration processes leading to mineral precipitation at these sites. Because of the moderate ionic strength of the waters (1.45mol/kg), we are restricted from modeling these processes with contemporary geochemical models such as WATEQ and PHREEQ. As a result, I have been working with a program developed at CRREL called FREZCHEM2, which models solutions of high ionic strengths under both (1) freezing and (2) evaporation conditions. The model is currently limited by the fact that carbonate species are not included in the program, thereby limiting calculation for changes in pCO2 and precipitation of CaCO3.

Therefore, I find myself at a loss for progress. I am interested in modeling how the system will respond to changes in temperature (simulating summer and winter seasons), but I do not understand (1) how models work with solutions with high ionic strengths and (2) how models would furthermore extend this situation to low temperatures (below freezing). I would also like to observe mineral precipitation as a result of freezing experiments. There is abundant literature on the theory of evaporation of brines and associated field/laboratory observations, however the same does not hold true for freezing environments, and I am unaware of any system that has led to the precipitation sequence of laminates that we observe in these springs. Do you know of any work in this area, or papers examining any of these topics?

You are caught between the ion association modeling approach and the Pitzer specific interaction modeling approach. I will describe two models that I have helped develop that may be applicable to your problems, PHREEQC (ion association) and PHRQPITZ (Pitzer), but both have their limitations. Other codes you should consider that implement the Pitzer specific interaction approach are the Geochemist's Workbench, EQ3/6, and SOLEMNEQ.

PHREEQC is a very capabable reaction/transport model that easily handles carbonate/sulfide precipitation, redox reactions, and version 2 has a very general kinetic reaction capability. PHREEQC would be able to model evaporation or freezing by removal of water. However, PHREEQC uses an ion-association model to calculate thermodynamic properties and this model starts to break down at higher ionic strengths, somewhere between .1 and 2 molal, depending on solution composition. The activity coefficients of the major cations ions are fit mostly from chloride salts, the model may do reasonably well for majors in chloride dominated waters. However, all of the ion pairs and complexes use the Davies equation for activity coefficients, which is purely speculation. So any calculation that involves large amounts of complexing at high ionic strength (>0.5) is questionable.

PHRQPITZ is a reaction model (no transport) that uses the Pitzer specific interaction approach for modeling thermodynamic properties. The parameters are fit to mixed electrolyte solution data and reproduce observations to relatively high ionic strength, frequently 2 molal or higher. PHRQPITZ data base includes carbonate species, major elements, and some trace elements. For your work, the deficiency is that it is not possible to model redox reactions with iron and sulfur. EQ3/6 and the Geochemist's workbench are other geochemical models that implement the Pitzer approach that you should check out, but as far as I know, they can not model redox reactions with the Pitzer approach either. PHRQPITZ does not do as good a job with water mole balance as PHREEQC, but it is probably adequate.

I believe FRZCHEM implements the Pitzer approach. The aqueous model for PHRQPITZ is defined in database files, so I think it may be possible to incorporate the parameters and temperature dependence from FRZCHEM into PHRQPITZ. I would think the PHRQPITZ parameters for carbonate would be adequate at 6.5C, but you may need to revise parameters for the carbonate system if you need to model temperatures below zero. Unless you are lucky and somebody has worked them out, devising Pitzer parameters is a lot of work.

For a start, my strategy would be to make some nonredox calculations at 6.5 with PHREEQC and use PHRQPITZ (or FRZCHEM where applicable) to assess the reliability of the calculations. I think it is possible that PHREEQC will be qualitatively adequate.

I never pursued the work, but I did publish an aqueous model in (D.L.Parkhurst, Ion-Association and mean activity coefficients of various salts, ACS Symposium Series 416, American Chemical Society, Washington DC, pp 30-43, 1990) that was a mixture of Pitzer and ion association approaches. I'm a little leary of what will happen with this model if you get outside the range or solution compositions over which the parameters were fit.
REDOX DISEQUILIBRIUM: I cannot find out how to 'tell' PHREEQC not to react nitrate with ammonium. Although this reaction is thermodynamically possible (favourable) it is not occuring in my aquifer (at least to any great scale). As my contaminant plume is principally composed of nitrate and ammonium (~ 500 mg/l of both) it makes carrying out any 'realistic' equilibrium, exchange or transport simulations impossible as PHREEQC always reacts nitrate with ammonium, which in my case has a large effect on the water chemistry. Please help! Thank you in advance and apologies if I have just been stupid and not been able to find the command in the manual.

You haven't been stupid, it is not obvious how to create redox disequilibria. You must redefine the database to remove the connection between the redox states. In effect, you have to define new elements that don't interact with each other. I've outlined the changes to phreeqc.dat below. The plan is to remove nitrite, dissolved N2, and ammonium from the definition of "N" species. They are replaced with separate elements "Nzero" and "Amm".

Now the danger of this strategy is that the three valence states (5, 0, -3) are essentially inert. There is no longer any way for homogeneous redox transformations to occur in the nitrogen system. You can define a REACTION that transfers N to Nzero or Amm, but that is rather awkward in version 1 because the reaction rate is constant. If you are using the version 2 (beta test version), it is possible to define kinetic reactions that perform the transformation that can take into account concentrations of oxidants and reductants and concentrations of the various nitrogen valence states.

I don't know which way your reactions are going, but a variation on the kinetic approach that slowly oxidizes the ammonium would be as follows. Retain the complete definition of N plus add a definition of "Amm". Initial conditions define nitrate as N(5) and ammonium as Amm. The N species are always at equilibrium, whereas Amm reacts only by a kinetic reaction. The kinetic reaction transfers Amm to the N, with the effect that eventually all the Amm is transferred to the N pool and redox equilibrium is attained.

I've eliminated nitrite in the following example, but if you leave in the definitions of N(5) and N(3) (SOLUTION_MASTER_SPECIES) and the definition of NO2- in SOLUTION species, you would have a partial equilibrium; nitrate and nitrite would react to equilibrium with each other, but would not react with Nzero and Amm.
```
SOLUTION_MASTER_SPECIES
N        NO3-           0.0     N               14.0067
Nzero   Nzero2  0.0     Nzero   14.0067
Amm		AmmH+	0.0	AmmH		18.0
#N(+5)    NO3-           0.0     N
#N(+3)    NO2-           0.0     N
#N(0)     N2             0.0     N
#N(-3)    NH4+           0.0     N

SOLUTION_SPECIES
Nzero2 = Nzero2
        log_k 0.0

AmmH+ = AmmH+
	log_k	0.0
	-gamma	2.5	0.0

AmmH+ = Amm + H+
	log_k	-9.252
	delta_h 12.48	kcal

AmmH+ + SO4-2 = AmmHSO4-
	log_k	1.11

#
#The following species definitions need to be removed
#

#NO3- + 2 H+ + 2 e- = NO2- + H2O
#        log_k           28.570
#        delta_h -43.760 kcal
#        -gamma    3.0000    0.0000

#2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
#        log_k           207.080
#        delta_h -312.130        kcal

#NH4+ = NH3 + H+
#        log_k           -9.252
#        delta_h 12.48   kcal
#        -analytic    0.6322    -0.001225     -2835.76

#NO3- + 10 H+ + 8 e- = NH4+ + 3 H2O
#        log_k           119.077
#        delta_h -187.055        kcal
#        -gamma    2.5000    0.0000

#NH4+ + SO4-2 = NH4SO4-
#        log_k           1.11

PHASES
Nzero2(g)
            Nzero2 = Nzero2
            log_k -3.26
Amm(g)
	Amm = Amm
	log_k	1.77
	delta_h -8.170	kcal

EXCHANGE_SPECIES
	AmmH+ + X- = AmmHX
	log_k	0.6
	-gamma	2.5	0.0


END
	
```
ANALYTICAL SOLUTIONS OF TRANSPORT: I am doing transport calculations and compare PHREEQC results with an analytical solution, but results do not compare well?

If PHREEQC fronts are more disperse: try a grid refinement and use more cells. A finer grid will reduce numerical dispersion.

If midpoint concentrations are displaced: Was the boundary condition {dc / dx = 0 at x = infinite} used for the analytical solution? And, do you compare concentrations at a point x which is the last cell of PHREEQC? If so, add one or more cells to make the PHREEQC flowtube 'infinite' to comply with the analytical solution.
IONIC STRENGTH LIMIT: I've been perusing the user's guide to PHREEQC and note that it works well at low ionic strength. What's the upper ionic strength limit you recommend?

You also say it MAY be reliable at higher ionic strength solutions for sodium chloride-dominated systems. How high do you mean by higher? (upper limit). Does seawater qualify as sodium chloride-dominated, or do you mean virtual absence of ions other than sodium and chloride?

Hedging, I say between .1 and 1 molal (seawater is ~0.7). A non-sodium- chloride solution would be toward the lower end and a sodium-chloride (or at least chloride-) dominated solution would be toward the higher end.

For the major cations, the activity coefficients as a function of ionic strength are derived from chloride salts. In a pure sodium-chloride solution, the activity coefficients of sodium and chloride are fit to perhaps 6 molal. For other strong electrolyte cations that are fit from chloride salts, the calculations are probably not too bad as long as chloride is the dominant anion. As anions other than chloride are introduced, things get progressively worse. In addition, most complexes have only the non-ion-specific Davies equation for the activity coefficient, so large amounts of complexing at high ionic strength is much less reliable. In general I would be very cautious (skeptical) with calculations at ionic strengths over 2 molal. Between .5 and 2 molal, I would consider the calculations as qualitiative and check to the extent possible with the specific interaction approach.

You can use other programs (EQ3/6, PHRQPITZ, Geochemists Workbench) that implement the specific interaction approach of Pitzer. This approach fits activity coefficients, frequently to 2-6 molal, in mixed electrolyte solutions. This is definitely the better approach for high ionic strengths, but the drawback is the data are not available for redox reactions or aluminosilicate reactions.

MODELING DIFFUSION: Diffusion with version 2 does not give results that compare well with the Ogata-Banks solution. We ran the file below, for diffusion through a 1 meter-thick layer with a constant chloride concentration at the upper boundary (C = 1). To match the Phreeqc result with the Ogata solution, we have to add some convection (i.e. a velocity of 5.2e-5 m/day). Have you come across similar difficulties, or is there an aspect we have overlooked ?
TITLE example.--Diffusion only (all time units in days) (20 cells) SOLUTION 0 CaCl2 units mmol/kgw pH 7.0 charge temp 25.0 Ca 0.5 Cl 1.0 SOLUTION 1-20 Initial solution for column units mmol/kgw pH 7.0 charge temp 25.0 TRANSPORT -cells 20 -shifts 20 -time_step 5.e3 -flow_direction diffusion_only -boundary_conditions constant -lengths 0.05 -diffusion_coefficient 2.59e-5 -punch 20 PRINT -reset false SELECTED_OUTPUT -file ex9dif.pun -totals Cl END

You have been defining an extremely large diffusion coefficient. Note that -diffc should be given in m2/s, and that 2.59e-9 is reasonable. With your diffusion coefficient of 2.59e-5 m2/s you have no longer the boundary condition at x = 1 that dc/dx = 0 for the analytical solution. It is better, still, to first calculate how far the diffusion front will come with the analytical solution, and then dimension the PHREEQC grid.

        
TITLE example.--Diffusion only (all time units in days) PHREEQC-2 WORKS IN
SECONDS (it is so fast...)
(20 cells) 
SOLUTION 0  CaCl2
        units   mmol/kgw
        pH      7.0     charge
        temp    25.0
        Ca      0.5
        Cl      1.0
SOLUTION 1-20  Initial solution for column
        units   mmol/kgw
        pH      7.0     charge
        temp    25.0
TRANSPORT
        -cells  20
        -shifts 20
        -time_step 5.e3   # in seconds
        -flow_direction diffusion_only
        -boundary_conditions constant
        -lengths 0.05
        -diffusion_coefficient 2.59e-9    # in m2/s
        -punch 20
PRINT
        -reset  false
SELECTED_OUTPUT
        -file ex9dif.pun
        -totals Cl
END

INSTALLING PHREEQCI: I cannot expand the PHREEQC-I download. When I execute it, it comes up saying that there is an error in zip!
You should have downloaded the file pci1_03.exe which has a size of 3885042 bytes. Check to make sure this is the right file and the right size. If you used ftp you need to use "type binary" or "type image". It may be that your browser made an ascii translation that it shouldn't have or there may have been some errors in transmission.
INSTALLING PHREEQCI: Once I finished inputting data, I return to the main screen and click on "Run Current File". After specifying the output file name (I'm using the default), I click on "OK". At this point, I get an error message saying that the output file cannot be found. I'm running Windows 3.11. I tried your suggestion of going to the directory in which phreeqci is installed (c:/phreeqci by default) and copying examples\ex1.in to the installation directory. Then I ran phreeqc with the following command: phreqc32 ex1.in ex1.out phreeqc.dat When I did this in a DOS environment, both with and without Windows running. I get an error message that says "This program must be run under Win32." What's going on?
I think you are missing the program called Win32s for Windows 3.1. (Windows95 and WindowsNT have this package as part of the operating system.) As far as I know, there should be a directory win32s in the \windows\system directory. Windows 3.1 doesn't have the capability to run 32-bit programs without this package. You can download Win32s from microsoft or from our site. Check out Win32s installation instructions at URL: https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqci
ELEMENT REDOX STATES: Using the SOLUTION keyword, as you suggested, I input my own data. I feel about 90% confident that I'm doing this correctly (that 10% uncertainty mainly has to do with how to enter SO4 data--must it be entered as S(6), or can I type in "SO4" myself?)

"SO4" will not work, element concentrations are always entered solely as the element symbol, optionally with a valence state in parentheses. S(6) is safest because it means specifically sulfate (sulfur of valence 6). You could also enter "S", but that would be interpreted as total sulfur and would be distributed between sulfate and sulfide according to the pe (or redox couple). If the pe is sufficiently low, the sulfur could end up as sulfide instead of sulfate in the distribution of species.
INSTALLING PHREEQCI: I am having difficulties running the setup program for PhreeqcI. It says there is a program conflict even when nothing is open to conflict with it. Have you had any complaints about installation of PhreeqcI on a NT 4.0 operating system?
First do you have administrator rights?

Try running the setup.exe program through Start/Run on the taskbar. Turn on the "Run in Separate Memory Space" click box.

If that doesn't work, check the task manager. Right click on a blank area of the task bar (or ctrl-alt-del) and select task manager. See if you have any old versions of phreeqci running or something that looks suspicious under the Processes tab. It will be listed as an item under NTVDM.EXE. I think under some conditions there can be a conflict with a previously running program called wowexec.exe.
PAPERS USING PHREEQC: I am teaching a new course winter quarter on Environmental Chemistry. It would be useful to review how the program has been used successfully in the past. Do you keep some record of publications that have used PHREEQC?
I need to start being more careful in keeping records. I can give you a few references. If anyone emails additional references, it would be appreciated. Click for a list of references. I'll try to keep the list as current as possible. I've included references to PHREEQM, which is a precursor of PHREEQC that had dispersive transport capability. PHREEQC version 2 contains all of the capabilities of PHREEQM.
INVERSE_MODELING: I keep getting an error message "Br is included in solution 1, but is not included as a mass-balance constraint". I get this message for most of the constituents in all of my solutions. It appears as though the inverse-calculation is being made based on pH,Alkalinity, C(4), C(-4), Ca and O concentrations only, is this proper? In order to get the model to run I had to set the uncertainty as high as 20 percent on some of the runs. What does the error message mean and is it ok to set the uncertainty so high?
You have to add some more input to do pure mixing properly. Normally, there are mineral sources and sinks for most of the major elements and PHREEQC automatically includes mole-balance equations for all of the elements included in the phases. If there are elements for which there should be mole-balance equations, but these elements are not in the phases, then you have to request a mole-balance equation using the -balances option. For your mixing case, where you have few or no minerals, then you have to request the mole-balances manually. You should include all of the major elements Ca, Mg, Na, Cl, S, C and probably K. You may want to include Br. (If you have isotopic data for D and O-18, the beta test version 2.0 will handle these as well. https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc.)

Ideally, mixing will account for all of the elements simultaneously. However, you may not be able to account for all of the elements by mixing alone. If this is the case, you will have to rethink your conceptual mixing model and consider whether the waters are related, you have all the endmembers, and if there are plausible reactants for the system.
INVERSE_MODELING: I am trying to interpret the inverse model solutions (there are up to 11 separate solutions for each one of the simulations). Do I pick out the simulation which has the lowest sum of residuals and lowest fractional error? or, If these results seem out of whack, Do I use my own judgement too and pick a model with reasonable results and as low an error as possible?
If you run the program without the -minimal option, the first model that is printed has the smallest sum of residuals. However, other models may have an equally small sum of residuals. The sum of residuals is essentially how much the analytical data have been fudged. Think of a residual of 1.0 meaning that one analytical datum has been changed by its maximum uncertainty; 10.0 meaning that 10 analytical data have been changed by their maximum uncertainty.

I usually include the -minimal option, which decreases the number of models. This mode of running generates larger residuals because it tries to eliminate phases by changing more of the analytical data. The first model will not necessarily have the smallest sum of residuals. However, it seems to give the simplest and most important reactions. Kind of an Occam's razor approach. At any time that you think the data have been unduly adjusted, you can decrease the uncertainty for any or all of the data. You generally can't use an uncertainty of 0 for all of the data because the program has to adjust things a little to obtain charge balance.
THERMODYNAMIC DATA: I was wondering if you have any recommendations for geochemical modelling for chromium. EPA apparently considers the thermodynamic data for chromium in MINTEQA2 to be inadequate. Do you know of a good source of reliable thermodynamic data for chromium?
You could try the EQ3/6 database, it used to be on the web, but I have been unable to find it recently. Ball and Nordstrom have been working on an evaluation of chromium thermodynamic data. Chromium is also in the MINTEQ database in the distribution of PHREEQC. I had some questions about some of the species and particularly about some of the chromium species. However, I never got a response from EPA.
SOLUTION: If I have water analysis with: HCO3- (x mg/l) and free gas content CO2 (y mg/l) How do I write the input file for PHREEQC?
Phreeqc has a few options for carbon, I'm not sure if any fit your data exactly:

(1) Enter pH and alkalinity. Total moles of carbon(4) is calculated by the program.

(2) Enter pH and total carbon(4). Alkalinity is calculated by the program.

(3) Enter alkalinity and total carbon. pH is calculated by the program.

(4) Enter pH, estimate of total carbon, and a partial pressure of CO2. Program will adjust total carbon until desired PCO2 is attained.

(5) Enter total alkalinty or total carbon, estimate of pH, and a partial pressure of CO2. If possible, the program will adjust the pH to attain the desired PCO2.

Now if by "free gas content" you mean the total of all the carbonate species (H2CO3, HCO3-, CO3-2, plus pairs) then you can use method 3 to obtain a pH. Otherwise, you probably want to enter pH and alkalinity and compare calculated CO2(aq) to your measured values. Let me know if you still have questions.
EXAMPLE 1 of manual--Uranium in Seawater: The input concentrations for components were provided in ppm. The output concentrations are reported in molality and moles, which are equal as 1 kg of seawater is assumed. How do the initial ppm concentrations relate to the reported total moles for each constituent? For example, Cl was input as 19353.0 ppm. The output indicates a total moles of Cl equal to 0.5657. The total moles equates to 20055.7 mg of Cl(-1). Using the input seawater density of 1.023 kg/L does not seem to solve my problem.
Ppm is mg/kg solution. Molality is mol/kg water. The conversion is as follows:

Convert mg of a solute to moles of a solute. Subtract the mass of the dissolved constituents from the kg of solution to get the mass of water. Divide moles of solute by mass of water to get molality. The initial solution is then assumed to have 1 kg of water and moles of solutes are scaled accordingly.

Density does not affect the calculation unless input units are per liter.
GUI, Graphical User Interface: I am most interested in: PHREEQC. I have been able to run the examples, but I cannot find the program that allows you to create the input file needed to run PHREEQC. Is there a interactive module of the program (such as phrqinpt for PHREEQE) which helps you input modeling information?
You're in luck. We have released PhreeqcI, an interactive version for PCs and I've put together a web site for downloading the codes. Try https://wwwbrr.cr.usgs.gov/projects/GWC_coupled. This program will be available soon at http://water.usgs.gov/software.
SOLID_SOLUTIONS: Are there any possibilities to model solid solutions directly in PHREEQC? Or is it still necessary to use some other code (eg. SOLISOL) additionally to PHREEQC?
Solid solutions were not included in the 1995 release. However, a beta test version of PHREEQC version 2 is available (https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc.). You can model binary, nonideal solid solutions and multicomponent ideal solid solutions. It is capable of dissolving and precipitating solid solutions in equilibrium with an aqueous phase. Unlike SOLISOL, the solid solution need not always be present initially. The main deficiency of the model is that the entire solid solution reequilibrates with each reaction, whereas it is more reasonable that only the outer layer should reequilibrate.
INVERSE_MODELING: Given that the uncertainty of all elements within a system are equal, how does the program decide which element(m) to adjust (solving for the unknown "del-m,q") in order to achieve a charge balance?
The program allows for the possibility that any of the ions could change within their uncertainties. Which ones it actually changes is in part to make the mole balances work out for a given set of minerals, in part determined by the magnitude of concentrations (given a choice between changing a small concentration or a large concentration by the same absolute amount, it will choose the large concentration), and in part arbitrary in that it may not be a unique choice of deltas for a given set of phases.
INVERSE_MODELING: How and why does the program opt to adjust only one element as opposed to adjusting a group of elements within their specified uncertainty limits? Given a net positive charge imbalance, per say, how and why does the program decide to either add cations or subtract anions to achieve a balance?
The objective function is to minimize the sum (delta/Uncertainty), so if using group would decrease the sum, it should do it. In practice, it probably adjusts the largest concentration, because that gives the smallest delta/uncertainty.

Pretty much the same answer as the first question. The choice of deltas is required to produce charge balance, the choice of deltas is required to have mole balance for all the elements with mole transfers of a set of minerals, and, in addition, the sum of the deltas/uncertainties is required to be a minimum for that set of minerals. The program solves for combinations of deltas and mineral mole transfers that satisfy all these criteria.

I published a paper that describes the approach, but it may seem a repeat of the PHREEQC manual: Parkhurst, D.L., 1997, Geochemical mole-balance modeling with uncertain data: Water Resources Research, v. 33, no. 8, pp. 1957-1970.
REDOX POTENTIAL: In the general frame of nuclear wastes storage, we have performed leaching experiments of uranium dioxide. I would to calculate solubility limit of uranium secondary phases under my experimental conditions. I have determined experimentally pH, Eh, temperature and water composition (Na, Cl, HCO3-, SO42-,K, Mg, Ca and U). After calculations, I obtain solubility limits for phases but the pe varies (-3.08 to -5.9) during calculations. Can you fix the pe for reaction calculations?
I generally discourage people from trying to fix pe during reaction calculations (it's fine for speciation calculation). The reactions should determine the pe just as the reactions determine the pH. Bad things can happen, for example it is possible to slip out of the stability field for water if pe is fixed, say at low pH initially, and then the pH rises.

If the modeled pe is not what you expect, I think there are two possible reasons (1) redox disequilibrium in your system, or (2) unexpected reactions.

Disequilibrium can be modeled with some effort. It is necessary to rewrite the database to split out each redox state as an "element". If all redox elements are split, then all redox reactions would have to be performed with REACTIONs that transfer elements from one redox state to another. Partial disequilibrium may be more reasonable.

You can fix pe the same way as pH. However, you must choose a reactant that enters or leaves the solution to cause the pe to be fixed. It is best if this is a plausible reactant. The example below assumes atmospheric oxygen entered your system? Oxygen could also leave the system, which would be implausible.
```
 PHASES
 Fix_pe
 e- = e-
 log_k 0.0
 EQUILIBRIUM_PHASES
 Fix_pe  4      O2(g)   # fixes pe to be -4 by adding or removing O2
```
INVERSE_MODELING: What is the relative importance of the fractional error value and the sum of residuals? I am not quite sure what they represent.
The fractional error should be less than the global uncertainty that you set (-uncertainty), unless you have overridden the global uncertainty with specific uncertainties with the -balances option, in which case the maximum fractional error should be less than the larger of the global and the specific uncertainties. I look at the sum of residuals as a measure of the number of analytical data that were adjusted by their maximum amount. The larger the sum of residuals the more data that were fudged by their maximum amount. Given two models, the one with the smaller sum of residuals is more consistent with the original data.
INVERSE_MODELING: I am using only two solutions, one is pure water, and the second the geochemical analysis of one flushed pore volume. The manual states that if only two solutions are used, somewhat equivalent amounts will be compared - ratio 1:1, but I get many models that uses zero of the first solution and 1 of the second. Could this be because it is pure water, and doesn't contain any of the aqueous phase contaminants?
You can generate water through chemical reactions. For example if you dissolve 25 moles of gypsum, you generate 50 moles of water (~1 kg). Usually models with very large mole transfers (> 1 mole) are not realistic. As a hint, make sure you charge balance the pure water in these problems. The program has trouble with small numbers, like the alkalinity of pure water, which is around 1e-10 if the water is not charge balanced. Otherwise, just ignore models that use 0 of the initial solution.
NONUNIQUE REACTION TO A PHASE BOUNDARY: Here are below two input files which should give the same results because the solution 1 is in equilibrium with a log pCO2(g) of -2.9 (as a result of simulation). However if we don't equilibrate it with log pCO2(g) of -2.9, it gives a wrong result, with a pH of 10 and a log pCO2 of -10. What's going on?
It's not that one answer is wrong and the other right, there are two answers to the problem you formulated. The one you expected with pH 7, pCO2 near -2, is what you get if you add a little CO2, driving the pH down and making the concentration of CO3-2 smaller until you reach calcite equilibrium. The other answer gets to calcite equilibrium by taking CO2 out; if nearly all the carbon is removed, the concentration of CO3-2 can be made small and thereby attain equilibrium with calcite.

Unfortunately, there is no way to know in general, which answer PHREEQC will find. It appears to be a random result of the numerical procedure. The slight difference between equilibrating with CO2 at -2.9 seems to be enough to switch from one answer to the other. In your case, it is possible to force the pH-10.5 answer by having no CO2 reactant initially (replace 10.0 with 0.0 in equilibrium_phases), but I don't know how to force the pH-7 answer.
THERMODYNAMIC DATA: I am interested in thermodynamic data base accompanied PHREEQC code, however, unfortunately only logK value and enthalpy of reaction are included in data base. I would like to know, Which thermodynamic data are used to calculate logK and enthalpy of reaction for each species and minerals?
The thermodynamic data in PHREEQE is included in a database file for PHREEQC called phreeqc.dat. This data is largely a subset of another database file for PHREEQC called wateq4f.dat. As the name implies, wateq4f.dat is derived from the thermodynamic data included in the speciation model WATEQ4F. Sources of thermodynamic data are included in the documentation for that program:

Ball, J.W. and Nordstrom, D.K., 1991, User's manual for WATEQ4F, with revised thermodynamic data base and test cases for calculating speciation of major, trace, and redox elements in natural waters: U.S. Geological Survey Open-File Report 91-183, 188 p.

Check the web site https://water.usgs.gov/software/lists/geochemical/ for more information on WATEQ4F and how to order a manual.

PHREEQC also has a database file derived from MINTEQA2, called minteq.dat, but I don't think that MINTEQA2 documents the sources for these thermodynamic data.
DISPERSIVE TRANSPORT: Presently I am using Phreeqc v1.5 to model our experimental breakthrough curves (BTC's) from column experiments. I have fitted breakthrough curves of bromide from column experiments using CXTFIT with very good results. However, because of its greater "geochemical versatility", I would like to use Phreeqc instead of CXTFIT to model the observed BTC's for Br and other ions (e.g. Na, SO4). I have previously used PhreeqM for this purpose but it doesn't do inverse modelling and to be consistent I want to use only one program. Is is there a Phreeqc code that does transport with dispersion?
A beta test version of PHREEQC version 2 is available that allows 1D advective/dispersive (or diffusive) transport (https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc.).
KINETICS: We have heard about a further version of the code which include diffusion and kinetics of dissolution for some minerals. If this is the case, I wonder whether it would be possible for us to get a copy of the new version of the code for PC?
A beta test version of PHREEQC version 2 is available that allows general rate expressions for reaction kinetics. (https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc.).
TESTING AND QUALITY ASSURANCE: We are currently using PHREEQC and were using PHREEQE in the past in our group. We are interested in how you tested these two geochemical codes. Did you compare them with other codes? did you use a particular kind of test? How are you sure that the data you get are 'correct'? Could you please send us the test procedure when convenient?
PHREEQE was tested largely by hand calculations. All of the equations are algebraic, so the results can be tested directly. Activity coefficients, ionic strength, mole balance, saturation indices and other quantities can all be tested from the printout.

PHREEQC was tested against PHREEQE results for the PHREEQE test problems. It was tested against WATEQ4F for the wateq4f.dat data base. PHREEQC results were very close to the other programs, with slight differences in the WATEQ4F results attributed to slightly different activity coefficient formulations.

The surface complexation modeling was tested against a couple of test problems in MINTEQA2.

We have currently been working on dispersion and stagnant zones. These parts of the code have been tested against analytical solutions. Kinetics have been tested against test problems developed by Valocchi and more recently Valocchi and Tebes. Solid solutions have been tested by hand calculations and against MBSS, a code for solid solution calculations.

The code has been used extensively, and any bugs have been tracked down and corrected. There have been other efforts, particularly by a company called Intera, to compare PHREEQE results to EQ3/6. I know of two current efforts along these lines as well as with MINTEQA2. One final check is performed within the code to test the mole balance and mineral equilibrium algebraic equations before results are printed. In all, there have been few bugs that resulted in incorrect results. Most program failures are due to nonconvergence in the numerical method.
CORROSION CONTROL: I'm working with a municipality that uses the lime-soda ash softening process to treat their drinking water (groundwater source). This treatment process raises the pH to greater than 10 and reduces alkalinity, making the water extremely aggressive. In reviewing the literature, the American Waterworks Assn. recommends using PHREEQ for corrosion control. Are you aware of any other instances where the program is actually being used for this purpose? If so, do you know of any input files that have been created for this purpose or any contact names/phone numbers of those using the program for this function?
I haven't had much contact with people that are using it for this purpose other than talking to the authors of the AWA report, I think it was nearly 10 years ago. However, I expect PHREEQC is being used for this purpose. Sorry, I'm not much help here.

To get you started, a simple input file would look like this. Portlandite is in the wateq4f.dat database, which you can set in the Options menu of PHreeqcI.
```
SOLUTION  1     Pre-treatment water; replace with your analysis
        units   mg/L
        pH      7.0
        Na      15.
        Ca      60.
        Mg      40.
        SO4     10.
        Alkalinity      250. as CaCO3 
        Cl      10.
END
USE solution 1
REACTION 1
# this adds .01 mole of CaO and .005 moles of NaOH to solution 1
# replace with stoichiometry and amount of lime and ash added per liter
        CaO     1.0
        NaOH    .5
        .01 mole        
EQUILIBRIUM_PHASES 1
# This allows portlandite and calcite to precipitate if they
# become supersaturated
        Portlandite     0.0     0.0
        Calcite         0.0     0.0
SAVE solution 2
END     
# solution 2 is treated water which could be used for additional
# reaction modeling such as pH adjustment?
```
MAC VERSION: I'm interested in the geochemical model named PHREEQC which can run on Macintosh computers (PowerMac 7.1 or 8). If versions for such computers are available I would gratefully appreciate in receiving such a version. If not, is it possible to compile the source code for Macintosh, and, what compiler should be used in this case?
Brian Marshall of U.S. Geological Survey, Lakewood, CO (bdmarsha@usgs.gov) has created a Mac version that is available at https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc.
INVERSE_MODELING: What does it means if the output gives you no phase transfers as the minimum model to account for changes in water quality between two solutions?
It means the two solutions are essentially the same composition within the given uncertainty limits and no reactions are needed to account for the differences. However, if the uncertainty limit were 1.0 for each constituent of each solution, that is 100 percent uncertainty, then any two solutions would be essentially the same. If uncertainties are small, say less than .05, then it means concentrations between the two solutions differ by less than 5 percent. The key is whether you have used uncertainty limits that are reasonable, if you have, then the solutions are not significantly different.
INVERSE_MODELING: Is there any way that I might form an output file which lists only the sum of residuals and maximum fractional error for each test run?
There is no way with version 1 of PHREEQC to write inverse-modeling information to another file. The beta test version 2.0 has the capability to write the mole transfers to a selected-output file. I'll have to add the maximum fractional error and sum of residuals to this output. (https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc).
INVERSE MODELING: If I specify a list of minerals without using the minimal identifier, does this mean that the model will use every mineral given and try to equilibrate these?
The inverse modeling does not do ANY equilibration, it is purely an exercise in accounting for moles of elements. If you do not put in -minimal, the program attempts to find every combination of minerals that can satisfy the mole-balance constraints, subject to changes in analytical data that are within the specified uncertainty limits. The -minimal option limits the models to those sets of minerals for which no model can be found with a subset of that set of minerals.

In general -minimal gives fewer models at the cost of greater "sums of residuals". The superset models that are eliminated with -minimal will have smaller sums of residuals. However, all models satisfy the uncertainty limits, so implicitly even models with greater sums of residuals are consistent with your description of uncertainty.
INVERSE MODELING: Do you know what may be causing error messages of Exiting due to signal 'SIGFPE' or 'SIGSEGV'?
The problem with your run is with the number of minerals in your list. Though never stated in the documentation, the maximum number of mineral in inverse modeling is 30 to 32.

On the other hand, even with 30 minerals, I expect the execution time will be slow if not infinite and the number of models you get will be very large. You need to simplify your problem.

Consider eliminating some of the trace elements, Se, Cd, Pb, As, U? and related minerals. I think you want to reduce the list to at most 20 minerals, and probably fewer, that you think are reasonable. Some reductions are easy, because polymorphs or phases that differ only in the amount of water in the formula are normally indistinguishable in mole-balance modeling. Thus you should include only gypsum or anhydrite; hematite, goethite, or lepidocrocite; one CdSO4; and probably others. In addition, eliminate any of the minerals forced to dissolve that are absent or present in trace quantities.

Many of the minerals are probably not likely to form, for example, copper sulfate and iron sulfate minerals are very soluble. Anyway, I think if you pare your problem down to a tractable definition, you will have less than 30 minerals and the program will function fine.
TRANSPORT: I was hoping you could explain the way in which PHREEQC treats the TRANSPORT function when used with the MIX function. I would like to simulate the influx of surface water into groundwater and show how the chemistry changes. Does the input that I've written work?
The input should run, but I don't think it does what you were trying to do. For transport, the program calculation first does advection, which is simply moving the water from cell 0 to cell 1, cell 1 to cell 2, ... cell n-1 to cell n. Following advection, mixing and equilibration is performed. You have defined no other reactants so only mixing and aqueous equilibration is performed.

Now if you consider your first cell for an advective step: solution 0 is moved into cell 1, then by the MIX definition .98 of solution 0 is mixed with .02 of the new solution in cell 1 (which equals solution 0) and the result is placed in cell 1. Thus, from advective step 1 onwards, cell 1 is always equal to solution 0. At step two, the same logic will apply to cell 2 and so on. Even if you have no mixing, the surface water will propagate through the system in 50 time steps. The mixing fractions in addition to the transport don't make sense to me. I think mixing solution 1 and solution 2, 2 and 3, and so on would make more sense in terms of simulating dispersion.

Note the beta version of 2.0 has a diffusive/dispersive plus advective transport with boundary conditions. The model also includes a temperature retardation factor that allows you to account for the heat capacity of the matrix (https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc).
SPECIES DEFINITIONS: We are using Phreeqc1.6 to model bytownite feldspar dissolution in aqueous gluconate solutions and we are having some difficulty in figuring out how to put the stoichiometry of the Al:gluconate complexes into the SOLUTION_SPECIES database. The complexes are given in Smith and Martell (Vol 6) as:

Al+3 + L- ---> AlL+2 log K = 1.98 AlL+2 ----> AlH-1L+ + H+ log K = -2.87 AlH-1L- ----> AlH-3L- + H+ log K = -9.29

All constants were determined at 25oC and 0.1M KNO3 by Motekaitus and Martell (1984, J. Inorg Chem 23:18-23). The problem arises from the negative subscripts in the latter two species. According to 84MM, these result from a tridentate interaction between gluconate and Al3+ involving the hydroxyl groups of the gluconate rather than the carboxylic group of the ligand. Normally gluconate

H H H OH H | | | | | H-C--C--C--C--C--C=O | | | | | | HO OH OH H OH OH

is treated as a mono-carboxylic acid and therefore there are no defined deprotonated species in solution equivalent to H-1L-2 or H-3L-4.

How can I define these species in PHREEQC for the purpose of modeling their impact on chemical affinity and on the stability of secondary dissolution products? Will PHREEQC be able to handle the negative subscript for H in the complex formulation without otherwise impacting the calculations?

I think you should define the ligand master species as containing 3 hydrogens or actually the ligand "L" as the deprotonated form. PHREEQC won't let you have negative coefficients, but I don't see any reason you can't define the ligand as whatever state of deprotonation you choose. I've roughed in the input below. Not sure about the last one, your equation has only one hydrogen, but a net change of 2 hydrogen in the ligand. So check the details, but I think it should work.
```
SOLUTION_MASTER_SPECIES
#Alkalinity should be defined relative to stable 
#complex at pH ~4.5

L       LH3-    -3.0        C6H12O7     ~200.?

SOLUTION_SPECIES

        LH3- = LH3-
        log_k   0.0

        Al+3 + LH3- = AlLH3+2
        log_k   1.98

        AlLH3+2 = AlLH2+ + H+
        log_k   -2.87

#Not sure about stoichiometry for this one.

        AlLH2+ = AlL- + 2H+
        log_K = -9.29
```
INVERSE MODELING: Just wanted to send you a file that results in some curious model results. All the models are somewhat nonsensical from my point of view (although computationally legit) since they all dilute solution 1 to make solution 2. What's computationally interesting is that about five models come out with 0 amount of solution 1 to make solution 2.

Reading your manual, I see in the notes on the inverse modeling section that evaporation/dilution can be modelled by adding the phase H2O to the list of phases considered. However, in my experience PHREEQC also occassionally does evaporation and/or dilution even when the H2O phase is not specified (In addition to the above dilution example, one of my files also resulted in spontaneous evaporation).

It is not really evaporation or dilution, but looks kind of like it. The water-balance equation allows minerals to generate or consume water. This is possible even if there are no waters of hydration; the program simply looks and the coefficient for water in the dissolution equation for the mineral and adds or removes this much water per mole of mineral that reacts. You generally need 10's of moles of mole transfer of minerals to generate all the water in solution 2 (55 moles), so not only is the mixing fraction of solution 1 near zero, but the mole transfers of the minerals are unreasonably large.

I just ignore these models. I did add a switch to ignore the water generated by the minerals in version 2 (-mineral_water false). I am a little concerned that this switch might mask some legitimate water effects in some models, so I usually do not use it.
REDOX DISEQUILIBRIUM: When I discovered PHREEQC some time ago, I thought it may be a valuable tool to calculate the speciation in aqueous electrolytes. In order to check the capabilities of the program and to learn its usage, I posed the following simple question: What will be the pH of pure water, if exposed to pure, dry air under ambient conditions? Besides program handling, this question also seemed to have some implications to natural waters. Why is the model model predicting pH ~ 1 due to N2 oxidation?
It is the correct equilibrium calculation, but the confusion is that the natural world does not reach this equilibrium.

As I see it, you must define your reactions as equilibrium or kinetic. You have a choice, does nitrogen react with oxygen to equilibrium or not. You have so far modeled it as equilibrium and do not like the results. Therefore, you are arguing that there is a kinetic reaction. I have modeled the end member kinetic case below in which nitrogen and oxygen simply do not react; nitrogen inert relative to oxygen. The results are as you expected.

To consider the case of kinetic reaction of nitrogen with oxygen with non-zero rate, it would be possible to slowly transfer mass from Nzero (inert relative to oxygen) to N (equilibrium reactions relative to oxygen) through REACTION calculations and the pH would decrease from 5.8 as nitrate was generated by reaction of nitrogen with oxygen.
The following example produces the results you expected. As I said before, your concern in your previous emails is that nitrogen gas does not react to equilibrium with oxygen at low temperature and pressure. Thus, we must remove the default equilibrium among the nitrogen species; SOLUTION_MASTER_SPECIES and SOLUTION_SPECIES are used to do that. In the following example, nitrogen gas is assumed to be inert relative to oxygen. It does obtain equilibrium between the gas and aqueous phase. The problem is defined as follows: pure water reacts with a gas volume that is initially 1.0 L, 25C, and pressure 1.0. The gas is a bubble so the pressure remains 1.0 atm, but the volume will vary as the gas reacts with the water. Initially the partial pressures of dry air, CO2, N2, O2 and H2O--10^-3.5, 0.78, .2, and 0 atm.
```
SOLUTION_MASTER_SPECIES
Nzero   Nzero2  0.0     Nzero   14.0067
SOLUTION_SPECIES
        Nzero2 = Nzero2
        log_k 0.0
PHASES
Nzero2(g)
        Nzero2 = Nzero2
        log_k -3.26
H2O(g)
        H2O = H2O
        log_k  1.5
END
SOLUTION 1
        temp    25.
        pH      7.0
END
GAS_PHASE 1
        -pressure  1.0
        -volume    1.0
        -temp      25.
        CO2(g)     0.00032
        Nzero2(g)  0.78
        O2(g)      0.20
        H2O(g)     0.0
END
USE solution 1
USE gas_phase 1
END
        
```
The results follow. The pH after equilibrium is 5.8 under the specified conditions. The log partial pressure of CO2 has decreased to -3.81 due to dissolution of about half the gas into water. The partial pressure of water in the gas has increased from zero to .03 atm, which is saturation. A small amount of N2 and O2 gases have dissolved. No reaction occurs between the N2 and the O2.
THERMODYNAMIC DATABASE: In your database file - PREEQC.DAT - there is a disagreement in values of equilibrium constants(log K). e.g. in solution_species section, the value of log K in equation 'Pb+2 + CO3-2 = PbCO3' is 7.240 and in phases section, the value of log K in cerrusite is -13.13. the equation is same as former except change in right and left. Why are they different 7.24 and 13.13? and which value is used in the simulation?
The two equations refer to two different things. The equation in SOLUTION_SPECIES refers to an aqueous complex, that is a dissolved ion pair of lead and carbonate. The equation in PHASES refers to a solid mineral (cerrusite). The equilibrium constants are based on the free energy of formation of the aqueous species and the mineral respectively, so they should not be numerically equal. In the simulation, the log k for the aqueous species is used whenever lead and carbon are present in the aqueous solution. The log k for cerrusite is used to calculate the saturation index for cerrusite and whenever cerrusite is included as a mineral in EQUILIBRIUM_PHASES.
CONVERGENCE AND ALKALINITY: I am working on some uranium pore fluids and I have one set of results that will not speciate. The arsenic concentration is to high and the model won't converge. The arsenic concentration was measured at two labs one in the US and one in Canada. The Canadian value was 71,000 ug/L and the US value was 55,523 ug/L total As or 611 ug/L As3+ and 54912 ug/L As5+. The problem is that the file wont converge because there are to many cations I guess. (charge balance is positive 10-3). The file will converge if I drop the As concentration down to 25,000 ug/L. What are my other choices can I adjust the KNOBS to get the model to converge or should I be evaluating the data set to see whether it should be modeled at all?
I ran your input file and the problem appears to be in the specified alkalinity. I removed the alkalinity input and calculated the total non-carbonate alkalinity (> 1. meq), which is greater than the specified total alkalinity (52 mg/L, ~ 1 meq). The program tries to calculate total carbonate species such that carbonate_alk + non_carbonate_alk = total_alk. In this case carbonate_alk has to be negative, which is not possible and the program fails. I think what you found is that the arsenic species account for a lot of the noncarbonate alkalinity and if you decrease arsenic concentration you decrease the noncarbonate alkalinity to the point where carbonate alkalinity can be nonnegative.

You've got a few choices (1) ignore carbonate species by removing alkalinity, (2) specify a total concentration of C(4) instead of total alkalinity, (3) specify total carbon in some indirect way, like equilibrium with atmosphere (but at this pH, that would calculate a large concentration of carbon) or equilibrium with calcite if that makes sense.
EXCHANGE COMPOSITION: One of our students is trying to simulate a batch test, where a known volume of solid (with a known number of adsorption sites) is put in contact with a certain volume of solution of known composition. He needs to calculate concentrations of species adsorbed and in solution. But PHREEQC does not ask volumes as input. If we specify just SOLUTION and EXCHANGE, PHREEQC considers an infinite supply of solution. How should we treat this problem with PHREEQC? (On the other hand this problem is easily solved using MINTEQ).
A solution defined with a solution keyword has 1 kilogram of solution. It sounds to me like you are specifying EXCHANGE (or SURFACE for sorption sites?) in equilibrium with a solution with the -equilibrium option of EXCHANGE. This does assume an infinite volume of the solution and provides the exchange composition in a column after a large volume of water has been eluted.

It also sounds like this is not the model that you want, rather you want to take the an amount of solution 1 and equilibrate it with a fixed amount of the exchanger. In that case, you need only remove the -equilibrium identifier from EXCHANGE. In that case however, what is the initial exchange composition when it is added to the water? You must specify the initial exchange composition explicitly. If you have pre-equilibrated the exchanger with an ion, say calcium, you could specify the initial exchanger as CaX2.

The following is an outline of the type of calculation you could make (and I think it is simpler than MINTEQ). I have put some arbitrary concentrations just for discussion. An exchanger with .1 mmol of sites initially filled with sodium ions is brought together with approximately 200 mL (actually a solution with .2 kg of water) of 1 millimolal CaCl2 solution and the exchanger and solution come to equilibrium. The final solution is sodium-calcium-chloride and the exchanger is nearly all calcium filled.
```
SOLUTION 1
        units mmol/kgw
        pH      7.0
        Ca      1.
        Cl      2.
END
MIX
        1       .2
EXCHANGE 1
        NaX     .0001
END
        
```
COMPILATION: I've been trying to compile the source code of phreeqc using both an old copy of Turbo C and a newish version of Borland C++. Compilation is fine, but no matter which options I try the program fails at run time. Some memory options simply cause a crash, presumably because the program writes/tries to allocate protected memory. Some options are simply too large to load, e.g. program size over 600K with huge memory model and all debugging turned on.

I can get a reasonably-sized executable (around 300 K or less up to about 500K) if I turn off all debugging and use medium or large memory models. When I do this there is a run time error:

Initializing .. ERROR: Null pointer returned for malloc or realloc

Can you advise? What options should be used with Turbo C compiler and linker to ensure a successful result? Should I be using command line rather than the integrated environment? (Just for compilation, I always drop out to DOS before trying to run the program.)

The error message is printed by PHREEQC and means that the program could not get more memory to store data; it ran out of RAM (and swap space if it is available).

Here's what I use with Borland, I would go into the TargetExpert by right clicking on the executable file name in the list of project files. The target type should be "application[.exe]", the platform "Win32", and the Target Model "Console". My checked standard libraries include "Class Library", "Runtime" (disabled), and "BWCC". "Static" makes a portable code that will run on other machines.

The 32-bit compiler does not have any choice of memory models under the options/project menu.
THERMODYNAMIC DATA: I am trying to add some new thermodynamic data I have a table of Dissociation constants for 1:1 metal (Mg+2, Ca+2, Mn+2, Fe+2,Fe +3, Ni+2, Al+2) arsenate complexes (H2AsO4-, HAsO4-2 and AsO4-3) (From Don Langmuir) Do I have to go to the data files and add the species there or can I use the solution species command. If I have to modify the thermodynamic database file where do I add these items at the end of the file and what about the consecutive numbers are they important? this was added to the file as a test

SOLUTION_SPECIES 1 # Do not change this number it is used by the program! MgH2AsO4+ = Mg+2 + H2 AsO4- -log_k 1.52

I get this error

MgH2AsO4+ = Mg+2 + H2 AsO4- log_k 1.52 ERROR: Elements in species have not been tabulated, MgH2AsO4+. ERROR: Reaction for species has not been defined, MgH2AsO4+. ERROR: Calculations terminating due to input errors.

The only thing wrong is that it needs to be an association reaction. In the balanced equation, the species being defined must be the first species to the right of the equal sign. (The program interprets your equation as defining Mg+2.) Sorry it's pretty subtle, but that's the convention I used.
```
Mg+2 + H2 AsO4-  = MgH2AsO4+ 
       -log_k   -1.52
        
```
REDOX DISEQUILIBRIUM: We are doing some redox titrations (using Fe metal) and would like to be selective for the redox species that equilibrate. We would like the sulphate to be ignored by the redox system. I notice that the input files have three ways to include sulphur (as S, S(6) and S(-2)) and assumed that if it was input as S then the redox would be ignored. This does not seem to be the case. We have tried a number of different approaches withou sucess.

Do you have a suggestion. I seem to remember from experiences with MINTEQA2 that you could select those species that you wanted to be ignored for redox calculations.

Here's the user's solution, which is much simpler than my reply.

Thank you for your reply and suggestions. For your information, we found a short cut approach for this example. We did not want the Sulphate to Sulphide reaction to occur in a solution containing Fe metal. We therefore altered the equilibrium constant just for that reaction so that the sulphide stability field was lower than the Hydrogen-Water stability field and the net effect is the absence of sulphate reduction. This seems to be a quick and dirty bandaid approach that can work for specific examples.

Here was my reply:

Redox disequilibrium in PHREEQC applies only to the initial speciation calculation; all reaction calculations include redox equilibrium. Using S or S(6) and S(-2) may affect speciation or saturation indices in the initial calculations, but they only affect the reaction calculations by providing different starting points. Reactions still go to redox equilibrium.

It is a little bit of work with the database to develop redox disequilibrium. Essentially you must define each redox state as a separately named element i.e. "Sulfate", "Sulfide", redefine all aqueous species, and redefine all phases to be written in terms of these new elements. If you do this, then the Sulfate is isolated from the Sulfide just like calcium is isolated from magnesium. I worry that you will eliminate reasonable reactions that occur kinetically with this approach.

I've got a couple items in a FAQ on how to separate redox states in the nitrogen system--https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc If you really need disequilibrium, I would add it to the database but retain the sulfur system as currently defined. Eventually, I think you want to be able to transfer sulfur from the "Sulfate" pool (disequilibrium) to the "S" pool (equilibrium). PHREEQC version 2 allows for kinetic reactions that would allow this type of transformation to occur.

REDOX REACTIONS: I'm trying to simulate the trace metal chemistry of the deep Black Sea. These waters are quite well characterized and are anoxic with total sulfide concentrations around 400 uM. When I use the input files below, forcing equilibrium on the system by allowing some minerals to precipitate, the pe rises to too high a level (and the SO4(2-) data, too), despite several attempts to keep pe low and SO4(2-) non-existent. I tried forcing the pe by specifying a specific pe and according to the S(6)/S(2-) couple (with very high S(2-) and very low (negligible) S(6)). I tried eliminating virtually all of the mineral precipitation (to see if some mineral precip drove the pe up). Finally, I tried increasing the bottom-water Mn(2+) and Fe(2+). In every attempt sulfate and pe always reach unrealistically high values (pe=-4.1 or so and S(6)/S(2-) = 0.25 or so).

My question is simple. Why does sulfate insist on reaching these high levels at equilibrium? Put another way, what drives the initial pe to a substantially higher value at equilibrium? Phrased yet a different way, how can I specify high sulfide and low sulfate concentrations in the initial solution and have them be maintained at equilibrium, after a reaction calculation?

If sulfide is oxidizing to sulfate, something is being reduced. In your case it is C(4) reduced to methane. The calculations indicate that there should be some methane (1e-4) present in the equilibrium solution, so sulfide is oxidized making methane until equilibrium among carbon and sulfur species is attained.

(1) You can add some methane to the initial solution, in which case less sulfide needs to be oxidized. If you have a measurement of methane, you could add it to the solution composition. (2) You can add a little bit of CH2O in the reaction step to reduce carbonate to generate methane (or depending how you look at it, reduce the sulfate back to sulfide). (3) You can "pre-generate" the methane by entering total carbon ["C" instead of C(4)] and using the sulfate/sulfide couple you will generate an appropriate ratio of carbonate/methane in the initial solution. However, I did a little experimenting and the sulfate/sulfide and C(4)/methane boundaries are pretty close together. That means that if you have 10^5 more sulfide than sulfate, you will have a similar preponderance of methane over carbonate.

There are probably other ways to skin the cat, but the key is looking at the oxidation and reduction halves of the reaction.

Check out the comments on your input file below.

 SOLUTION 1
         temp    2
         pH      7.7
         pe      -4.8
         redox   S(-2)/S(6)
         units   Mol/kgw
         C(4)    0.0022  S(-2)/S(6)

(1) The pe and -redox have no effect on the distribution of species for a single redox state [e.g. C(4)]. They do have an effect if a total for a redox element is entered (e.g. "C").

(2) Once -redox is set, it becomes the default, so you don't have to re-enter it for every element.

         P       8e-06   S(-2)/S(6)

(3) Redox choice has no effect on elements that exist in only one valence state (within the database).

         S(-2)   0.004   S(-2)/S(6)
         S(6)    1e-09   S(-2)/S(6)
         N(3)    0.0001  S(-2)/S(6)

(4) You probably mean N(-3), ammonia. This does have an effect on your problem since N(3) nitrite will be reduced to ammonia contributing somewhat to the oxidation of sulfide.

         Mn(2)   5e-06   S(-2)/S(6)
         Fe(2)   2e-08   S(-2)/S(6)
         Ni      6e-09   S(-2)/S(6)
         Cu(1)   3e-10   S(-2)/S(6)
         Zn      6e-10   S(-2)/S(6)
         Cd      1e-11   S(-2)/S(6)
         Pb      4e-10   S(-2)/S(6)
         Na      0.3     as Na+  S(-2)/S(6)
         K       0.0064  as K+   S(-2)/S(6)
         Ca      0.0065  as Ca++ S(-2)/S(6)
         Mg      0.0333  as Mg++ S(-2)/S(6)
         Cl      0.36    as Cl-  charge

(5) Defaults are "as" K, Ca, Mg, and Cl. Since you are actually defining a gram formula weight with the "as", the charge is not necessary and is ignored.

 TITLE S(2-) 10x higher than Black Sea...final pe is?
 
 EQUILIBRIUM_PHASES 1
         Anilite         0.0             0
         Aragonite       0.0             0
         BlaubleiI       0.0             0
         BlaubleiII      0.0             0.0
         Calcite         0.0             0
         Chalcocite      0.0             0.0
         Chalcopyrite    0.0             0
         Covellite       0.0             0
         Djurleite       0.0             0
         Dolomite        0.0             0
         Galena          0.0             0.0
         Greenockite     0.0             0
         Hydroxyapatite  0.0             0
         Magnesite       0.0             0
         Mackinawite     0.0             0
         Millerite       0.0             0.0
         MnHPO4          0.0             0.0
         Pyrite          0.0             0
         Sphalerite      0.0             0

(6) Normally would have an END here, but maybe you just concatenated a couple files. Without the END, SOLUTION 1 and EQUILIBRIUM_PHASES 1 are overwritten by the following data.

 
 SOLUTION 1
         temp    2
         pH      7.7
         pe      -4.8
         redox   S(-2)/S(6)
         units   Mol/kgw
         C(4)    0.0022  S(-2)/S(6)
         P       8e-06   S(-2)/S(6)
         S(-2)   0.004   S(-2)/S(6)
         S(6)    1e-09   S(-2)/S(6)
         N(3)    0.0001  S(-2)/S(6)
         Mn(2)   5e-06   S(-2)/S(6)
         Fe(2)   2e-08   S(-2)/S(6)
         Ni      6e-09   S(-2)/S(6)
         Cu(1)   3e-10   S(-2)/S(6)
         Zn      6e-10   S(-2)/S(6)
         Cd      1e-11   S(-2)/S(6)
         Pb      4e-10   S(-2)/S(6)
         Na      0.3     as Na+  S(-2)/S(6)
         K       0.0064  as K+   S(-2)/S(6)
         Ca      0.0065  as Ca++ S(-2)/S(6)
         Mg      0.0333  as Mg++ S(-2)/S(6)
         Cl      0.36    as Cl-  charge
 TITLE S(2-) 10x higher than Black Sea...final pe is?
 
 EQUILIBRIUM_PHASES 1
         Anilite         0.0             0
         Aragonite       0.0             0
         BlaubleiI       0.0             0
         BlaubleiII      0.0             0.0
         Calcite         0.0             0
         Chalcocite      0.0             0.0
         Chalcopyrite    0.0             0
         Covellite       0.0             0
         Djurleite       0.0             0
         Dolomite        0.0             0
         Galena          0.0             0.0
         Greenockite     0.0             0
         Hydroxyapatite  0.0             0
         Magnesite       0.0             0
         Mackinawite     0.0             0
         Millerite       0.0             0.0
         MnHPO4          0.0             0.0
         Pyrite          0.0             0
         Sphalerite      0.0             0

REDOX SYMBOLISM; SULFIDE and ALKALINITY: I was wondering if you could help me out with a problem I'm having. I am trying to run your program to speciate a well. I entered the sulfur concentration as SO4 for a well that has a pH around 11.4 and pe of -3.849. All the sulfur was in the form of sulfate. To convince myself that this is true and I should not get any sulfides I decided to figure out the equivalent concentration as S-2 and run the well. I thought the sulfide would all be converted to sulfate. Generally speaking, if I convert this, shouldn't my results be the same?

On an Eh/pH diagram, the sulfate dominates at the given parameters. When I tried to run the well with the sulfide equivalent, the program would not run with the same alkalinity valueI had used for the sulfate run.

I then tried to rationalize this by the idea that if I reduce SO4-2 to H2S, there is an increase in alkalinity that creates precipitation of calcite. So I thought if the reverse reaction occurs H2S to SO4-2, calcite would dissolve and increase alkalinity. I had actually messed with the alkalinity on the sulfide run to get the program to run but I decided that I probably shouldn't have to do that.

If you enter the sulfur as sulfate you would use "S(6)" for the element name. All of the sulfur defined this way will necessarily end up in sulfate species in the speciation. Any sulfide would have to be entered as "S(-2)". However, if you enter the total sulfur as "S", then the pe (or redox couple if you use -redox) will be used to speciate the total sulfur between sulfate and sulfide; high pe will put all sulfur in sulfate, pe -4 at pH 11 should put at least some of the sulfur into sulfide, but the transition from sulfate to sulfide is very sharp.

Most of the sulfide is HS- which should be titrated in an alkalinity titration. You must have had more alkalinity from the sulfide species alone than the total Alkalinity that you entered.

I don't know enough about your data, but it sounds like you have measured values for a well water. Generally, total sulfur is not measured; sulfate is the common analyte and occasionally sulfide is measured, but someone has to take some effort to do it. So usually, you would enter the sulfate with S(6) and, if you had it, sulfide with S(-2). Where did the pe come from? pe is a pretty nebulous, but could come from a platinum electrode or a redox couple.

If you have sulfate analysis and pe (from somewhere) there is not an easy way to calculate the sulfide. You can put the sulfate as total sulfur ("S") and see how much ends up in sulfide.

If you are using a measured alkalinity, I think you probably overestimated the amount of sulfide that was in the sample. Even if H2S gas evaded from the sample, the alkalinity should not change. However, the alkalinity would change if some sulfide or carbonate minerals precipitated from the sample.
REDOX AND SURFACE DEFINITION: I am trying to run a PHREEQC simulation that (1) speciates a solution and (2) equilibrates the solution with an iron hydroxide surface. I would like to keep pe after equilibration at -1, but it changes to -3.78. I am wondering why it changes and how I might keep it from changing. I run the following input file with the wateq4f.dat data base.

TITLE Example 1.--Add arsenic and speciate well water SOLUTION 1 MW96-9R sample units ppm pH 7.01 pe -1 temp 10.44 As 620 ppb Ca 120 Mg 103 Na 612 K 140 Fe 33000 ppb Mn 1000 ppb B 1.33 Si 27 as SiO2 Cl 870 charge Alkalinity 2613 as CaCO3 S(6) 2 as SO4 N(-3) 417 END SURFACE 1 Hfo_w 0.105 600 45 SURFACE_SPECIES # used to stop Ca and Mg use of sites Hfo_wOH + Mg+2 = Hfo_wOMg+ + H+ log_k -15. Hfo_wOH + Ca+2 = Hfo_wOCa+ + H+ log_k -15. USE SURFACE 1 USE SOLUTION 1 PHASES # Fix_e- # e- = e- # log_k 0.0 Fix_H+ H+ = H+ log_k 0.0 # Fix_As # H3AsO3 = H3AsO3 # log_k 0.0 EQUILIBRIUM PHASES 1 # Fix_As -5.082 Fix_H+ -7.01 NaOH 10 # Fix_e- -1 END

What is magic about pe -1? Rather than look just at the pe, look at the amounts of each valence in redox couples. With your problem you start with As(3), C(4), Fe(2), Mn(2), N(-3), and S(6). After reaction you've got As(3), 14 mmol of C(-4)?, 50 C(4), Fe(2), Mn(2), 27 N(-3), 2 N(0), and S(-2). The main question is how did we make 14 mmol of C(-4)? If C is reduced, something must be oxidized. The answer is pretty subtle. The surface you defined was Hfo_w, uncharged.
```
Hfo_w + H2O = Hfo_wOH  + 1/2 H2
4H2 + CO2 = CH4 + 2H2O
       
```
So the first thing is to define the surface as the uncharged form "Hfo_wOH", since Hfo_w alone should have a charge of +1.

Once you do this, the pe is still -3.6, but you no longer have such a large amount of methane. If you go through the same process as above, you find that some of the ammonium has oxidized to N2(aq) generating a couple of millimoles of methane. So not worrying about what the pe is, does this amount of methane or N2(aq) bother you? If it does, there are a couple things you can do.

(1) You can define a new element Amm so that ammonia can not oxidize, see below. Enter the ammonium as "Amm" in the solution input. This way you will not generate either methane or N2(aq). (You would have to add some more input to allow Amm sorption (SURFACE_SPECIES) or exchange (EXCHANGE_SPECIES).

(2) If the N2(aq) is ok, but you don't like the methane, you can add an oxidant to generate the -1 pe just like you did for the pH. You have it roughed in in your input file, simply uncomment the phases and include the following in the equilibrium_phases.
```
	Fix_e-    -1	O2	10
       
```
Here are the basics for making ammonia a separate element.
```
SOLUTION_MASTER_SPECIES
	Amm		AmmH+	0.0	AmmH		18.0
SOLUTION_SPECIES
AmmH+ = AmmH+
	log_k	0.0
	-gamma	2.5	0.0
AmmH+ = Amm + H+
	log_k	-9.252
	delta_h 12.48	kcal
AmmH+ + SO4-2 = AmmHSO4-
	log_k	1.11
	
```
REQUIRED SPECIES IN DATABASE: I am trying to use a database management program I used for PHREEQE to write files for PHREEQC. When I try to run PHREEQC using these files, I get an error message that "H2(aq) [and other things] not defined in solution_species". The difficulty is that H2(aq) [and the other things] are not in the thermodynamic data file read by PHREEQC and I have no clue where they are coming from, so it would be most helpful if there was some way to see exactly what PHREEQC is reading.

Though not well documented, a minimum set of elements and aqueous species must be defined. SOLUTIION_MASTER_SPECIES must contain H, H(0), O, and O(0). SOLUTION_SPECIES must contain H+, H2, H2O, O2, e-, and OH-.
THERMODYNAMIC DATA FOR Sb (Antimony): I'd like to add Sb species data to the database because I need to model Sb and there are no such data.

By the way, is it possible to get (download?) an annotated version of the PHREEQC database (or WATEQ4F), complete with references with the sources of the data?

For the most part, the data in phreeqc.dat and wateq4f.dat are taken from the wateq4f documentation (USGS Open-File Report 91-183) and a paper by Nordstrom and others, 1990 ACS volume Chemical Modeling of Aqueous Systems II. There are a few numbers from other sources that may not be documented. Every once in a while, Jim Ball called and said I had the wrong number (which was in WATEQ4F) and asked me to put in a different one.

Obviously, I should go through and add the documentation as comments in the data file. I've tried to avoid getting caught up in thermo data and haven't done it. Unfortunately, no one else has paid much attention to the thermo data either.

I'm a little reluctant to start putting out new thermodynamic database files unless they are published. I felt ok about translating the MINTEQ, PHREEQE, and WATEQ4F database files, because I didn't have to make any decisions (except to try to fix MINTEQ where I could) and could point people to other sources. I would consider adding a "contrib" directory to my web page making it available, but also making it clear that the USGS and I are not responsible for what is in there.

Additional information provided by the questioner: The EQ3/6 database has a lot of Sb data and is the best possible source that I'm aware of, although I don't have any documentation (yet) for their data. The guardian of the EQ3 database, I am told, is James W. Johnson at jwjohnson@llnl.gov. I've sent him an e-mail already to try to track down the source of the data so that it doesn't just get added blindly.

I'm not sure if the EQ3 database is still on the net at LLNL. But it exists in adopted form in the CHESS database which is on-line at http://www.cig.ensmp.fr/~vanderlee/chess/index.html.
ALKALINITY & TDIC: I'm working with some very dilute groundwaters with low pH values (4.8 - 5.0) for which there is no measured alkalinity or carbon data. Since the pH is close to the extended alkalinity endpoint I guess the carbon concentrations are below detection limits. Is it possible to set up PHREEQC to calculate an alkalinity value, assuming all of the alkalinity is due to carbon species?

The alkalinity is probably very small. However, the total carbon (virtually all dissolved CO2) may be significant, depending on what you are modeling. I would adjust the total carbon to be in equilibrium with a soil PCO2, maybe log PCO2 of -2. to -3 atm, in the SOLUTION input, as follows:
```
C	1	CO2(g)	-2.5	
	
```
ELECTRICAL BALANCE: A quick 2 questions: 1) Does PHREEQC calculate the electrical balance using the speciated solution, rather than the just the converted (eq/kg) "raw" input data? In other words do neutral species, such as H2CO3, not partake in the neutrality calculation? 2) Does PHREEQC use the following to

(1) You are correct. The charge balance is sum(m(i)*z(i)) for all species in the distribution of species.

(2) No. It is simply sum(m(i)*z(i)), which is equivalents. It is not normalized by the sum of cations and anions, nor is it multiplied by 100. In the next version 2 (>2.0.30), I've added the percent error to the printout.
SPECIFIED pH IN REACTION: One possible upgrade that I would appreciate with PHREEQC is a way to titrate a solution to a specified pH. (or can this be done presently?)

It's a little cryptic, but example 8 in the manual for version 1 calculates compositions at fixed pH values. First you must define a phase as follows:
```
PHASES
        Fix_H+
        H+ = H+
        log_k  0.0
	
```
Then you can include the following:
```
SOLUTION
   your solution composition here
EQUILIBRIUM_PHASES 1
        Fix_H+   -4.5   HCl    10.0
END
	
```
One reaction solution will be calculated that adds enough HCl to obtain a pH of 4.5.
SPECIFIED ACTIVITY: 1. Can you draw log C vs. pH or logC vs. pe? This is can be down in Mineql by fixing pH or pe and calculating the concentrations of various species in the solution.

2. How to fix activity (not the total concentration, Ctot) of a component? For example, I want to find in a solution with HCO3-=2E-3, how much CO3(total) to put in.

You can't directly fix an activity quantity in PHREEQC, there must be a reaction that causes that quantity to stay constant. Example 8 in the manual has surface sorption as a function of pH.
You can use a similar strategy for HCO3. The following would make a solution with activity of 2e-3.
```
SOLUTION 1
# pure water
PHASES
Fix_HCO3
HCO3- = HCO3-
log_k 0
EQUILIBRIUM_PHASES
# -2.699 is log10(2e-3)
Fix_HCO3     -2.699   NaHCO3   10
END
	
```
UNREACTIVE MINERALS: How do you ignore certain minerals (such as dolomite) so that they do not participate in the calculation?

Unless you include dolomite in the reactive phase assemblage (EQUILIBRIUM_PHASES), it will not react. You must select which phases are reactive. All are entered through EQUILIBRIUM_PHASES keyword data block. Those present have positive numbers of moles. Those that are not present initially, but could precipitate have zero for the number of moles.
CO-PRECIPITATION/SOLID SOLUTIONS: Correct me if I'm wrong, but I don't think PHREEQC directly considers co-precipitation (e.g. coprecipitation of trace metals by FeS in anoxic waters). Can you think of a way to pseudo-simulate this (adsorption to a precipitated FeS, some sort of ion exchange, etc.)?

Version 1.6 does not consider coprecipitation. Version 2 does consider solid solutions. It is possible to model nonideal solid solutions with two components or ideal with any number of components. The main drawback in the formulation is that complete thermodynamic equilibrium is assumed. At each point in reaction the entire solid is in equilibrium with the solution. An onion-skin type model that allows equilibrium between the surface and the solution may be more reasonable.
ANOTHER SOURCE OF THERMODYNAMIC DATA: For Minteqa2, Dr. Willard Lindsay made many additions to the databases that made them more pertinent to geochemical solutions in soils. (His changes can be downloaded from http://www.colostate.edu/Depts/SoilCrop/mnteq.html). Has anyone incorporated these additions into a Minteq database in the Phreeqc format? If no one has, then I will probably do so, but i hate to duplicate effort that has already been made.

I translated the MINTEQA2 database to PHREEQC format. The file "minteq.dat" in the distribution is my best guess at the translation. I found several (maybe 20) errors in the reactions, but I did not do any checking of the actual log K data. I have not done anything with an updated version of the MINTEQ databases, but would welcome any additions or modifications to the minteq.dat file. I haven't heard of anyone else pursuing this work, though there are frequent requests for expanded databases.
I am learning Java and I was considering writing a Java front-end to Phreeqc for teaching purposes, but I have heard that Dr. Parkhurst is working on a Java version of Phreeq. Is this a Java front-end or a complete rewrite of the program? What is the timetable for this tool? Are alpha or beta testers needed?

A Visual Basic version called PhreeqcI is available for DOS for version 1.6. It is available at the web page https://wwwbrr.cr.usgs.gov/projects/GWC_coupled.

A revised version 2 of PHREEQC is in review and has many additional capabilities, most notably reaction kinetics. We considered writing a JAVA interface for version 2, but got bogged down in undocumented features and by the need for other features that were only in beta release. Though a JAVA version was preferable, we dropped back to Visual C++. This interface is in the initial phases, but we have pretty well settled on the overall look and feel. The main screen that is a tree of keywords, file handling, line editor, the SOLUTION keyword, running, and viewing output should be at the alpha stage in the near future. I'll try to put it on the net for comments.
SORPTION WITH FREUNDLICH AND LANGMUIR ISOTHERMS: Can I model sorption according to Freundlich or Langmuir isotherms with PHREEQC?

Yes. The derivation and examples are given here. The Freundlich equation is:
```
     q = Kf * C^n                                           (1)
	
```
For PHREEQC, The mass-action equation is derived from the chemical reaction equation that defines a species. The following surface complexation reaction generates the correct mass-action equation for the Freundlich equation:
```
     Sites + n * C = SitesC                                     (2)
	
```
The mass-action equation for reaction (2) is:
```
     K = [SitesC] / ([Sites] * [C]^n)                           (3)
	
```
The brackets indicate activity, which for a sorbed species is the fraction of sites the sorbed species occupies. Now q = m(SitesC), where m(SitesC) is the number of moles of C that is sorbed. [SitesC] = m(SitesC)/TOT(Sites), [Sites] = m(Sites)/TOT(Sites), m(Sites) is the number of moles of unoccupied sites, and TOT(Sites) is the total number of sorption sites. Substituting into equation (3) gives the following equation:
```
     K = {q/TOT(Sites)} / ({m(Sites)/TOT(Sites)} * [C]^n)         (4)
	
```
Canceling TOT(Sites) and rearranging 4 gives:
```
     q = (K * m(Sites)) * C^n.                                (5)
	
```
Equation (1) and (5) are identical when K = Kf / m(Sites). The trick is to keep m(Sites) (the number of unoccupied sites) constant throughout the calculations. This can be arranged by making TOT(Sites) large relative to the amount of C that sorbs. In that case, the unoccupied sites, m(Sites), will stay nearly equal to the total number of sites, TOT(Sites). The value for the association constant of the SURFACE_SPECIES is then K = Kf / TOT(Sites).

Note in equation 2 that the mass-action coefficient for C is n, but for SitesC it is 1. This equation is not balanced in C. PHREEQC-2 allows unbalanced equations by defining SURFACE_SPECIES with the option -no_check, which disables the element- and charge-balance checking of an equation. However, when an unbalanced equation is used for mass-action, it is necessary to define the explicit stoichiometry of the product with the option -mole_balance. In this case, the option should be used as follows:
```
	-mole_balance SitesC 
	
```
The Langmuir equation is:
```
     q = Smax * C / (Kl + C).                               (7)
	
```
The equation can written as:
```
     q = (Smax - q) * C / Kl.                               (8)
	
```
This is the mass action equation for:
```
     S + C = q;     K,
	
```
since S = (Smax - q), and when K = 1/Kl. (Notice that mole fractions are used for the activities of the surface species in PHREEQC-2).

The following input set for PHREEQC version 2 has two pollutants: Polf sorbs according to a Freundlich isotherm and Poll sorbs according to a Langmuir isotherm. The parameters for the Freundlich isotherm are Kf = 10, n = 0.8, and the parameters for the Langmuir isotherm are Smax = 30 and Kl = 2. The input file prints a selected output file with 6 columns. Column 1 is the dissolved concentration of Polf; 2 the sorbed concentration of Polf; 3 the amount of Polf sorbed as calculated from the dissolved concentration and the Freundlich isotherm; 4 is the dissolved concentration of Poll; 2 the sorbed concentration of Poll; 3 the amount of Poll sorbed as calculated from the dissolved concentration and the Langmuir isotherm. Columns 2 and 3 are equal and columns 5 and 6 are equal, which indicates the Freundlich and Langmuir isotherms are being calculated correctly.
```
SOLUTION_MASTER_SPECIES
        Polf  Polf    0.0     Polf    1.0
        Poll  Poll    0.0     Poll    1.0
SOLUTION_SPECIES
        Polf = Polf;  log_k 0.0
        Poll = Poll;  log_k 0.0
SURFACE_MASTER_SPECIES
        Sites  Sites
        Smax Smax
SURFACE_SPECIES
        Sites = Sites
                log_k 0.0
        Smax = Smax
                log_k 0.0
# Freundlich: SitesPolf = Kf * Polf^0.8
        Sites + 0.8Polf = SitesPolf   
                -no_check
                -mole_balance SitesPolf
                log_k -99.0             # log ((Kf = 10) / TOT(Sites))
# Langmuir: SmaxPoll = (tot_Smax - SmaxPoll) * Poll / Kl
        Smax + Poll = SmaxPoll    
                log_k -0.30103          # log (1 / (Kl = 2))
END
SOLUTION 1
	-units mmol/kgw
        Polf 1e3                  # All concentrations 1 mol/l
        Poll 1e3
SURFACE 1
        Sites 1e100 1.0 1e100
        Smax 30   1.0 30
                -equil 1
                -no_edl true
REACTION 1 Removes 11 moles of Polf and Poll in 5 steps
        Polf -1.0  Poll -1.0
        11 in 5 steps
SELECTED_OUTPUT
        -file freundl.sel
        -reset false
USER_PUNCH
        -heading diss_Polf_  sorb_Polf_ _q_Freund_ diss_Poll_  sorb_Poll_ _q_Lang___
10 Kf = 10
20 n = 0.8
30 punch mol("Polf"), mol("SitesPolf"), Kf*mol("Polf")^n
40 Kl = 2
50 Smax = 30
60 punch mol("Poll"), mol("SmaxPoll"), Smax*mol("Poll")/(Kl + mol("Poll"))
END
	
```
CONVERTING SOLID CONCENTRATIONS FOR PHREEQC: I am modeling the surface complexation of various arsenic species. I am using the Dzomback and Morel HFo_w and HFo_s parameters, geothite and gibbsite parameters from Manning and Goldberg, and other surfaces that I have found data for. I am trying to see how well the model matches partial extraction data for stream sediments. My problem lies in interpretation of the molal concentration of surface species, for instance, Hfo_wOHAsO4-3.

Here is the rub, for some simulations I have assumed that the 1 kg (same as 1 L at my ionic strength) of water is in a .3 porosity system with the sediments having a density of 2.5. This results in the 1 kg or 1 L of water being in contact with 3.3 L or 8.3 kg of sediment. Now, I assume that like everything else in the model, and the definition of molal, that the concentration reported by the model for Hfo_wOHAsO4-3 is with respect to the 1 kg of water. This in reality is moles sorbed per 8.3 kg of seds. Correct? So that to compare the reported molal sediment concentrations to my observed partial extraction data (mol/kg) I need to divide the PHREEQC results by 8.3? Or have I really missed the concept somewhere?
The number of surface sites is defined as moles, not moles/kg(water), but since there is usually about a kilogram of water it is numerically about the same. Think of it as a surface that is defined in moles and put in a beaker. A solution is then added that usually is about a liter (kg water), but not required. A 0.5-kilogram solution could be added, in which case the moles of surface would be the same, but the molality of the species would be double what you would have if 1 kilogram of solution were added.

.7 L(rock)/.3L(water) = 2.33 L(rock)/L(water)

2.33L(rock)*2.5 kg(rock)/L(rock) = 5.825 kg(rock)/L(water)

So the number you want is probably (moles of Hfo_wOHAsO4-3)/5.825.

A more direct way to get at this number is to use the number of sites per gram of sediment. Then the number of sites you use divided by the number of sites per gram of rock gives you directly the number of grams of rock. So moles sorbed divided by grams of rock gives you the correct number. This way you don't have to work through the volume of water at all.
CONVERTING SOLID CONCENTRATIONS FOR PHREEQC: I am going to use the default amount of water 1 kg. 30% porosity resulting in 2.33 L of rock in contact with 1 L of water. Since I spreadsheet the following calculations, and am going to vary porosity as a fitting variable, I think that I need to keep the water/rock ratio in for now. I also am working with 1 L of water to limit the number of brain tremors that I have.

For example, I take the laboratory determination of amorphous FeO surface, a number in g/kg of sediment. Since I want to react 1l of water with the FeO For example, I take the laboratory determination of amorphous FeO surface, a number in g/kg of sediment. Since I want to react 1l of water with the FeO surface on 5.825 kg of seds (your correct number) I determine the grams of surface in contact with that 1 L of water as FeO(g/kg)* 5.835 = 16.9 and use this in the mass of surface block in PHREEQC. I convert this to moles of FeO surface available (.30) to the 1 L water and then apply the .005 mole-sites/mole Fe correction from Dzomback and Morel, resulting in .0015 moles of FeO sites available to the 1 L water. This is input into the surface sites (moles) block. I use default D&M surface area of 600 m2/g sorbant. I then allow the surface to complex with many pore volumes of the solution (TRANSPORT) to emulate a constant concentration and flow through saturated sediments. Bingo! Out comes a final molal concentrations of several nasty ole' arsenic species on the sites. I sum the molal concentrations to get the total molal arsenic sorbed on the seds. But wait a second, the molal concentration is defined w.r.t. a sediment total mass of 5.825 kg. So to get to the "real" molal arsenic concentration (for comparison to the laboratory measured molal concentration) I divide the PHREEQC sum of molals by 5.825.

Whaddaya think?

I'd keep the mass of water constant too. However, as you change porosity, the volume of media needed to contain a liter of water gets bigger.

If you change porosity, you start all over with a different number of kg of sediment [(1-Por)/(Por)]*2.5. Lower porosity greater mass of sediment.

D&M use gram formula weight of about 100 g/mol (check to be sure) for hydrous ferric oxides, which might screw things up a little. Your FeO is about 72 g/mol. Seems like the 600 m2/g would be off by 30%, but may not be a big effect. You may want to inflate your g/kg sediment by 100/72 to get g Hfo/kg sediment.

I prefer moles, but molal should equal moles.
LEACHING URANIUM: What I'd like to really do is to simulate the leaching of Fe-oxide (containing U) by DI water. In the latest experiments, we actually leached contaminated soil with DI water and the solution1 concentrations in the file I sent were the final solution concentrations. What I had planned to do was this.

1. Add U to pure water as UNO3 to about 1000 ppm (simulate High U waste stream)

2. allow this solution to equilibrate with the Fe-oxides to put the U on the solids

3. allow Fe-oxides in (2) to equilibrate with pure water in the presence of calcite.

As an alternative: From ICP I know the soil is about 1% Fe or about 1.9% Fe(OH)3. The U concentration in the soil is about 250 ppm with most of this (from fission track work) in the Fe-oxides. In the experiments we had 20 g soil and 200g water--thus on a per liter basis there is 1.9 g Fe(OH)3 of 0.001M U concentration in contact with the DI water. How do I fix the initial U concentration in the Fe(OH)3 before I equilibrate it with the DI water?

Here's an outline:
```
SOLUTION 1
#  high uranium water composition
END
SURFACE 1
	Hfo_s 	sites 	specific_area 	grams
# this next line gives you the composition of the surface
# that is in equilibrium with solution 1,
# but does not allow solution 1 to vary
# Adjusting U in solution 1 will put
# different amount of U on surface 1
	-equil 	1
END
SOLUTION 2
# DI water here
END
USE surface 1
USE solution 2
END
     
```
To flush 10 pore volumes you can use TRANSPORT (version 1.6 or ADVECTION in version 2) replace starting at solution 2:
```
SOLUTION 0
# DI water	
TRANSPORT
	-cells 1
	-shifts 10
END
    
```
SELECTED OUTPUT: In example 8, specifically how can I produce from the SELECTED_OUTPUT keyword the Figure 4 of the user's guide?

When example 8 is run, a file ex8.pun is produced. Columns of data are printed to ex8.pun for the following items:
```
1. sim 
2. stat 
3. soln 
4. exch 
5. surf 
6. pure  
7. gas step           
8. temp             
9. ph             
10. pe 
11. ionic_strength            
12. H2O         
13. step_x         
14. Zn+2_m    
15. Hfo_wOZn+_m    
16. Hfo_sOZn+_m
      
```
The SELECTED_OUTPUT keyword causes the file ex8.pun to be printed with 13 default columns (1-13 above) and 3 additional columns as specified by the -molalities identifier. The extra 3 columns are the molalities of aqueous Zn2+ and molality of zinc sorbed on weak surface sites (Hfo_wOZn+) and stong surface sites (Hfo_sOZn+)

From the input file, the first 13 simulations including and following the selected_output keyword use solution 1 (Zn .0001 mmol/kgw, 1e-7 m) and fix the pH in increments between 5 and 8. The first 13 lines of ex8.pun are used to generate the top panel of figure 4. X value is pH taken from column 9 and y values are weak-site sorption, column 15, and strong-site sorption, column 16.

Lines 14-26 are printed for simulations that use solution 2 (Zn .1 mmol/kgw, 1e-4 m). Data from these lines are used to generate the lower panel in figure 4. Again pH is in column 9 and weak- and strong-site sorbed zinc are in column 15 and 16.
DOWNLOADING PROBLEMS WITH NETSCAPE: I went to the site you indicated: (https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/), and tried to download the DOS beta version, but what happens is that I just get a new screen with a large collection of ASCII binhex characters. When I try to download the Mac version I get a file and progress dialog box, but the DOS files come across as ASCII characters.

Is there a way to get the package as a file that I can then transfer to the PC at home and unpack??

Sounds like you are using Netscape. I think you need to click the right button and pick "Save link as ..." from the menu.
SORPTION SITES: I was wondering about the actual concentration of U in the Fe-Oxide from my simulations. I entered 1 g, 0.01 mol sites, 600 m2/g area and have about 1x10-4 moles of U on the solid. Is this 10-4 moles U in 1 gram (e.g. 238.029 g/mol x 1,0000 mg/g x 10-4) = 23 ppm. or is this divided by the .01 mole of sites? The ICP data shown several thousand ppm U in the solid heavy mineral-oxide phase.

You have 1e-4 moles of U sorbed on sorption sites and the total number of sorption sites is 0.01 m. Now, I'm not sure about the concentration in the solid. You have defined a gram of solid, but is that the total amount of solid, or the iron oxide portion of the solid? Anyway, 1e-4*238. = .023 g of Uranium is sorbed. You need to divide by the total mass of solid in g times 1e6 to get ppm. If there is only 1 g solid, then it's 23,000 ppm or 2.3 percent.

If I did the calculation right, with 1 g of solid, you have 23,000 ppm U in the solid.
EXAMPLE 9 IS REALLY SLOW ON A MAC: I ran the ex9 transport problem on my PPC mac (Mac 8500) and it took approximately 1/2 hour to solve.

Try using "-status false" in a PRINT keyword data block. The status line is written many many times to the screen. I think Macs do not buffer the screen output so the program waits each time for the status line to be written and it seems to take a very long time.
EXAMPLE 10 BOMBED: Ex10 bombed complaining that it didn't have data for As.

Example 10 requires the wateq4f.dat database file which contains definitions for arsenic. In PHREEQCI, the database file is changed with the "Options" menu. In batch mode, the third argument to phreeqc is the name of the data base file
```
phreeqc  input_file  output_file  data_base_file
      
```
ALKALINITY: Does PHREEQC calculate carbonate alkalinity after it 'corrects' for silica and other ions that act as alkalinity in a titration?

Yes. The alkalinity that is entered as an input parameter is assumed to include carbonate species and all other species that would be titrated in an alkalinity titration.

PHREEQC calculates alkalinity by a sum(alk(i)*m(i)), that is sum of alkalinity contribution of each aqueous species times its molality. For example, H3SiO4- has an alkalinity contribution of 1. So alkalinity is not simply carbonate alkalinity, but the sum of all aqueous species that will accept a proton in an alkalinity titration.

One possibility that sometimes occurs is that the alkalinity you specify is less than the alkalinity of silica, borate, or other noncarbonate species. The program is trying to adjust the carbon species to sum account for the difference between the noncarbonate alkalinity and the total alkalinity. If this difference is negative, no mathematical solution exists in which carbon concentrations are positive. You'll get a nasty-gram about Noncarbonate alkalinity exceeding total alkalinity.
ALKALINTITY: If I allow alkalinity to charge balance is it only balancing carbonate species?

You can't balance with alkalinity, it isn't allowed by in the program. You can balance with carbon [C or C(4)], in which case, the total carbon in solution will be adjusted such that the charge on the carbon species is exactly what is needed for charge balance.
EH, REDOX POTENTIAL: My Eh measurments are mV from a platinum electrode that is referenced to an Ag/AgCl pair (in the pH electrode). The Orion meter, according to Orion, does not correct the display value for the reference volatage Sooo, how do I correct for the reference electrode?

Eh system (input to PHREEQC)= EMF measured - Eh reference.

I think Ag/AgCl has a potential of about 220 mV, but depends on temperature. You should subtract the half-cell potential of the reference from the measurement to get the Eh relative to the hydrogen half cell, which is what should be reported and used in PHREEQC. Eh in air saturated water is not well poised, but is usually in the 300 - 500 mV range.
CONCENTRATION UNITS, MEQ/L: A student of mine is trying to model water concentrations given in the literature as meq/L, but the program is not converting to molality. What is wrong?

It's a bug or perhaps a feature. The program only allows equivalents for the definition of alkalinity, no other element concentrations can be input as equivalents, it defaults back to mmol/L for everything but alkalinity. You'll have to convert to mmol/L for calcium and other non-monovalent elements.
LEAKAGE IN 1D TRANSPORT SIMULATIONS: I'm a phd-student working on a semi-confined aquifer. I'm studying the relation between groundwater dynamics and the resulting chemical variation in the groundwater quality. I want to use PHREEQC for 1D transport coupled to chemical reactions as dissolution / precipitation, redox reactions and cation exchange. I use the version 2.0 so I can consider dispersion and diffusion.

The problem is: along the flowpath water is lost through the overlaying confining unit. As I am trying to model freshening of a once saline aquifer I can not match the chlorine concentrations in the downstream part of the system which rise fairly steeply. When I use the transport keyword, the total volume of water is shifted from cell i to cell i+1. In this way the solutes are transported too fast out of the column, especially in the downstream part where leakage through the overlaying confining unit becomes important and where flow velocities are as 10 times lower than in the recharge area.

I tried to use extremely big dispersion (1000m for cells of 150m) in the downstream part but freshening is to fast (after three pore volumes all Cl has disappeared).

Is it possible to shift only a fraction of the volume from cell i to cell i+1 in the downstream part ? And can this fraction be changed from cell to cell (from big to small) ?

PHREEQC is not really set up to handle this sort of problem easily. However, I think you should be able to come close by carefully defining some MIX keywords for PHREEQC version 2.0. The following example effectively transports only 20% of the solution from cell to cell. It makes use of stagnant zone cells, which are numbered n+2 to n+1+n (where n is the number of cells in the column) essentially as a means of storing information about the cell.

In the case of cell 1, this is what happens: (1) transport moves solution 0 into cell 1, (2) the "MIX 1" keyword causes cell 1 to be mixed with cell 7 (which is stagnant) to make a new composition for cell 1, (3) before this new composition is saved for cell 1, "MIX 7" causes cell 1 (that is solution 0 that has been transported in) to be mixed with cell 7 to generate a new composition for cell 7. Thus, as solution 0 is transported into cell 1 at the next time step, cell 7 contains composition of cell 1 from the previous time step to be used for mixing. Kind of messy, but may be of some use to you. There may be other ways to use MIX to get the desired result. Any minerals or other reactants should be defined the same for each column cell and its associated stagnant cell.
```
SOLUTION 1-20
	units mol/kgw
	Na	1
	Cl	1
END
SOLUTION 0
	units mmol/kgw
	Na	1
	Cl	1
END
MIX 1
	1  .2
	7  .8
MIX 2
	2  .2
	8  .8
MIX 3  
	3  .2
	9  .8
MIX 4
	4  .2
	10 .8
MIX 5
	5  .2
	11 .8
MIX 7
	1  .2
	7  .8
MIX 8
	2  .2
	8  .8
MIX 9  
	3  .2
	9  .8
MIX 10
	4  .2
	10 .8
MIX 11
	5  .2
	11 .8
TRANSPORT
	-cells    5
	-shifts   5
	-length	  1
	-disp	  0
	-diffc    0
	-stag     1
SELECTED_OUTPUT
	-reset false
	-solution
USER_PUNCH
20 punch TOT("Na")
30 punch TOT("Cl")
END	
       
```
INVERSE MODELING WITH GASES: I am modeling denitrification along a flowpath and examining pyrite oxidation as a possible mechanism. When using netpath for the simulation I use O2(g) as a phase (among others). My question is that when output states that O2(g) is lost from solution, would this include O2 which is reduced to O(-2) or just degassing?

The mole transfer of O2(g) is the amount of oxygen that must enter or leave the solution. If, for example, the only reaction is the oxidation of organic matter with dissolved oxygen, and not all of the dissolved oxygen is consumed, the reaction would have no mole transfer of oxygen gas.
MODELING PYRITE OXIDATION BY OXYGEN: If I use the reactions from netpath, which require oxygen to enter solution, and do forward modeling using PHREEQC is it necessary then to have the removal of O2 as one of my reactions?

You must include a source of O2 in the forward modeling, either with REACTION keyword datablock that has O2 as a reactant or with EQUILIBRIUM_PHASES that specifies a fixed partial pressure of O2(g). If you use REACTION, the amount of O2 should be the mole transfer calculated by NETPATH. If you use EQUILIBRIUM_PHASES, an appropriate partial pressure needs to be specified.

Note that (depending on how you do the forward modeling) charge balance errors make it possible that PHREEQC forward modeling does not match the final solution pH or other solution composition data, so check whether the results of PHREEQC match the original data.

INVERSE MODELING WITH REDOX: Would you happen to have an example input/output for a Phreeqc Inverse calculation involving redox reactions?

Here's an example that is similar to a NETPATH example for a coastal plain aquifer. It doesn't have a lot of redox reactions, but it does have the input necessary for modeling organic decomposition and iron reactions.

#
#   Coastal Plain, problem 3
#
SOLUTION 1 
	units   mg/l
	temp	17
	pH 	7.6
	Ca	41
	Mg	13
	Na	3.8
	Alkalinity	178 as HCO3
	S(6)	19
SOLUTION 2 
	units   mg/l
	temp	19
	pH 	8.4
	Ca	3
	Mg	1.9
	Na	140
	Alkalinity	367 as HCO3
	S(6)	13

INVERSE_MODELING 1
	-solutions 1 2
	-uncertainty 0.03
	-phases
		CH2O
		Calcite
		Goethite
		Pyrite
		CaX2
		MgX2
		NaX
		Gypsum
PHASES
	CH2O
	CH2O + H2O = CO2 + 4H+ + 4e-
	-log_k	0.0

END

SURFACE COMPLEXATION WITH EXPLICIT DIFFUSE LAYER: Why am I getting a "Negative sum in g_function" nastygram?

If you use the diffuse_layer calculation, you should charge balance all your solutions. Use "charge" with one of the more inert ions in each solution.
PE, EH, REDOX POTENTIAL: Most field redox measurements I see are expressed as millivolts. I've tried to find a formula for how to convert redox measurements expressed as millivolts to a pe value, but have not had any success. Is there a conversion formula for this?

pe = 16.9*Eh(Volts) at 25C or, pe = .0169*Eh(mV)

more generally,
```
pe = F/(2.303RT)*Eh(Volts)

F = 96.5 kJ/(V-eq)
R = .0083147 kJ/(deg-mol)
T = Kelvin
2.303 ~ log(10)
        
```
TRANSPORT SIMULATIONS: My computer is not fast enough to run transport simulation with PHREEQC-2. 120 cells with 120 stagnant cells takes hours to run.

If you are using a pentium II processor, Normally this size simulation doesn't take that long unless (1) you are using -diffuse_layer surface calculations or (2) you have a very large dispersivity relative to cell size.
SATURATION INDICES: I want to calculate the saturation index with respect to calcite for great amount of samples (more than 1000). Using Netscape Database, I can do this for single chemical analysis. What may be a procedure for simultaneous calculation for many samples?

The easiest way would be to use the beta version of PHREEQC 2. It has a new keyword called SOLUTION_SPREAD that allows each analysis to be entered in a row. Headings of columns indicate the quantities being entered. Headings and data values are tab delimited. A subheading line allows for all of the optional input that follows concentration values in SOLUTION data block. Columns for column headings that are not understood are ignored. You should be able to cut columns out of a spread-sheet and paste them into an ascii file, add the column headings as necessary, and run.

Beta test of Version 2 and description of input is available at https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc
```
SOLUTION_SPREAD
Number	pH	Si	Ca	Mg	Na	K	Alkalinity	S(6)	Cl
							gfw 50			
1	6.2	0.273	0.078	0.029	0.134	0.028	0.328      	0.01	0.014
2	6.8	0.41	0.26	0.071	0.259	0.04	0.895    	0.025	0.03
        
```
LINEAR Kd AND pH-DEPENDENT SORPTION: I am trying to write some PHREEQC tranpsort code that will adsorb HAsO4-2 dependent on the pH of the cell solution. I would like to use a linear equation for Kd and pH with following form: Kd = -.175*log(H+)+1.925.

Surface sorption gets pretty complex, and Kd's are often not the best way to go, because they apply over such a narrow range of solution composition. I would be more comfortable if a surface-complexation Log K were fit rather than a Kd, but it is much more work.

The key is defining the right mass-action equation for your Kd. Assuming you are sorbing HAsO4-2 and you want:
```
Kd = [HAsO4-2sorbed]/[HAsO4-2], 
        
```
then I think the following chemical reaction will give you the desired mass-action equation:
```
HAsO4-2 + Sor = SorHAsO4-2 + .175 H+
        
```
The mass-action equation is as follows:
```
K = ([SorHAsO4-2][H+]^.175)/([Sor][HAsO4-2])
        
```
or
```
LogK = log[SorHAsO4-2] + .175log[H+] - log[Sor] - log[HAsO4-2] 
        
```
rearranging gives the following:
```
(LogK + log[Sor]) - .175log[H+] = log[SorH2AsO4-] - log[HAsO4-2]
        
```
So if you define Sor = 1e100, and you want log Kd = log(1.925), then the log K for the reaction should be
```
log(K) = -100 + log(1.925) = -99.72
        
```
Here are the surface data blocks that define the the pH dependent Kd reactions:
```
SURFACE_MASTER_SPECIES
    Sor Sor
SURFACE_SPECIES
Sor = Sor 
    log_k 0.0
Sor + HAsO4-2 = SorHAsO4-2 + 0.175H+ 
    log_k -99.72
    -no_check
    -mole_balance SorHAsO4
SURFACE 1-15 linear arsenic adsoprtion
    Sor 1e100 1.0 1e100
     
```
KINETIC RATE EXPRESSIONS FOR PRECIPITATION: I tried to specify kinetic rate laws for the precipitation of secondary uranyl phases by setting moles (m)=0. It did not seem to have any effect on the results that I obtained compared to equilibrium precipitation of these phases with no precip rate law. Can PHC use a precipitation rate law?

Precipitation can be modeled, but you need to know the sign convention. If you SAVE the variable "moles" in the rate expression (defined in RATES) and the reaction is UO2 with a coefficient of "coef" (defined in KINETICS), then the abs( moles*coef) is the amount of UO2 that reacts. However, the sign of moles*coef determines dissolution or precipitation. If the sign is positive, the reaction increases solution concentrations (mineral dissolution); if the sign is negative, the reaction decreases solution concentrations (mineral precipitation). You can not dissolve if moles of mineral is equal to zero (m = 0), but you can precipitate. The new amount of kinetic reactant is m - moles*coef.

SPREAD SHEET INPUT OF SOLUTION DATA: I have used the following input file for the spreadsheet saturation-indices calculation:

TITLE Saturation-index calculation SOLUTION_SPREAD Number pH Ca Mg Na K Cl S(6) Alkalinity N(5) 1031080 18.0 2.655 1.094 1.222 0.026 1.300 0.183 2.947 1.451 ..... PRINT -saturation_indices true -selected_output true SELECTED_OUTPUT -file SI.res -saturation_indices Calcite Dolomite END

Will the calculation of SI for a temperature of 15 C be right if I add a column "temp" filled with 15?

That will work. It is also possible to insert a line within the SOLUTION_SPREAD data block.

-temp 15

and any solution compositions that follow will have the temperature of 15 C.

JUNK IN THE SELECTED-OUTPUT FILE (SPREAD-SHEET FILE): In selected results file I get a lot of waste columns. Can I set a required structure of this file?

Put "-reset false" at the top of the SELECTED_OUTPUT data block. It will turn off all default printing. You can then select the columns you want by setting the appropriate identifiers in SELECTED_OUTPUT to "true". In addition, you can use USER_PUNCH to select data columns to be printed to the selected output file. Basic functions in USER_PUNCH allow access to all data items controlled by identifiers in SELECTED_OUTPUT, thus it is possible to use USER_PUNCH exclusively to tailor your selected-output file.
SATO RELATION FOR DISSOLVED OXYGEN AND REDOX: We have a question about methods of determining pE in WATEQF. We'd like to know if there are significant differences between using the regular DO method and using the Sato DO method.

The Sato relation gives lower pe values. It's an ad hoc approach to try to produce pe values that are more consistent with Eh values that are measured with a platinum electrode in oxygenated waters. The relavence of the Sato pe to other redox processes is not clear.

The real question is what you are doing with the pe. Usually, the only thing the pe really affects is the speciation of iron and the saturation indices of ferric oxyhydroxide minerals. In general, the pe will affect the speciation of any redox element for which the TOTAL concentration of the element is entered, which may include Fe, U, As, and maybe a couple others. It doesn't affect speciation of redox elements that are entered in specific redox states: sulfur entered as sulfate, nitrogen entered as nitrate, etc. So if you don't have total concentrations of a redox element or you don't care about its speciation and saturation indices, don't worry about pe, Eh, Sato, or redox in general.

If you are worried about iron, then you're kind of out of luck anyway. The concentration is probably pretty low if it is an oxygenated water unless you have a pretty low pH. Between uncertain measurements of iron and the range in stability of Ferrihydrite, Goethite, hematite, and others, you can almost always speculate you are near equilibrium with some ferric oxyhydroxide.

I guess my bottom line is that 90% of the time it doesn't matter about the pe and the choice between the Sato and thermodynamic oxygen couple is simply a statement of some of the uncertainties.

One last option is to assume the water is initially in equilibrium with N2(g) at atmospheric pressure (0.8 atm). The N(5)/N(0) couple can then be used to calculate a pe.
AMOUNTS OF SOLIDS IN SIMULATIONS: How do I calculate the appropriate amount of solid (CEC, mineral amount) for PHREEQC calculations?

(1) Recalculate the concentration from 'per kg solid' (=s) to 'per liter pore water' (=q):
```
q = s * rho_b/eps_w,
       
```
where rho_b is bulk density of the sediment (kg/l), eps_w is water filled porosity (unitless).

(2) Recalculate to moles,

(3) enter moles of exchanger under EXCHANGE, moles of mineral under EQUILIBRIUM_PHASES, etc.

EXAMPLE: CEC = 1 meq/100 g; rho_b = 1.8 kg/l; eps_w = 0.3.
```
(1) CEC = 10 meq/kg * 1.8 / 0.3 = 60 meq / liter pore water
(2) X- = 0.06 mol/l
(3) EXCHANGE 1; X 0.06; -equil 1
       
```
AMOUNTS OF SOLIDS IN TRANSPORT SIMULATIONS: How are solid concentrations entered in transport simulations? Are they dependent on the discretization (no. of cells)?

Amounts of solid are entered for PHREEQC-2 in moles for each cell. This is not strictly a concentration, but usually the number of moles is nearly equal to 'moles per kg H2O' because generally, each cell contains a kilogram of water.

The amounts of solid are independent of discretization, the 'invisible' dimensions y and z make up for any changes in cell-length. It helps to visualize the column as extending in the y/z direction to make up for 1 kg H2O in each cell (flow is in the x direction).

When you change the discretization in TRANSPORT in PHREEQC-2, adapt:
```
-cells
-lengths
-time_step
-shifts
-print/punch options for cells
      
```
Nothing changes in the definitions of SOLUTION, EXCHANGE, SURFACE, etc. And, dispersivity is a medium property, which is not affected by the discretization.

EXAMPLE: 1 m sediment, 1 mM CaCl2 solution; CEC = 1 meq/100g; flush 4 times with 5 mM Mg(NO3)2; discretization in 5 and 10 cells, respectively.
```
SOLUTION 1-10; Ca 1; Cl 2; EXCHANGE 1-10; X 0.06; -equil 1
SOLUTION 0; Mg 1; N(5) 2
TRANSPORT; -cells  5; -length 0.2; -shifts 20; -time_step 1e6; -disp 0.1
#TRANSPORT; -cells 10; -length 0.1; -shifts 40; -time_step 2e6; -disp 0.1
END 
      
```
UNITS OF ALKALINITY, ESPECIALLY AS CaCO3: Why is it that when I specify alkalinity as HCO3 that I get the same number of moles of alkalinity as if I just give the concentration without the qualifier?

Also, I am confused about entering alkalinity. I have a reported alkalinity as mg/l of CaCO3. When I specify the concentration and use the qualifier "as CaCO3", the correct number of moles are calculated (assuming 50 for gfw). If I just enter the number, I get nothing that resembles the number of moles expected. I thought "as CaCO3" was the default?

The default gram formula weight for alkalinity Ca0.5(CO3)0.5 confuses me even more on this issue. I get the correct number of moles on the assumption that 50 is the gram formula weight, but the gram formula wt of CaCO3 is 100. Why do I get these results and why do I get a warning message if I specify "as CaCO3".

First, you should get the same number if the default units are in terms of moles (mol, mmol, or umol). The "as HCO3" would have no effect. Default units are mmol/kgw or whatever is specified with the "-units" identifier. Whenever default units are in terms of mol, mmol, or umol, Alkalinity is interpreted as equivalents, meq, or ueq.

Second, if the default units are mass (g, mg, ug) per (kgw, kgs, or L), then for alkalinty, they must be converted to equivalents and the "as HCO3" would be used to calculate a gram equivalent wt of 61; the entered concentration in, say, mg/L would be converted by dividing by 61 to obtain meq/L and further conversion would lead to eq/kgw.

Third, by default the gram equivalent wt for Alkalinity is defined in the 4th column of SOLUTION_MASTER_SPECIES in the database that you are using. If you are using minteq.dat, that entry is 61.0173 (if this field is a number, for alkalinity it is treated as the gram equivalent wt), which is the gram equivalent wt for HCO3. In this case including "as HCO3" in the input file does nothing because it is the same as the default value in the minteq.dat database.

In phreeqc.dat, the 4th column is Ca0.5(CO3)0.5 (if this field begins with a character, then it is treated as a chemical formula from which the gram equivalent wt is calculated), which is the correct formula for the gram equivalent wt for CaCO3; the gew is about 50. If phreeqc.dat is used then "as HCO3" has an effect; the mg/L are divided by 61 from the input file instead of 50 from the database file to calculate the number of equivalents of alkalinity.

Finally, if you have 100 mg/L of alkalinity reported as CaCO3, you have 2 meq/L of alkalinity. That is because there are 2 equivalents per mole of CaCO3. The problem is the gram formula wt (to convert to moles) is 100 mg/mmol, but the gram equivalent wt is 50 mg/meq. All other concentrations, except for Alkalinity, are converted to moles, only Alkalinity is converted to equivalents. For Alkalinity, the 4th column of SOLUTION_MASTER_SPECIES is interpreted as the gram equivalent wt or a formula by which a gram equivalent wt can be calculated; all other elements, the entry is for the gram formula weight.

email: dlpark@usgs.gov

[Home]
[Previous]

URL: https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/faq.html
Contact:dlpark@usgs.gov

Last modified: Feb 3, 1999
Visitor number 15527 since Sept 1, 1998