Seawater speciation modeling with Phreeqc coupled with Python and Pandas – Tutorial


The speciation modeling allows to calculate the distribution of aqueous species in a solution. Phreeqc is capable to simulate this speciation calculation and we are going to demonstrate this capability on a study case of aqueous species in seawater.

We have done a tutorial for the speciation modeling of seawater with Phreeqc that runs under Python in a Jupyter Lab enviroment. The code can run the Phreeqc executable, define the databases and stablish the output files. Results from simulation are available as Pandas dataframes and plots are made for the main components and the distribution of saturation indices.



# This tutorial works on Windows and might work on Linux althought it have not been tested. 
# You need to have the batch version of Phreeqc installed
# Download Phreeqc from:

# import required packages and classes
import os
from workingTools import phreeqcModel
from pathlib import Path

Create a Phreeqc object, define executables and databases

#define the model object
chemModel = phreeqcModel()

#assing the executable and database
phPath = "C:\Program Files\USGS\phreeqc-3.6.2-15100-x64"
chemModel.phBin = os.path.join(phPath,"bin\phreeqc.bat")
chemModel.phDb = os.path.join(phPath,"database\phreeqc.dat")
'C:\Program Files\USGS\phreeqc-3.6.2-15100-x64\database\phreeqc.dat'

Set the input and output files

chemModel.inputFile = Path("../In/")
chemModel.outputFile = Path("../Out/ex1.out")

Run model and show simulation data

Input file:

Output file: ..Outex1.out

Database file: C:Program FilesUSGSphreeqc-3.6.2-15100-x64databasephreeqc.dat

                           * PHREEQC-3.6.2 *              

                    A hydrogeochemical transport model    

                     D.L. Parkhurst and C.A.J. Appelo     

                            January 28, 2020              

Simulation 1. Initial solution 1.        /                                             

End of Run after 1.666 Seconds.
Parsing output file: 1 simulations found
Simulation 1:  Example 1a.--Speciate seawater. from line 18 to 261

Show simulation components

simObj = chemModel.getSimulation(1)
Simulation content:
        Initial solution calculation: True
        Batch reaction calculations: False
        Number of reactions steps: 0
# In case you need more information about the line interval of the output file

Get initial solution

initSol = simObj.getInitialSolution()
dict_keys(['solutionComposition', 'descriptionSolution', 'redoxCouples', 'distributionSpecies', 'saturationIndices'])
iS = initSol['solutionComposition']
Molality Moles Description
Alkalinity 2.406000e-03 2.406000e-03 None
Ca 1.066000e-02 1.066000e-02 None
Cl 5.657000e-01 5.657000e-01 None
Fe 3.711000e-08 3.711000e-08 None
K 1.058000e-02 1.058000e-02 None
Mg 5.507000e-02 5.507000e-02 None
Mn 3.773000e-09 3.773000e-09 None
N(-3) 1.724000e-06 1.724000e-06 None
N(5) 4.847000e-06 4.847000e-06 None
Na 4.854000e-01 4.854000e-01 None
O(0) 4.376000e-04 4.376000e-04 Equilibrium with O2(g)
S(6) 2.926000e-02 2.926000e-02 None
Si 7.382000e-05 7.382000e-05 None
iS[['Molality']].plot(kind='bar', grid=True)
dS = initSol['descriptionSolution']
pH 8.220000
pe 8.451000
Specific Conductance (µS/cm, [0-9]5°C) 52630.000000
Density (g/cm³) 1.023230
Volume (L) 1.012820
Activity of water 0.981000
Ionic strength (mol/kgw) 0.674700
Mass of water (kg) 1.000000
Total carbon (mol/kg) 0.002182
Total CO2 (mol/kg) 0.002182
Temperature (°C) 25.000000
Electrical balance (eq) 0.000794
Percent error, 100*(Cat-|An|)/(Cat+|An|) 0.070000
Iterations 7.000000
Total H 111.014700
Total O 55.630540
pe Eh(volts)
Redox couple
N(-3)/N(5) 4.6750 0.2766
O(-2)/O(0) 12.4062 0.7339

Working with the distribution of species

dS = initSol['distributionSpecies']
Molality Activity Log Molality Log Activity Log Gamma mole V cm3/mol
H+ 7.984000e-09 6.026000e-09 -8.098 -8.220 -0.122 0.00
H2O 5.551000e+01 9.806000e-01 1.744 -0.009 0.000 18.07
C(4) 2.182000e-03 NaN NaN NaN NaN NaN
HCO3- 1.485000e-03 1.003000e-03 -2.828 -2.999 -0.171 26.98
MgHCO3+ 2.560000e-04 1.610000e-04 -3.592 -3.793 -0.201 5.82
# Show molalities for all species without H+ and H20
dS[['Molality']].iloc[2:].plot(kind='bar', figsize=(24,4))

Working with the saturation indices

sI = initSol['saturationIndices']
SI log IAP log K(298 K, 1 atm) Description
Anhydrite -0.93 -5.20 -4.28 CaSO4
Aragonite 0.61 -7.73 -8.34 CaCO3
Calcite 0.75 -7.73 -8.48 CaCO3
Chalcedony -0.52 -4.07 -3.55 SiO2
Chrysotile 3.36 35.56 32.20 Mg3Si2O5(OH)4
sI[['SI']].plot(kind='bar', figsize=(24,4), grid=True)

Input data

You can download the input data from the following link.

Level 6 Domain Names are FREE – Register Now.

Related Posts