Documentation: Solution.I and Solution.mu (Ionic Strength) unit not mol/L (likely mol/kgw)

[j-bl]

Issue Description

The documentation of the Solution.I property and the Solution.mu method both state that the ionic strength is returned in $\mathrm{mol/L}$. However, the returned unit seems to differ and is likely molality ($\mathrm{mol/kgw}$).

Steps to Reproduce

Here is a minimal reproducible example using a 4 $\mathrm{mol/kgw}$ NaCl solution:

from phreeqpython import PhreeqPython

pp = PhreeqPython(
    database="pitzer.dat",
)

# Create a concentrated 4 mol/kgw NaCl solution.
sol = pp.add_solution({
    'units': 'mol/kgw',
    'Na': 4.0,
    'Cl': 4.0
})

# 1. Get total element moles.
moles_na = sol.species_moles['Na+']
moles_cl = sol.species_moles['Cl-']

# 2. Calculate ionic strength manually using water mass (kg).
# We use:
#   I = 0.5 * sum(m_i * z_i^2),
# where m_i is the molality of ion i and z_i is the charge of ion i.
# The resulting unit is mol/kgw.
mass_water_kg = sol.mass # Mass of water
m_na = moles_na / mass_water_kg
m_cl = moles_cl / mass_water_kg
i_molal_calculated = 0.5 * (m_na * (1**2) + m_cl * ((-1)**2))

# 3. Calculate ionic strength manually using solution volume (liters).
# We use:
#   I = 0.5 * sum(c_i * z_i^2),
# where c_i is the molarity of ion i and z_i is the charge of ion i.
# The resulting unit is mol/L.
vol_liters = sol.volume # Volume of the total solution in liters
c_na = moles_na / vol_liters
c_cl = moles_cl / vol_liters
i_molar_calculated = 0.5 * (c_na * (1**2) + c_cl * ((-1)**2))

print(f"Reported sol.I:\t\t\t{sol.I:.8f}")
print(f"Reported sol.mu():\t\t{sol.mu():.8f}")
print(f"Calculated Molal (mol/kgw):\t{i_molal_calculated:.8f}")
print(f"Calculated Molar (mol/L):\t{i_molar_calculated:.8f}")

Output

Reported sol.I:			4.00000013
Reported sol.mu():		4.00000013
Calculated Molal (mol/kgw):	4.00000000
Calculated Molar (mol/L):	3.68640569

PhreeqPython's output matches the ionic strength as calculated in molality, but not that calculated in molarity (which the documentation suggests). This suggests that the documented unit of $\mathrm{mol/L}$ may be incorrect.

Suggested Fix

I believe this is only a documentation error caused by a misinterpretation of the unit returned by the VIPhreeqc backend. Thus, only the documentation for solution.I and solution.mu should be updated to state the correct unit to prevent confusion (particularly when working with concentrated solutions).

Environment

PhreeqPython 1.6.2 with Python 3.10 on Windows.