Specific Conductance calculation

[Isil]

Hi,
I am trying to reproduce the calculation of Specific Conductance that PHREEQC performs.
On PHREEQC, I have calculated the speciation of different solutions: from a solution with a single salt dissolved to a solution that reproduces seawater. Then, on an excel file, I implemented the calculation of specific conductance with the formula present at the link http://www.hydrochemistry.eu/exmpls/sc.html (from PHREEQC manual version 3).
A and B parameters are calculated by the formula present in the following paper: https://core.ac.uk/download/pdf/84005019.pdf at page 38.

I use the molalities provided by the PHREEQC speciation and the diffusion coefficient (Dw(TK), d, a1, a2) in the phreeqc database.

My problem is that the specific conductance calculated is higher than that calculated by PHREEQC of 15-20%, depending on the solution composition.

I have attached a figure with the formula used.

Does anyone know where I am wrong? In the calculation of parameters A and B?
Thank you

[dlparkhurst]

Here is a script that calculates the specific conductance of a 1 molal CaCl2 solution. There are minor contributions from other ions in solution (H+, OH-, CaOH+), but they are negligible. You'll have to figure out the difference.

Code: [Select]

SOLUTION
-units mol/kgw
Ca 1
Cl 2
USER_PRINT
-start
10 DATA "Ca+2", 0.793e-9,  97,  3.4,  24.6, 2, \
        "Cl-", 2.03e-9,  194,  1.6,  6.9, -1
20 RESTORE 10
30 FOR i = 1 TO 2
40   READ spec$, dw, dw_t, dw_a, dw_a2, z
50   l_z = ABS(z)
60   dw = dw * exp(dw_t/TK - dw_t/298.15)
70   ka = DH_B * dw_a2 * MU^0.5 / (1 + MU^0.75)
80   ff = exp(-dw_a * DH_A * l_z * MU^0.5 / (1 + ka))
90   dw = dw * ff
100  dw_t_SC = MOL(spec$) * l_z * l_z * dw
110  cond = cond + dw_t_SC
120 NEXT i
130 F_C_MOL = 96493.5
140 R_KJ_DEG_MOL = 0.00831470
150 cond = cond * 1e7 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298150)
200 PRINT cond, SC
210 END
-end
END


[dlparkhurst]

For a temperature other than 25 C, I left out the viscosity correction. You need to add the following line to the previous script, where the numerator is the viscosity of pure water at 25 C.

Code: [Select]

160 cond = cond * 0.89002 / VISCOS


[Isil]

Thank you very much, now it works. The main problem was in parameters A and B.
I corrected the calculation of B, but I still don't understand how the parameter A is calculated. If I set A equal to 0.51002 (mol/kg)^-0.5, the calculation of specific conductance is correct, but I would like to understand the origin of this value.

[dlparkhurst]

Here are the methods where DH_A, DH_B, and rho_0 are calculated.

Code: [Select]

LDBLE Phreeqc::
calc_dielectrics(LDBLE tc, LDBLE pa)
/* ---------------------------------------------------------------------- */
{
/* Relative dielectric constant of pure water, eps as a function of (P, T)
   Bradley and Pitzer, 1979, JPC 83, 1599.
   (newer data in Fernandez et al., 1995, JPCRD 24, 33,
  and Fernandez et al., 1997, JPCRD 26, 1125, show its correctness)
   + d(eps)/d(P), Debye-Hueckel A and B, and Av (for Av, see Pitzer et al., 1984, JPCRD 13, p. 4)
*/
if (llnl_temp.size() > 0) return OK;
if (tc > 350.)
{
tc = 350.;
}
LDBLE T = tc + 273.15;
LDBLE u1 = 3.4279e2, u2 = -5.0866e-3, u3 = 9.469e-7, u4 = -2.0525,
u5 = 3.1159e3, u6 = -1.8289e2, u7 = -8.0325e3, u8 = 4.2142e6,
u9 = 2.1417;
LDBLE d1000 = u1 * exp(T * (u2 + T * u3)); // relative dielectric constant at 1000 bar
LDBLE c = u4 + u5 / (u6 + T);
LDBLE b = u7 + u8 / T + u9 * T;
LDBLE pb = pa * 1.01325; // pa in bar
eps_r = d1000 + c * log((b + pb) / (b + 1e3)); // relative dielectric constant
if (eps_r <= 0)
{
eps_r = 10.;
warning_msg("Relative dielectric constant is negative.\nTemperature is out of range of parameterization.");
}

/* qe^2 / (eps_r * kB * T) = 4.803204e-10**2 / 1.38065e-16 / (eps_r * T)
   = 1.671008e-3 (esu^2 / (erg/K)) / (eps_r * T) */
LDBLE e2_DkT = 1.671008e-3 / (eps_r * T);

DH_B = sqrt(8 * pi * AVOGADRO * e2_DkT * rho_0 / 1e3);  // Debye length parameter, 1/cm(mol/kg)^-0.5

DH_A = DH_B * e2_DkT / (2. * LOG_10); //(mol/kg)^-0.5

/* A0 in pitzer */
if (pitzer_model || sit_model)
{
A0 = DH_B * e2_DkT / 6.0;
if (pitzer_model && aphi != NULL)
{
calc_pitz_param(aphi, T, 298.15);
A0 = aphi->p;
}
}

/* Debye-Hueckel limiting slope = DH_B *  e2_DkT * RT * (d(ln(eps_r)) / d(P) - compressibility) */
DH_Av = DH_B * e2_DkT * R_LITER_ATM * 1e3 * T * (c / (b + pb) * 1.01325 / eps_r - kappa_0 / 3.); // (cm3/mol)(mol/kg)^-0.5

DH_B /= 1e8; // kappa, 1/Angstrom(mol/kg)^-0.5

/* the Born functions, * 41.84 to give molal volumes in cm3/mol... */
ZBrn = (-1 / eps_r + 1.0) * 41.84004;
QBrn = c / (b + pb) / eps_r / eps_r * 41.84004;
/* dgdP from subroutine gShok2 in supcrt92, g is neglected here (at tc < 300)...
   and, dgdP is small. Better, adapt Wref to experimental Vm's */
dgdP = 0;

return (OK);
}


