Vitens / phreeqpython

Object-oriented python wrapper for the VIPhreeqc module
Apache License 2.0
68 stars 19 forks source link

Conductivity discrepancy at high molality #20

Open trinberg opened 2 years ago

trinberg commented 2 years ago

Hi there,

First of all, thank you so much for this fantastic tool! It's really helping us in our research.

Second, I wanted to see if you had a solution to this issue I'm seeing:

I would like to find the conductivity of high molality solution (1 - 3 molal) for potassium carbonate solution. When I use PHREEQC Interactive (the USGS software) directly I get results that match chemistry handbooks. So, for example, entering this input and running the Pitzer model:

SOLUTION 0
units mmol/kgw
pH 10 charge
temp 20                 
K 1000
C 500 
END

Gives 73.3 mS/cm and pH = 11.79

BUT, when I use your phreeqpy function

pp = PhreeqPython(database='pitzer.dat')
solution = pp.add_solution({'units':'mol/kgw', #set the units (moles per kg of water)
                                 'pH': '10 charge',
                                 'temp': 20, 
                                'K':2,
                                'C': 1
                           })
pH = solution.pH
cond = solution.sc
print('pH:',pH)    
print('Cond:',cond)

I get 50.5 mS/cm and the same pH as above (11.79)

There is a substantial difference in conductivity and I'd really like to be able to reconcile that. The phreeqpy function is outputting a lower conductivity than physical measurements and than PHREEQC. Do you have any suggestions for why this may be happening? Note that at molality below ~0.5, I don't see the discrepancy anymore.

Thanks so much, and again, really appreciate this tool!

AbelHeinsbroek commented 2 years ago

Hello!

Thanks for your kind words, glad to hear that you like the work we've done!

In your phreeqpython example you use a different input then in your example given above, adding 2000/1000 instead of 1000/500 mol K and C, which obviously results in a different pH and conductivity!

However, I think the main discrepancy is caused by a difference in database file. The pitzer.dat file included with PhreeqPython hasn't been updated since the project inception (we never use it ourselves..) and therefore has diverted from the one maintained by the USGS. If you direct phreeqpython to the same database as PHREEQC interactive uses you should get the same results:

import phreeqpython

pp = phreeqpython.PhreeqPython(database='/Users/Abel/pitzer.dat')
pp.ip.debug = True
solution = pp.add_solution({'units':'mmol/kgw', #set the units (moles per kg of water)
                                 'pH': '10 charge',
                                 'temp': 20, 
                                'K':1000,
                                'C': 500
                           })
pH = solution.pH
cond = solution.sc
print('pH:',pH)    
print('Cond:',cond)

resulting in:

SOLUTION 0
  units mmol/kgw
  pH 10 charge
  temp 20
  K 1000
  C 500
SAVE SOLUTION 0
END 

pH: 11.789155189291549
Cond: 73290.5280343256
trinberg commented 2 years ago

Thanks for your help here! And my mistake, I did try the correct values (1000/500), but accidentally copied the wrong values into the example above! Sorry for the confusion!

I've redirected the database to point to the USGS repository, as you suggested, but I'm still somehow seeing the discrepancy. pp = PhreeqPython(database='C:/Program Files/USGS/phreeqc-3.6.2-15100-x64/database/pitzer.dat')

Just figured I'd respond here while I continue debugging in case you have any suggestions. Thanks again!

trinberg commented 2 years ago

Ok, I just wanted to follow up and say that I'm still thoroughly stuck on this. I'm able to call different databases (SIT, minteq, etc.) from the USGS folder to confirm that the database loading function is working properly.

BUT, the following function in python which ostensibly should be using pitzer:

import phreeqpython

pp = phreeqpython.PhreeqPython(database='C:/Program Files/USGS/phreeqc-3.6.2-15100-x64/database/pitzer.dat') #initializes a database instance
pp.ip.debug = True
solution = pp.add_solution({'units':'mmol/kgw', #set the units (moles per kg of water)
                                 'pH': '10 charge',
                                 'temp': 20, 
                                'K':1000,
                                'C': 500
                           })
