I have the following in my phreeqc input file

PHASES
BaSO4 = 1.0000 SO4-2 + 1.0000 Ba+2
	-analytic	3.630937e+003	1.218154e+000	-1.393850e+005	-1.437233e+003	5.362526e+006	-4.566328e-004
	-Vm	 52.1000

is it possible to do this with phreeqpython, either as a python statement or as a startupfile?

Great tool BTW!

I want to calculate the DIC in an open carbonate system, assuming just NaHCO3 or NaOH is present. The input should be the pH, T, and pCO2.

So far I just modelled it step by step.

#shift the pH with NaOH (for alkaline range) or HCl (for acidic range)
 sol.change_ph(pH,'NaOH')
 
#concentration CO2
p=415 # ppm

#partial pressure pCO2 in atm
pCO2=p*1e-6            
            
#phreeqc always uses log10 values (NaOH will be used to increase and HCl used to lower the pH)
input_pCO2=np.log10(pCO2)
            
            
#equalize solution with CO2
#equilibrium with the atmosphere
sol.equalize(['CO2(g)'], [input_pCO2])

print(sol.pH)

The problem is that I basically want to change pH in equilibrium with CO2(g) directly. With the current approach, the resulting pH will always be much lower which makes sense.

Is there a function that I can use to get directly the water speciation for a given pH, pCO2 and T ? Assuming just pure carbonate system alkalinity is present?

Hi, professor,
when I learn this case (Kinetic Dissolution of Quartz), I don't know how to get the value of quartz_dissolved.
When I add the command "print(quartz_dissolved)" in the "ratefun()" function, I can get a lot of value (just one time_step). However, this value is not the same with the value I got from the customer software.
I hope to get your help.
best wish!!

%pylab inline
import phreeqpython
from scipy.integrate import odeint
pp = phreeqpython.PhreeqPython('phreeqc.dat')
def ratefun(sol, quartz_dissolved, m0, A0, V):
m = m0-quartz_dissolved
rate = (A0/V)(m/m0)**0.6710**-13.7*(1-sol.sr("Quartz"))
print(quartz_dissolved) # I try this method, but it is wrong.
return rate * 1e3
solution1 = pp.add_solution({})

t = np.array([])
y = []

year = 365243600

for time, sol in solution1.kinetics('SiO2',
rate_function=ratefun,
time=np.linspace(0,5*year, 2),
m0=158.5,
args=(23.13,0.16)):
t = np.append(t, time)
y.append(sol.total_element('Si', units='mmol'))

Hello,
I am modeling the sorption of dissolved metals to ferrihydrite surfaces. Is there a way to use this code to access the SURFACE functionality of PHREEQC? Is there a method to create a surface or a work-around function? Thank you for time!

In the 'Features' section of the README.md file the wrong 'PhreeqPython' method is called:

...
# add a solution consisting of 1 mmol CaCl2 and 2 mmol NaHCO3
solution = pp.add_solution({'CaCl2':1.0,'NaHCO3':2.0})
...

should be:

...
# add a solution consisting of 1 mmol CaCl2 and 2 mmol NaHCO3
solution = pp.add_solution_simple({'CaCl2':1.0,'NaHCO3':2.0})
...

and

solution2 = pp.add_solution({'KCl':1.0})
# create mixture of 50% solution and 50% solution2
...

should be:

solution2 = pp.add_solution_simple({'KCl':1.0})
# create mixture of 50% solution and 50% solution2
...

Curious if you could provide some insight on how you can retrieve the amount of moles of a phase dissolved or precipitated. I can see how to do this for individual elements or species, but it is unclear how to do this for a phase.

for example, given the code, where I have water with a pH of 1 and I add Cu metal to the solution. If I equalize to a SI =0 for all phases, how can I see how much of each phase exists in the solution?

from phreeqpython import PhreeqPython
from pathlib import Path


db_dir = Path('C:/Users/Craig/source/repos/eWasteSim/Database')
db = Path('minteq.dat')

pp = PhreeqPython(database= db, database_directory= db_dir, debug= False)

# this is for simulation purposes, in the model this will be a copy of the block from the left
solidsInlet = {
            'CuMetal'           : 10.0, 
            'Atacamite'         : 0.0, 
            'Cu(OH)2'           : 0.0, 
            'Cuprite'           : 0.0, 
            'Melanothallite'    : 0.0, 
            'Nantokite'         : 0.0, 
            'Tenorite'          : 0.0
}