Code: [Select]

LDBLE Phreeqc::
calc_rho_0(LDBLE tc, LDBLE pa)
/* ---------------------------------------------------------------------- */
{
/* Density of pure water
        Wagner and Pruss, 2002, JPCRD 31, 387, eqn. 2.6, along the saturation pressure line +
interpolation 0 - 300 oC, 0.006 - 1000 atm...
    */
if (llnl_temp.size() > 0) return OK;
if (tc > 350.)
{
if (need_temp_msg < 1)
{
std::ostringstream w_msg;
w_msg << "Fitting range for dielectric constant of pure water is 0-350 C.\n";
w_msg << "Fitting range for density along the saturation pressure line is 0-374 C,\n";
w_msg << "                         for higher pressures up to 1000 atm    0-300 C.\n";
w_msg << "Using temperature of 350 C for dielectric and density calculation.";
warning_msg(w_msg.str().c_str());
need_temp_msg++;
}
tc = 350.;
}
LDBLE T = tc + 273.15;
//eqn. 2.6...
LDBLE Tc = 647.096, th = 1 - T / Tc;
LDBLE b1 = 1.99274064, b2 = 1.09965342, b3 = -0.510839303,
b4 = -1.75493479, b5 = -45.5170352, b6 = -6.7469445e5;
    rho_0_sat = 322.0 * (1.0 + b1 * pow(th, (LDBLE) 1./3.) + b2 * pow(th, (LDBLE) 2./3.) + b3 * pow(th, (LDBLE) 5./3.) +\
               b4 * pow(th, (LDBLE) 16./3.) + b5 * pow(th, (LDBLE) 43./3.) + b6 * pow(th, (LDBLE) 110./3));
//pressure...
LDBLE p0 =  5.1880000E-02 + tc * (-4.1885519E-04 + tc * ( 6.6780748E-06 + tc * (-3.6648699E-08 + tc *  8.3501912E-11)));
LDBLE p1 = -6.0251348E-06 + tc * ( 3.6696407E-07 + tc * (-9.2056269E-09 + tc * ( 6.7024182E-11 + tc * -1.5947241E-13)));
LDBLE p2 = -2.2983596E-09 + tc * (-4.0133819E-10 + tc * ( 1.2619821E-11 + tc * (-9.8952363E-14 + tc *  2.3363281E-16)));
LDBLE p3 =  7.0517647E-11 + tc * ( 6.8566831E-12 + tc * (-2.2829750E-13 + tc * ( 1.8113313E-15 + tc * -4.2475324E-18)));
/* The minimal pressure equals the saturation pressure... */
if (ah2o_x <= 1.0)
p_sat = exp(11.6702 - 3816.44 / (T - 46.13)) * ah2o_x;
else
p_sat = exp(11.6702 - 3816.44 / (T - 46.13));
//ah2o_x0 = ah2o_x; // for updating rho in model(): compare with new ah2o_x
if (pa < p_sat || (use.Get_solution_ptr() && use.Get_solution_ptr()->Get_patm() < p_sat))
{
pa = p_sat;
}
if (!use.Get_gas_phase_in())
patm_x = pa;
pa -= (p_sat - 1e-6);
rho_0 = rho_0_sat + pa * (p0 + pa * (p1 + pa * (p2 + sqrt(pa) * p3)));
if (rho_0 < 0.01)
rho_0 = 0.01;

/* compressibility, d(ln(rho)) / d(P), 1/atm... */
kappa_0 = (p0 + pa * (2 * p1 + pa * (3 * p2 + sqrt(pa) * 3.5 * p3))) / rho_0;

return (rho_0 / 1e3);
}


[Isil]

Thank you, it is very useful.

[dlparkhurst]

Here is a graph for anyone following the thread.