Lithium speciation calculation on salt lake brines with Phreeqc and Python - Tutorial

Lithium interaction with the aqueous phase on brines phase is not well understood or at least not well covered by the scientific publications. To evaluate brine composition, brine dynamics and the interaction of the different components is really a challenge with the available tools for geochemical modeling as Phreeqc since it doesn’t appear in all the databases and in most it is only associated with sulfates.

This is an applied example of chemical speciation of brines from a salt lake on the Tibetan Plateau in China that has Lithium in concentrations close to 250 mg. A Phreeqc input file is written with the concentration values and the model is run and evaluated under a Python script in a Jupyter Notebook.

Concentrations were taken from this publication:

Hydrochemistry, Distribution and Formation of Lithium-Rich Brines in Salt Lakes on the Qinghai-Tibetan Plateau.

Tutorial

Input data

You can download the input data from this link.

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.7.3-15968-x64"
chemModel.phBin = os.path.join(phPath,"bin\\phreeqc.bat")
chemModel.phDb = os.path.join(phPath,"database\\llnl.dat")
chemModel.phDb

'C:\\Program Files\\USGS\\phreeqc-3.7.3-15968-x64\\database\\llnl.dat'

Set the input and output files

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

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

Run model and show simulation data

chemModel.runModel()

Input file: ..\In\exLithium.in

Output file: ..\Out\exLithium.out

Database file: C:\Program Files\USGS\phreeqc-3.7.3-15968-x64\database\llnl.dat

              █▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀█
              ║                                            ║
              ║             * PHREEQC-3.7.3 *              ║
              ║                                            ║
              ║      A hydrogeochemical transport model    ║
              ║                                            ║
              ║                    by                      ║
              ║       D.L. Parkhurst and C.A.J. Appelo     ║
              ║                                            ║
              ║             December  2, 2021              ║
              █▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄█


Simulation 1. Initial solution 1.        /                                             

End of Run after 1.699 Seconds.

chemModel.showSimulations()

Parsing output file: 1 simulations found
Simulation 1:  Example 1a.--Speciate salt lakes. from line 20 to 239

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	0.628900	0.628900	None
B	0.107700	0.107700	None
Cl	0.124300	0.124300	None
K	0.034870	0.034870	None
Li	0.037020	0.037020	None
Na	0.246400	0.246400	None
O(0)	0.000451	0.000451	Equilibrium with O2(g)
S	0.046960	0.046960	None

iS[['Molality']].plot(kind='bar', grid=True)

<Axes: xlabel='Elements'>

dS = initSol['descriptionSolution']
dS

	Value
Parameter
pH	7.00000
pe	4.00000
Activity of water	0.97900
Ionic strength (mol/kgw)	0.54330
Mass of water (kg)	1.00000
Total carbon (mol/kg)	0.70650
Total CO2 (mol/kg)	0.70650
Temperature (°C)	25.00000
Electrical balance (eq)	-0.52880
Percent error, 100*(Cat-|An|)/(Cat+|An|)	-52.23000
Iterations	5.00000
Total H	111.99720
Total O	58.07579

initSol['redoxCouples']

	pe	Eh(volts)
Redox couple
O(-2)/O(0)	13.6053	0.8048

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+	1.248000e-07	1.000000e-07	-6.904	-7.000	-0.096	0.0
H2O	5.553000e+01	9.791000e-01	1.744	-0.009	0.000	NaN
B(-5)	0.000000e+00	NaN	NaN	NaN	NaN	NaN
BH4-	0.000000e+00	0.000000e+00	-223.712	-223.881	-0.169	NaN
B(3)	1.077000e-01	NaN	NaN	NaN	NaN	NaN

# Show molalities for all species without H+ and H20
dS[['Molality']].iloc[2:].plot(kind='bar', figsize=(24,4))

<Axes: xlabel='Species'>

# Show Carbon species 
dS[['Molality']].iloc[9:15].plot(kind='bar')

<Axes: xlabel='Species'>

# Show Chloride species 
dS[['Molality']].iloc[15:26].plot(kind='bar')

<Axes: xlabel='Species'>

# Show Lithium species 
dS[['Molality']].iloc[39:44].plot(kind='bar')

<Axes: xlabel='Species'>

# Show Sodium species 
dS[['Molality']].iloc[44:52].plot(kind='bar')

<Axes: xlabel='Species'>

Working with the saturation indices

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

	SI	log IAP	log K(298 K, 1 atm)	Description
Phase
Aphthitalite	-6.44	-10.32	-3.89	NaK3(SO4)2
Arcanite	-3.69	-5.53	-1.84	K2SO4
B	-107.82	1.74	109.56	B
B(g)	-199.11	1.74	200.84	B
B2O3	-7.46	-1.92	5.55	B2O3

sI[['SI']].plot(kind='bar', figsize=(24,4), grid=True)

<Axes: xlabel='Phase'>