Davies equation in PHREEQC

[Wonwoo Yoon]

Dear all,

I'm trying to build a reactive transport model using the concept of a sequential non-iterative approach. In my model, PHREEQC is working as a reaction step.

Since my model is including mineral precipitation and dissolution, I'm trying to verify my model with the results from other reactive transport models including precipitation and dissolution.

However, I realized that the verification problem I choose utilizes a slightly different (modified, I guess) Davies equation to calculate the activity coefficient.

Therefore I'm trying to find out how to modify the equation setting for the calculation of activity coefficient in PHREEQC.

Thanks.

[dlparkhurst]

Hopefully, you are using PhreeqcRM. It is a module designed to link PHREEQC with transport models. It will save you a lot of time.

However you are using PHREEQC as your reaction module, the simplest testing would be to compare to 1D reactive-transport modeling with PHREEQC itself. You could easily use the same aqueous models simply by your choice of database.

Databases determine the activity coefficient model that is used; however, none of the databases use exclusively the Davies equation. The equations for activity coefficients are described in the Version 2 manual, which is distributed with all versions of PHREEQC. For databases other than sit.dat, pitzer.dat, and llnl-type databases, activity coefficient parameters are defined with the -gamma identifier in the SOLUTION_SPECIES data block. If -gamma is not defined for a species, then the Davies equation is used. Note that a parameter of 0.3 is hard-coded for the Davies equation. To make a Davies database, I would make a new database by removing all of the -gamma definitions from the phreeqc.dat database.

[Wonwoo Yoon]

Thanks for your reply, Dr. Parkhurst.

Yes, I'm aware of the -gamma parameters. To be specific, the paper I'm trying to benchmark is using the exactly same method as that of PHREEQC.

They specified the parameters to calculate the activity coefficients using the Debye-Huckel theory just as PHREEQC. So, it was possible to modify the -gamma parameters in PHREEQC database.

They also described that if the parameters are missing for a specific reaction, they used the Davies equation to calculate the activity coefficients.

However, the only thing I'm not sure about is whether I can change the parameter of 0.3 IN Davies equation to 0.24 since they used 0.24.

If there is no way to modify it, I would rather compare the results to 1D reactive-transport modeling with PHREEQC itself as you suggested.

Thanks a lot for your kind response!

[dlparkhurst]

You could change the source code, but I think you could get the activity coefficients that you need by careful definition of the -gamma parameters.

If you compare the Davies and extended Debye Huckel equation, they will be the same for a species if

a0 = 1/B, where B is given by the Basic function DH_B, and

bi = Az^2*0.3, where A is given by the Basic function DH_A.

The following shows a successful comparison of single and double charged species using Davies and -gamma.

For your case, you would use 0.24 instead of 0.3 in the calculation. Note that DH_A and DH_B are temperature dependent, so the activity coefficients would agree at only one selected temperature, or would require different definitions for each temperature. Hopefully, the example you are testing is constant temperature.

Code: [Select]

SOLUTION
-temp 25
USER_PRINT
10 PRINT "a0:  ", 1/DH_B
20 PRINT "bi+:  ", DH_A*0.3
30 PRINT "bi+2: ", DH_A*4*0.3
END
PRINT
-user_print false
SOLUTION_MASTER_SPECIES
    X_davies      X_davies+        0     1               1
    X_dh          X_dh+            0     1               1
    Xa_davies     Xa_davies+2      0     1               1
    Xa_dh         Xa_dh+2          0     1               1
SOLUTION_SPECIES
#a0:     3.0442e+00
#bi+:     1.5301e-01
#bi+2:    6.1203e-01 = 4*1.5301e-01
X_davies+ = X_davies+
X_dh+ = X_dh+
-gamma 3.0442e+00 1.5301e-01
Xa_davies+2 = Xa_davies+2
Xa_dh+2 = Xa_dh+2
-gamma 3.0442e+00 6.1203e-01
END
SOLUTION 2
-temp 25
X_davies 1
X_dh     1
Xa_davies 1
Xa_dh     1
REACTION
NaCl 1
2 in 20 steps
USER_GRAPH 1
    -headings               Mu X_davies+ X_dh+ Xa_davies+2 Xa_dh+2
    -axis_titles            "Ionic strength" "Activity coefficient" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X MU
20 GRAPH_Y GAMMA("X_davies+"), GAMMA("X_dh+")
30 GRAPH_Y GAMMA("Xa_davies+2"), GAMMA("Xa_dh+2")
  -end
    -active                 true
END


[Wonwoo Yoon]

Yes, the example is constant temperature, so it will be fine to apply the method you suggested.

I'll try that way, Thanks for your help!

[dlparkhurst]

I still hope you use PhreeqcRM. It simplifies tasks such as converting units between transport and chemistry, populating the solutions and reactions in the model cells, output, saving state, running chemistry on all of the model cells, parallelization, and many other tasks.