PhreeqcUsers Discussion Forum
Processes => Reactive transport modelling => Topic started by: ahockin on October 21, 2021, 04:27:59 PM

Hi!
I am trying to understand how PHREEQC corrects the diffusion coefficient for temperature for each solute in multicomponent diffusion.
Using the formula from the help manual:
Dw' = Dw_298 × (T/298) ×(n_298/n_T),
where Dw' is the (temperature corrected) tracer diffusion coefficient; Dw_298 is the tracer diffusion coefficient at 25 C, T is the temperature in K and n is the viscosity.
I define the tracer diffusion coefficient using SOLUTION_SPECIES and tried to manually calculate the resulting diffusion coefficient corrected for temperature and viscosity. However I did not get the same result as output from PHREEQC. What is odd is that if I use the new viscosity, but not the temperature [Dw' = Dw_298 ×(n_298/n_T)], I get the same value as PHREEQC gives. See attached code.
Any help on why this is would be really helpful, thanks in advance!

First, part of the problem may be that the documentation for the function DIFF_C is incorrect. The function returns the diffusion coefficient at the current temperature, not at 25 C.
The diffusion coefficient at temperature for a species depends on two definitions in SOLUTION_SPECIESDw and Dw_t. It also depends on the viscosity of pure water at 25 C and the viscosity of water at the working temperature. (Unfortunately, no correction of the viscosity for solutes is made at this time.) Note only phreeqc.dat, Amm.dat, and pitzer.dat have the necessary SOLUTION_SPECIES definitions.
The attached script calculates the diffusion coefficient of Na+ as a function of temperature with USER_PRINT. The key equation is
dw_temp = dw * EXP(dw_t/TK  dw_t/298.15) * viscos_0_25 / VISCOS
SOLUTION
temp 25
Na 1
Cl 1
REACTION_TEMPERATURE
10 20 30 40 50
USER_PRINT
10 REM Calculation of temperature dependence of Na+ diffusion coefficient
20 dw = 1.33e9 # SOLUTION_SPECIES
30 dw_t = 122 # SOLUTION_SPECIES
40 viscos_0_25 = 0.8900239182946
50 dw_temp = dw * EXP(dw_t/TK  dw_t/298.15) * viscos_0_25 / VISCOS
60 PRINT "DIFF_C(Na+): ", DIFF_C("Na+")
70 PRINT "Dw(TK): ", dw_temp
END

Hi Mr. Parkhurst,
Thanks for replying so quickly.
One clarification, in the documentation for SOLUTION_SPECIES the formula for calculating the diffusion coefficient is given as:
Dw(TK) = D * exp(damp / TK  damp / 298.15) * TK * 0.89 / (298.15 * viscos)
But the formula you gave above does not include the bolded items, namely the working temperature (TK) and temperature of 25 deg (in Kelvin, 298.15).
I'm wondering why these are not included in your formula? This is where my initial confusion arose in trying to manually calculate the change in diffusion coefficient with temperature and match it to what PHREEQC outputs.
Thanks again!

I'll have to ask Tony. I was copying from the code, but I'll check more carefully.

The function DIFF_C is returning the wrong value. The correct value would be multiplied by the bolded ratio. Also, the return value of SETDIFF_C needs to be multiplied by the bolded ratio.
The values used in TRANSPORT calculations are consistent with the formula in SOLUTION_SPECIES. The function DIFF_C and the return value from SETDIFF_C are not used directly in any calculations; they are only for user output.
I'll check with Tony, but make sure the next release fixes any bugs.

Super, thank you for letting me know. That clears up my confusion!