NaCl solubility

[Konjkov]

Hello pyEQL developers, I wrote a code that calculates the solubility of NaCl at different temperatures:

#!/usr/bin/env python3

from math import exp, sqrt, nan
import periodictable
from pyEQL import Solution
from scipy.optimize import root

def K_sp(temperature):
    R = 8.314
    # https://webbook.nist.gov/cgi/cbook.cgi?ID=C7647145&Mask=2#Thermo-Condensed
    dH0_NaCl = -411_120      # Дж/моль
    S0_NaCl = 72.13          # Дж/моль/К
    # https://www.mrbigler.com/misc/energy-of-formation.PDF
    dH0_Na_plus = -240_340   # Дж/моль
    dS0_Na_plus = 59.0       # Дж/(моль*К)

    dH0_Cl_minus = -167_080  # Дж/моль
    dS0_Cl_minus = 56.5      # Дж/(моль*К)

    dG0_solubility = dH0_Na_plus + dH0_Cl_minus - dH0_NaCl - (dS0_Na_plus + dS0_Cl_minus - S0_NaCl) * (temperature + 273.15)
    return exp(-dG0_solubility/(R * (temperature + 273.15)))

def fun(concentration, temperature):
    sol = Solution({'Na+': f'{concentration} mol/l', 'Cl-': f'{concentration} mol/l'}, temperature=f'{temperature} degC', engine='native')
    return (
        sol.get_activity_coefficient('Na+')
        * sol.get_amount('Na+', 'mol/L')
        * sol.get_activity_coefficient('Cl-')
        * sol.get_amount('Cl-', 'mol/L')
    ).magnitude - K_sp(temperature)

# Solubility NaCl g/100g H2O
data = (
    (-21.2,  nan),
    (-10,  nan),
    (  0,  35.65),
    ( 10,  35.72),
    ( 20,  35.89),
    ( 30,  36.09),
    ( 40,  36.37),
    ( 50,  36.69),
    ( 60,  37.04),
    ( 70,  37.46),
    ( 80,  37.93),
    ( 90,  38.47),
    (99.9,  38.99),
    (100,  38.99),
    (108.7, nan)
)

for temp, solubility in data:

    concentration = root(fun=fun, x0=sqrt(K_sp(temp)), args=(temp))['x'][0]
    sol = Solution({'Na+': f'{concentration} mol/l', 'Cl-': f'{concentration} mol/l'}, temperature=f'{temp} degC', engine='native')
    amount_NaCl = sol.get_amount('Na+', 'g/L') + sol.get_amount('Cl-', 'g/L')
    amount_H2O = sol.get_amount('H2O', 'kg/L')

    print(f'Solubility {amount_NaCl / amount_H2O:.2f}, {solubility*10:.2f} at {temp} degC')

Unfortunately, it does not work at temperatures of 100C and above, although it is known that a saturated solution of NaCl boils at 108.7С.

the following RuntimeWarning is returned instead
/home/vladimir/.local/lib/python3.10/site-packages/pint/facets/plain/quantity.py:1271: RuntimeWarning: invalid value encountered in scalar power magnitude = new_self._magnitude**exponent

Best wishes, Vladimir.

[rkingsbury]

Thank you for reporting @Konjkov . Can you post the full stacktrace (the complete error message, not just the last line?). I have an idea of where this might be coming from but it's hard to be sure without more context.

[Konjkov]

this is the full traceback, but I noticed that the activity coefficients sol.get_activity_coefficient('Na+') and sol.get_activity_coefficient('Cl-') is NaN at 100 C and further calculations become meaningless

[rkingsbury]

Hmm OK. I suspect that the problem originates in the model we use for water properties (IAPWS97), which may not be defined above 100C.

In this line, we create an IAPWS97 instance that is subsequently used to retrieve things like density, viscosity, and dielectric constant of pure water. These quantities enter into activity coefficient calculations.

So the first troubleshooting step I would suggest is to manually create an IAPWS97 instance and see if you can retrieve those properties (rho, mu, and dielectric_constant) at temperatures above 100 C.