In this post, I'm going to show how we can call Phreeqc from Julia to calculate the solubility of minerals in a solution. The workflow is very useful specially when one wants to conduct a sensitivity analysis and plot the results. Here is my workflow:
- Create a Phreeqc input file in Julia
- Save it to a text file
- Call Phreeqc from Julia to run the created input file
- Read the results from the output file created by Phreeqc
- Plot the results
If you are on windows, you need to add phreeqc to the path so you can call it from the command line.
create a Phreeqc simulation file¶
The following string is a phreeqc input file that mixes 10.0 moles of anhydrite in 1 kg of demineralized water at 50 degree Celsius and 250 atm, and writes the molality of Ca ion in the aqueous solution to an output file named
SOLUTION 1 Pure water pH 7.0 temp 50.0 pressure 10 EQUILIBRIUM_PHASES 1 Anhydrite 0.0 10.0 SELECTED_OUTPUT -file anhydrite_solubility.csv -reset false -totals Ca END
You can save this code in a file names e.g.,
anhydrite_sol.phr, and run it by
if you have phreeqc installed and in your path. In windows, you can download phreeqc for windows and run it.
My purpose here is to study the effect of thermodynamic models and databases, temperature, and pressure on the predicted solubility of anhydrite in water. It can be done manually, but I like to do it automatically from a Julia script. You ca do it in Matlab or any other tools of your choice.
We first create the Phreeqc input file as an string in Julia. Then, we save that script into a file and call phreeqc to run it. We save the results in an output file and read it in Julia. Finally we plot the results. The following script shows how this can be done in Julia. We need to use the packages
PyPlot for visualization and
DataFrames for reading the output file.
using PyPlot, DataFrames
T_min = 25.0 # degC minimum temperature T_max = 100.0 # degC maximum temperature p_atm = 250 # atm pressure T_range = collect(linspace(T_min, T_max, 20)) data_base = ["phreeqc.dat", "pitzer.dat"] for db in data_base anhydrite_sol = zeros(0) for T in T_range phreeqc_input = """ DATABASE /usr/local/share/doc/phreeqc/database/$db SOLUTION 1 Pure water pH 7.0 temp $T pressure $p_atm EQUILIBRIUM_PHASES 1 Anhydrite 0.0 10.0 SELECTED_OUTPUT -file anhydrite_solubility.csv -reset false -totals Ca END """ # save it to a file ph_file=open("phreeqc_input_file", "w") write(ph_file, phreeqc_input) close(ph_file) # run it run(`phreeqc phreeqc_input_file`) # read the results ph_res=readtable("anhydrite_solubility.csv", separator = ' ') # plot the results push!(anhydrite_sol, ph_res[:Ca][end]) end IJulia.clear_output() plot(T_range, anhydrite_sol) end legend(data_base) xlabel("Temperature(C)") ylabel("Anhydrite solubility [mol/kgw]");
You probably need to change the path to the Phreeqc databeses in line 11 of the above code. Play with the value of pressure and see what happens. It is truely fascinating that increasing pressure can affect the solubility that much. Don't forget that it can only be captured by the Pitzer model in Phreeqc.
There is another way of calling Phreeqc from Julia that I will discuss later.