SIT formulation

[pfdsantos]

Hello Mr Parkhurst,

I have been a phreeqc user for some time, exploring solubility calculations using the Pitzer and SIT formalisms.
In the case of SIT, phreeqc versions allow us to work with two ionic interaction parameters - epsilon and epsilon1, using a correction equation through ionic strength: epsilon + epsilon1*I.

I have been trying to change this equation to be consistent with what NEA-Group suggests: epsilon + epsilon1*log(I), which could further improve the performance of the code for ionic interactions other than 1:1 (2:1, 3:1, etc.).

I tried to change the equation directly in the sit.cpp file (I believe on line 413), but I was unsuccessful. So I imagine I am forgetting to change something else as well...

Do you have any suggestions or tips on what I should do?

Thanks a lot!

All the best,

Pedro

[dlparkhurst]

I need to see the equations that you want to solve to be confident, but I think you need the next code block after 413--case TYPE_SIT_EPSILON_MU. I think you need to change the I to log10(I). Note that it is a calculation for log10(gamma), not gamma. You also need to correct the equation for the osmotic coefficient that follows.

If you succeed in making the changes, I think I should add another parameter in addition to epsilon and epsilon1, something like epsilon_log_mu so that the program is backward compatible, but you can implement your calculations with the new parameter. I'll send a message so that we can continue the conversation offline.

Code: [Select]

case TYPE_SIT_EPSILON:
sit_LGAMMA[i0] += sit_M[i1] * param;
sit_LGAMMA[i1] += sit_M[i0] * param;
if (z0 == 0.0 && z1 == 0.0)
{
OSMOT += sit_M[i0] * sit_M[i1] * param / 2.0;
}
else
{
OSMOT += sit_M[i0] * sit_M[i1] * param;
}
break;
case TYPE_SIT_EPSILON_MU:
sit_LGAMMA[i0] += sit_M[i1] * I * param;
sit_LGAMMA[i1] += sit_M[i0] * I * param;
OSMOT += sit_M[i0] * sit_M[i1] * param;
if (z0 == 0.0 && z1 == 0.0)
{
OSMOT += sit_M[i0] * sit_M[i1] * param * I / 2.0;
}
else
{
OSMOT += sit_M[i0] * sit_M[i1] * param * I;
}
break;