# This is for simulations purposes, in the model this will be a copy of the block from above
liquidsInlet = solution = pp.add_solution({
                            'units' :   'mol/kgw',
                            'pH'    :   1,
                            'temp'  :   50.0,
                            'redox' :   'pe',
                            'Cu'    :   0,
                            'Cl'    :   '1 charge',
                            'O(0)'  :   0
})

liquidsOulet = liquidsInlet.copy()
liquidsOulet.change_temperature(60)

phasesInlet = [k for k in solidsInlet.keys()]
inletMoles = [v for v in solidsInlet.values()]
SI_Target = [0]*len(phasesInlet)

liquidsOulet.equalize(phases= phasesInlet, to_si= SI_Target, in_phase= inletMoles)

This can be more complicated when more than one metal exists and the pH is changing. Can you recommend a method to obtain the moles of a given phase exist?

Thanks

Do you have an additional example of using surface complexation modeling with your current package?

I am trying to model pyrite dissolution and set up a rate function that does increase the total Fe, however, the S-concentration does not change.
Using the odeint approach, it is evident that I only return the rate for Fe, but I am unclear on how to map a second return value so that it affects the S-concentration. I use the following code

import phreeqpython
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import pprint


def rate_pyrite(pyrite_dissolved, time, sol, A0, m0):
    # acid solution parameters
    a1 = 2.82e02  # A
    e1 = 56900
    n1 = -0.500
    n3 = 0.500
    # neutral solution parameters
    a2 = 2.64e05
    e2 = 56900
    n2 = 0.5
    r = 8.31451

    temp = sol.copy()  # make a temporary copy of the solution
    temp.add("FeS2", pyrite_dissolved[0], "mol")  # add Pyrite to solution
    m = m0 - pyrite_dissolved[0]

    # get solution properties
    pyrite_sr = temp.sr("pyrite")
    hplus_ac = temp.activity("H+", units="mol")
    o2_ac = temp.activity("O2", units="mol")
    fe3plus_ac = temp.activity("Fe+3",  units="mol")
    tk = temp.temperature + 273.15  

    if m0 <= 0:
        sa = A0
    else:
        sa = A0 * ((m / m0) ** 0.67)

    if m < 0:
        rate = 0
    elif m == 0 and pyrite_sr < 1:
        rate = 0
    else:
        acid_rate = a1 * math.e ** (-e1 / (r * tk)) * hplus_ac**n1 * fe3plus_ac**n3
        neutral_rate = a2 * math.e ** (-e2 / (r * tk)) * o2_ac * n2  #
        rate = (acid_rate + neutral_rate) * (1 - pyrite_sr) * sa

    temp.forget()  # cleanup the no longer needed temporary solution
    return rate


# ----------------- code starts below  ------------------------#
pp = phreeqpython.PhreeqPython()

solution = pp.add_solution(
    {
        "units": "mol/kgw",
        "temp": 30,
        "O(0)": "1 O2(g) -0.7",
        "pH": 1,
        "Fe": 0.1,
        "S" : 0,
    }
)
solution.equalize(
    ["O2(g)", "CO2(g)"],
    to_si=[-0.7, -3.5],
)

hours = 60 * 60
time = np.linspace(0, 30 * hours, 100)
surface_area = 50
mass = 0.05  # mol

# solve differential equation
yy = odeint(
    rate_pyrite,
    0,
    time,
    args=(
        solution,
        surface_area,
        mass,
    ),
)

pprint.pp(solution.species)
plt.plot(time / hours, yy * 1e3)
plt.xlabel("Hours")
plt.ylabel("mmol/l")
plt.title("Pyrite Dissolution")
plt.grid()
plt.show()

Hello, is it possible to use the phreeqpython interface to retrieve the species diffusion coefficients (or mobilities) used in the specific conductance calculation? I have tried a bit with, e.g., solution.pp.ip.get_selected_output_XXX() methods but have not been able to understand their call signatures. Thank you.

Hi,
Is there a method to get the charge balance in equivalents (and/or percent error) for a solution. I want to know if my solutions are electrically balanced and possibly manipulate an element to achieve balance. Thank you!

