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.

Tutorial

Code

# 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: https://www.usgs.gov/software/phreeqc-version-3

# 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")
chemModel.phDb

'C:\\Program Files\\USGS\\phreeqc-3.6.2-15100-x64\\database\\phreeqc.dat'

Set the input and output files

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

WindowsPath('../In/ex1a.in')

Run model and show simulation data

chemModel.runModel()

Input file: ..\In\ex1a.in

Output file: ..\Out\ex1.out

Database file: C:\Program Files\USGS\phreeqc-3.6.2-15100-x64\database\phreeqc.dat



                           * PHREEQC-3.6.2 *              

                    A hydrogeochemical transport model    

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

                            January 28, 2020              



Simulation 1. Initial solution 1.        /                                             

End of Run after 1.666 Seconds.

chemModel.showSimulations()

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
#simObj.getSimulationDict()

Get initial solution

initSol = simObj.getInitialSolution()
initSol.keys()

dict_keys(['solutionComposition', 'descriptionSolution', 'redoxCouples', 'distributionSpecies', 'saturationIndices'])

iS = initSol['solutionComposition']
iS

	Molality	Moles	Description
Elements
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)

<AxesSubplot:xlabel='Elements'>

dS = initSol['descriptionSolution']
dS

	Value
Parameter
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

initSol['redoxCouples']

	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']
dS.head()

	Molality	Activity	Log Molality	Log Activity	Log Gamma	mole V cm3/mol
Species
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))

<AxesSubplot:xlabel='Species'>

Working with the saturation indices

sI = initSol['saturationIndices']
sI.head()

	SI	log IAP	log K(298 K, 1 atm)	Description
Phase
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)

<AxesSubplot:xlabel='Phase'>

Input data

You can download the input data from the following link.