out = solution.species #output variable
pH = solution.pH
cond = solution.sc
print('pH:',pH)    
print('Cond:',cond)

results in:

SOLUTION 0
  units mmol/kgw
  pH 10 charge
  temp 20
  K 1000
  C 500
SAVE SOLUTION 0
END 

pH: 11.789155189291549
Cond: 50457.00068057522

Note that the pH is exactly the same as yours, but the conductivity is still quite off. Do you have any further suggestions here?

If you run phreeqc.dat, what conductivity do you get for the same inputs? (I get: pH: 11.70069773360222 Cond: 55380.14209795132). Maybe that will help us narrow in on whether it's something with the databases or the functions that process them?

AbelHeinsbroek commented 2 years ago

I see that you use an older version of Phreeqc (3.6.2). What happens if you try the example with the most recent pitzer.dat? Would you mind sending me your pitzer.dat and phreeqc.dat files so I can check if I get the same results?

trinberg commented 2 years ago

databases to try.zip

Let me know what the outcome is for you! And thank you again for your help with this - it would be a great resource if we can get this to work properly.

AbelHeinsbroek commented 2 years ago

I get the proper results with your databases. I've updated the database files supplied with phreeqpython to the latest versions from USGS and added a case specific test:

https://github.com/Vitens/phreeqpython/blob/00a4248e1a4b57356ae3bec640e0819da52eccab/tests/test_phreeqpython.py#L204

The test seems to run ok on all platforms, including windows:

[00:00:26] %PYTHON%/python.exe -m nose -v
[00:00:27] tests.test_phreeqpython.TestPhreeqPython.test10_use_non_default_database_directory ... ok
[00:00:27] tests.test_phreeqpython.TestPhreeqPython.test11_pitzer ... ok
[00:00:27] tests.test_phreeqpython.TestPhreeqPython.test1_basiscs ... ok
[00:00:27] tests.test_phreeqpython.TestPhreeqPython.test2_solution_functions ... ok
[00:00:27] tests.test_phreeqpython.TestPhreeqPython.test3_mixing ... ok
[00:00:27] tests.test_phreeqpython.TestPhreeqPython.test4_solution_listing ... ok
[00:00:27] tests.test_phreeqpython.TestPhreeqPython.test5_addition ... ok
[00:00:27] tests.test_phreeqpython.TestPhreeqPython.test6_misc ... ok
[00:00:27] tests.test_phreeqpython.TestPhreeqPython.test7_dump_and_load ... ok
[00:00:27] tests.test_phreeqpython.TestPhreeqPython.test8_raw_solutions ... ok
[00:00:27] tests.test_phreeqpython.TestPhreeqPython.test9_gas_phases ... ok
[00:00:27] tests.test_utility.TestUtility.test_00_convert_units ... ok

Would you mind installing the database_update branch to test whether that one gives you the correct results?

trinberg commented 2 years ago

Fantastic! I just ran this branch and got the exact same values as you! pH: 11.789155189291549 Cond: 73290.5280343256

I'm not sure what was going on there, but this seems to have addressed it! Thank you! Do you think you'll update the master?

AbelHeinsbroek commented 2 years ago

I think the problem stems from PhreeqPython originally using an old version of IPhreeqc (3.3.x), which may have had some Pitzer related bugs. The branch you tested uses the latest version of IPhreeqc (3.7.3) and Pitzer.dat

At my company we only use PhreeqPython for calculating drinking water at low conductivity so have never noticed this problem.

I need to run a few more integration tests, but given that the new branch passes all unit tests on all platforms I am fairly sure I can safely replace the master and PyPy package soon.

trinberg commented 2 years ago

Sounds good - thanks for working through this. I'm working on the branch now, but will update to master once you've replaced it.

Much appreciated!

trinberg commented 1 year ago

Hey there - just following up on this thread and wondering if you all pushed this update to the master? My colleague just pip (re)installed PhreeqPython and seems to be getting different outputs for the same underlying code as me.