Question about the -equilibrate command

[Pak]

When defining a Surface data block in Phreeqc there is the option to use the command:
-equilibrate number

In the manual is written:
"-equilibrate—Indicates that the surface assemblage is defined to be in equilibrium with a given solution composition. Optionally, equil, equilibrate, -e[quilibrate], equilibrium, or -e[quilibrium].
number—Solution number with which the surface assemblage is to be in equilibrium. Any alphabetic characters following the identifier and preceding an integer (“with solution” in Line 1) are ignored."

So, for instance, if I have a surface with the cell number 4 I can equilibrate it with the solution of the cell 1. But there is something that I do not grasp.

If I have a  Solution 1 and a Surface 1 (So both are in the cell 1), and I run two simulations (e.g. Simplification of example 14 of Phreeqc Manual).
                                              First simulation: I use the command -equilibrate 
                                              Second simulation: I do not use the -equilibrate. I got different results.

In the first Phreeqc proceeds with a) Solution speciation, b) Surface-composition calculations, and c) Batch-reaction calculation.
The second Phreeqc proceeds with a) Solution speciation, and b) Batch-reaction.

I was expecting to get the same results in the Batch-reaction of the first and second simulation, but that is not the case. Actually, the results of the Surface-composition calculation are equal to the ones of the Batch-reaction for the first simulation.

What is the reason for these difference? Is the same algorithm applied ? Do the initial guess change? In one case, is something keep constant, that is not keep in the other case?

Thanks

[dlparkhurst]

In general, the two calculations are different because the initial condition of the surface is different.

When you use -equil, the surface composition is calculated to be in equilibrium with the solution specified. The surface definition now includes the moles of elements sorbed, plus possibly the diffuse layer elements.

When you do not use -equil the surface is initially SurfOH, with no additional elements.

So, when you go to the reaction calculation, say a solution and the surface, the moles of elements in the system are different. The first case has the solution plus the surface plus the sorbed elements; the second case has the solution plus the surface as SurfOH. So if the solution is the solution used to equilibrate the surface, the first case will generate no reaction because the surface is already in equilibrium with the solution (because of the -equil definition in SURFACE); in the second case, species that sorb are removed from the solution, so the solution composition changes and the surface composition changes from its original SurfOH state.

[Pak]

I am not sure if I have understood the answer properly. Briefly, I have understood that the calculation procedure is pretty much the same, but the difference is given by the initial conditions (So, initial guesses of the Newton-Rapshon algorithm?). Is that correct?

Then, in a practical case, should I expect better results using the -equilibrate command?


Thank you

[John Mahoney]

The difference is that with the -equilibrate option you have added mass to the system and that mass is on the surface.   This is similar to, but not identical to, fixing a pH in a solution (In that case you add or subtract protons from the solution).   Consider a beaker of solution with some concentration of a potentially sorbing element let's say lead (Pb+2).  If I add a surface, say pure gibbsite (AlOOH), to the beaker then the lead concentration in solution will decrease as is moves onto the gibbsite surface (we did not equilibrate the surface prior to the addition of the gibbsite, so there was no addition of lead to the system).  But now the solution and surface are equilibrated with each other.   

If we were to come back to this beaker some time later we would have something analogous to a ground water system.  Now we start fresh, imagine we just lost all of our notes about what happened earlier.  We could analyze the solution in the beaker and also the lead content on the gibbsite,  we could calculate the equilibrium relationship.   

Now consider a groundwater system (this is like our forgotten beaker).   We go out and collect a soil sample and a water sample; there is a measured concentration of dissolved lead in the solution and an amount of lead on the gibbsite in the soil.  They are at equilibrium and their concentrations are in some part related through some sorption reaction.   The composition of the solution over time is static it does not change.  There is a distribution coefficient (a Kd) relating the lead in the soil to the lead in the solution.  If we want to model this we have two options.   We can start with the solution concentrations and equilibrate the surface with the measured concentration of lead. The solution composition does not change but we add lead to the surface and we should get the same Kd.   OR we can measure the solution concentration, and the amount of gibbsite surface ( but in this case we do not measure the lead concentration with the gibbsite)  but we know the  amount of surface sites. so we set up a model and equilibrate the surface with the solution; more  lead is added  to the system ( it all goes onto the surface) and in a perfect model we get the same Kd. 

So it is a matter of the process you are modeling that dictates whether you equilibrate the surface or not.  Also remember that each component will have its own unique equilibrium between the solution and the surface so the Kd's will all be different.  This makes it harder to figure out but if you think about it it makes sense. 

[dlparkhurst]

No, it is not related to initial guesses of the Newton Raphson. As John said the best approach depends on the system you are modeling, but the conceptual model between the two approaches is different.

Consider surface 1, which has 1 mol of SurfOH and solution 1, which has 0.1 mol/kgw Pb. If you put a liter of solution with the surface and let them react, there is 0.1 mol of Pb that is distributed between the surface and the solution. The concentration of Pb in the reacted solution will be less than 0.1 mol/kgw.

Alternatively, equilibrate surface 1 with solution 1 with an initial surface calculation, and say the surface is calculated to have 0.05 mol of Pb. Now put the equilibrated surface together with solution 1 and let them react. The system now has 0.15 mol Pb that is distributed between the surface and the solution. However, the composition of the solution does not change and the composition of the surface does not change because you have defined the surface to be in equilibrium with the solution. So the solution still has 0.1 mol/kgw and the surface has 0.05 mol of Pb.

[Pak]

Lot of thanks for your explanations. I get it now !!!

Thank you.