I am using miniconda on PopOS x64. I have tried installing the package using conda but it is not there, however phreeqpy is available [https://anaconda.org/conda-forge/phreeqpy]

Are these x2 different?

Hi,

I have a question about the run_string() function. I would like to use some features which are not yet in PhreeqPython (e.g. Exchange and Transport) and will need to use pp.ip.run_string() to run them.

However, when testing using it on a built in function like solution.equalize, I do not get the same results as when I try and run 'Equilibrium Phases' in the use pp.ip.run_string(). The solution is missing phases and I am not sure what I am doing wrong?

Thanks for any help!

##################################

#Github example 2

from phreeqpython import PhreeqPython

pp = PhreeqPython()

pp.add_master_species( element = 'Doc',
master_species = 'Doc',
alkalinity = 0,
egfw = 12)

pp.add_species(reaction = 'Doc = Doc',
log_k= 0)

#Add Jurbanite to the database
str_to_run = """
PHASES # Taken from WATEQ4f database as not present in phreeqc.dat
Jurbanite
AlOHSO4 + H+ = Al+3 + SO4-2 + H2O
log_k -3.230
"""

pp.ip.run_string(str_to_run)

solution0 = pp.add_solution({'temp' : 8.9,
'pH': 5.75,
'Ca': 0.9,
'Mg': 0.16,
'Na':0.65,
'K':0.05,
'Al':0,
'Fe(+2)': 0.19,
'Cl':0.76,
'C(+4)':3.5, #TIC
'S(6)':0.6,
'Si':0.1,
'DOC':0,
})

str_to_run ="""

EQUILIBRIUM_PHASES 0
Jurbanite 0 1

"""
pp.ip.run_string(str_to_run)

print(solution0.phases)

solution1 = solution0.copy()

#EQUILLIBRIUM PHASES

solution1.equalize(['Jurbanite'], [0], [1])

print(solution1.phases)

Hello,

There seems to be a discrepancy between the saturation indices returned by Solution.si(), and those returned by PhreeqPython.ip.get_selected_output.

If I run:

from phreeqpython import PhreeqPython
pp = PhreeqPython("/hpcdata/fijmjl/iphreeqc-3.7.3-15968/database/phreeqc.dat", debug=True)
soln = pp.add_solution(        
    {    
        "units": "ppm",
        "Mg": 516, 
        "Mn(2)": 2.6, 
    }, 
)
minerals = ["Hausmannite", "Manganite", "Pyrochroite", "Pyrolusite"]
{k: soln.si(k) for k in minerals}

I get the following SIs:

SOLUTION 0
  units ppm
  Mg 516
  Mn(2) 2.6
SAVE SOLUTION 0
END 

{'Hausmannite': -18.90189294986606,
 'Manganite': -8.963964316622022,
 'Pyrochroite': -5.823964316622019,
 'Pyrolusite': -18.003964316622024}

But using the selected output array:

from phreeqpython import PhreeqPython
pp = PhreeqPython("/hpcdata/fijmjl/iphreeqc-3.7.3-15968/database/phreeqc.dat", debug=True)
phases = soln.phases.keys()
pp.ip.run_string(
"""
SOLUTION 1
  units ppm
  Mg 516
  Mn(2) 2.6
SAVE SOLUTION 2

SELECTED_OUTPUT 1
-file selected_output.sel
-reset false
-saturation_indices Hausmannite Manganite Pyrochroite Pyrolusite
END
"""
)
{k: v for k, v in zip(*pp.ip.get_selected_output_array())}

I get different SIs for the minerals:

SOLUTION 1
  units ppm
  Mg 516
  Mn(2) 2.6
SAVE SOLUTION 2

SELECTED_OUTPUT 1
-file selected_output.sel
-reset false
-saturation_indices Hausmannite Manganite Pyrochroite Pyrolusite
END

{'si_Hausmannite': -10.902521584503347,
 'si_Manganite': -4.964278633940662,
 'si_Pyrochroite': -5.824278633940661,
 'si_Pyrolusite': -10.004435792599985}

Worth noting that the latter results match what I get from running the phreeqc desktop application. Any ideas on what might be going on?

Thank you in advance for the help!

Hi!

I would like to use a feature of a the latest version of PHREEQC ('implicit' in the TRANSPORT keyword). As far as I can tell this feature is in PHREEQC version 3.6, but phreeqpython uses Iphreeqc 3.3.7.

Is there anyway to update to the latest version of PHREEQC? to make this feature available?

Thank you!

Alex

I want to couple PHREEQC with a 3D groundwater flow and solute transport model. The 3D groundwater flow and solute transport model has API that should make this possible.

Basically I want to do this:
for each reaction/transport step:

1. calculate the geochemistry in all model cells of the 3D groundwater model using PHREEQC
2. pass these results from PHREEQC to the 3D groundwater flow model using the APIs of the models
3. Simulate non-reactive transport of the species in the 3D model
4. pass these results from the 3D groundwater flow model back to PHREEQC to the using the APIs of the models
5. repeat

Can I do this with phreeqpython? I did not see any examples of modifying variables during a run.

in phreeqc I can write

 EQUILIBRIUM_PHASES
Goethite 0 1.04e-3

which states that Geothite should be in equilibrium but that at best only 1.04e-3 mol can be dissolved in order to attain equilibrium. I am unclear on how to achieve this with the solution.equalize method. Thx for any pointers!

Hi,

I am simulating the evaporation of a lake (by removing moles of H20) in order to follow the evolution of brine element concentrations and the type of salt precipitated.

lake.desaturate enable me to allow the precipitation of salts under over-saturation sate.

I am however struggling to know how much salt has precipitated after each desaturation state.

As several minerals sharing similar elements can be simultaneously precipitating (e.g. CaSO₄ ; K₂MgCa₂(SO₄)₄•2H₂O and MgSO₄•7 H₂O ), it is very challenging to calculate the amount of precipitated salts from the evolution of the elements content in the brine, (by comparing element contents before and after the desaturation stages).

Is there a way to obtain the output of the moles of precipitated salt after each desaturation stage ?

Is there a way to retrieve the alkalinity of a solution? I examined the available properties and I cant see it.

In addition, if there is no clear method to do so, is it possible to pass a USER PRINT string and retrieve output to somehow capture it?

Thanks in advance!

Hi there! I just wanted to quickly ask: what's the function to pull out the activity of water? Clearly PHREEQC outputs this as you can see in the image below (at least with the PHREEQ database), but I just can't seem to find the right function with the python wrapper.

Thanks so much for your help with this - I am extremely grateful for this package!

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!

when I write a rate function where the addition of an element also changes the pH of the solution, the result will differ from what I get with phreeqc. A good example is the oxidation of pyrite which depends on the activity of H+ and Fe+3, and will also change both during the dissolution process.

Looking at the code, it appears that the activities will always be taken from the initial solution, and as such do not reflect the fact that the oxidation rate changes in concert with both activities.

Is there an easy way to calculate such a system with phreeqpython? I can provide code examples if need be.

Warning IPHREEQC official libs are same as VITENS ...change $PREFIX

You can compare sizes via (VItens are way bigger)
ls -sh libiphreeqc*

cd $PREFIX/local/lib [default official]

28M libiphreeqc-3.7.3.so 4.0K libiphreeqc.la
56M libiphreeqc.a 0 libiphreeqc.so

cd $PREFIX/lib [Vitens]

89M libiphreeqc-3.7.3.so 4.0K libiphreeqc.la
226M libiphreeqc.a 0 libiphreeqc.so

Extension of the IPhreeqc 3.3.7 module (Parkhurst&Appello). This extension aims to add more flexibility to the IPhreeqc module by exposing more information to the C extension. Instead of relying on SELECTED_OUTPUT this extensions enables users to directly gather information such as the pH, SC, speciation, element totals etc. of solutions.

git clone https://github.com/Vitens/VIPhreeqc
FC=gfortran-11 CC=gcc-11 CXX=g++-11 ./configure --prefix=$PREFIX
make -j4 
make install

https://github.com/Vitens/phreeqpython
pip install -U phreeqpython
It will install easily but default libs will not work as are not compiled for arm
Lib location, default names & replacement:

cd $PREFIX/lib/python3.9/sitepackages/phreeqpython/lib/
ls
rm *
ln -s $PREFIX/lib/libiphreeqc-3.7.3.so viphreeqc.so

Thus, first step is required (see above)

I ran the example included on the README.md file of your repository (also see #28) and got bad output from the Solution.total() method.

This file contains a jupyter notebook that demonstrates the discrepancy:
phreeqpython-readme.zip

This might be due to a problem with the https://github.com/Vitens/VIPhreeqc codebase.

Hi,
I'm having issues passing the correct solution numbers from PhreeqPython to Phreeq for the bits from Phreeq where I don't think there is built in functions in PhreeqPython e.g. EXCHANGE, TRANSPORTS etc.

For example, using EXCHANGE I can fix the problem by changing the solution number from 0 to 1, because in EXCHANGE I can specify what SOLUTION I want to use.

But for TRANSPORT, it requires a SOLUTION 0 as the influent and you cannot change the solution number that it takes as the influent.

Is there a way to fix this and make PhreeqPython pass a SOLUTIONS 0 to Phreeq when using the pp.ip.run_string(...) code?

Thanks you!

Code for EXCHANGE changing from SOLUTION 0 to SOLUTION 1;
from phreeqpython import PhreeqPython

pp = PhreeqPython()

pp.add_master_species( element = 'Doc',
master_species = 'Doc',
alkalinity = 0,
# gfw= 12,
egfw = 12)

pp.add_species(reaction = 'Doc = Doc',
log_k= 0)

#Add Jurbanite to the database
str_to_run = """
PHASES # Taken from WATEQ4f database as not present in phreeqc.dat
Jurbanite
AlOHSO4 + H+ = Al+3 + SO4-2 + H2O
log_k -3.230
"""

pp.ip.run_string(str_to_run)

#Create solution
solution = pp.add_solution({'temp' : 8.9,
'pH': 4.18,
'Ca': 0.3,
'Mg': 0.16,
'Na':0.65,
'K':0.05,
'Al':0.2,
'Fe(+2)': 0.11,
'Cl':0.76,
'C(4)':3.5, #TIC
'S(6)':1.03,
'Si':0.07,
'Doc':1.3,
})

CATION EXCHANGE

#Will not work if specify EXCHANGE 0 and -equilibrate 0 ?
str_to_run ="""

EXCHANGE 1
X 0.031 # See exercise manual. CEC must be in eq/L.
-equilibrate 1 # bring in equilibrium with solution 0
END

"""
pp.ip.run_string(str_to_run)

Hi,
when you calculate IS in phreeqc for a mineral you have the IS value, Log IAP and Log K. When you calculate the equilibrium phases you also have 'Moles in assemblage' (initial, final and delta). Do you know how to get those values in phreeqpython? could you help me? Thanks!

I've tried this version of phreeqpython and the version on pip. Both give this Error 126. Do you know how to fix this error?
It can run on command line and ipython in command line. However, it cannot run in an IDE like Spyder or Canopy.

import phreeqpython as phre

pp = phre.PhreeqPython()
Traceback (most recent call last):

File "", line 1, in
pp = phre.PhreeqPython()

File "C:\Users\Lohen\Miniconda2\lib\site-packages\phreeqpython\phreeqpython.py", line 15, in init
self.ip = VIPhreeqc()

File "C:\Users\Lohen\Miniconda2\lib\site-packages\phreeqpython\viphreeqc.py", line 51, in init
phreeqc = ctypes.cdll.LoadLibrary(dll_path)

File "C:\Users\Lohen\Miniconda2\lib\ctypes_init_.py", line 444, in LoadLibrary
return self._dlltype(name)

File "C:\Users\Lohen\Miniconda2\lib\ctypes_init_.py", line 366, in init
self._handle = _dlopen(self._name, mode)

WindowsError: [Error 126] The specified module could not be found

When I run

        -cells           20
        -lengths         0.00025 
        -shifts          1001		#total time = shifts x time_step
        -time_step       3.894488644294231e-07 year	#years
        -flow_direction  diffusion_only
        -boundary_conditions   constant closed
        -implicit true
        -punch_cells     1-20
        -punch_frequency 10
        -print_cells     1-20
        -print_frequency 10
	-multi_D true 2.29e-09 0.10865607270238731 0 2.7 
        -porosities 0.10865607270238731
	-dispersivities 0.001 
        END 

it gives me the error

ERROR: Unknown option.
ERROR:         -implicit true
ERROR: Unknown input in TRANSPORT keyword.
ERROR:         -implicit true
ERROR: Calculations terminating due to input errors.

while according to the manual it is (currently) supported: https://water.usgs.gov/water-resources/software/PHREEQC/documentation/phreeqc3-html/phreeqc3-56.htm#50593793_87317

This manual is for PHREEQC 3.56, I cannot deduce which phreeqc version is used within phreeqpython.

Is there any way we can use a newer version of phreeqc within phreeqpython

First of all,

Thank you for creating this package.

I was trying to get the density of a solution resulting from water-rock interaction process but it seems that there is no density class implemented in the solution.py file.

Could you please add it?

Regards,
Oscar

Rewritten self.total for faster performance (133 vs 83,6 µs ≈ 37% improvement).
Image below shows self.total_old vs. new self.total.