HPGeochemistry/HPx Short Course

Advanced Modelling of Water Flow and Biogeochemical
Transport in Porous Media Using the
HPx Software Packages
Prague, 03-04 September 2020

Authors: Diederik Jacques
Affiliations: EGA, W&D
Publication date: 2020-08-31

© SCK CEN – Publication date 2020-08-31
Stichting van Openbaar Nut - Fondation d'Utilité Publique ‐ Foundation of Public Utility
Registered Office:Avenue Herrmann Debroux 40 - 1160 Brussel – Belgium
Research Centres:
Boeretang 200 - 2400 Mol - Belgium
Chemin du Cyclotron 6 - 1348 Ottignies-Louvain-la-Neuve - Belgium
www.sckcen.be

Copyright rules

All property rights and copyright are reserved to SCK CEN. This document contains data, information and formats for
dedicated use only and may not be communicated, copied, reproduced, distributed or cited without the explicit written
permission of SCK CEN. Any infringement to this rule is illegal and entitles to claim damages from the infringer, without
prejudice to any other right e.g. in case of granting a patent or registration in the field of intellectual property.

Other Research Output SCK CENSCK CEN/39560003

p. 2 / 2

ISC: Public
Last saved on 31-08-20

Contents

Contents ........................................................................................................................................................................................ i
List of Figures ............................................................................................................................................................................. v
List of Tables .............................................................................................................................................................................. ix
List of Code Extracts ................................................................................................................................................................ x
PREFACE..........................................................................................................................................................................................xiii
INTRODUCTION .............................................................................................................................................................................1
Goals .........................................................................................................................................................................................1
Goals of the course .............................................................................................................................................................1
Reactive Transport Modelling ..............................................................................................................................................1
The Modelling Flowchart ..................................................................................................................................................1
Coupled Reactive Transport Code ................................................................................................................................2
Environmental and Engineered Applications............................................................................................................5
Course software .........................................................................................................................................................................7
Detailed course outline ..........................................................................................................................................................9
About the notes...................................................................................................................................................................... 11
Text font and colouring .................................................................................................................................................. 11
Code examples .................................................................................................................................................................. 12
Syntax examples ................................................................................................................................................................ 12
Session 1 – Aqueous speciation calculations ................................................................................................................... 15
Learning outcomes................................................................................................................................................................ 16
Geochemical System............................................................................................................................................................. 16
Geochemical system ........................................................................................................................................................ 16
Aqueous Speciation .............................................................................................................................................................. 19
Problem Definition ........................................................................................................................................................... 19
Conceptual Model ............................................................................................................................................................ 19
Mathematical world – Aqueous speciation ................................................................................................................. 19
Principles of thermodynamic equilibrium ............................................................................................................... 19
Activity correction models for aqueous species ................................................................................................... 24
Saturation state and saturation index ...................................................................................................................... 25
Computer world - Interface – HPGeoChemistry ........................................................................................................ 27
Learning Outcomes .......................................................................................................................................................... 27
HPGeochemistry graphical user interface ............................................................................................................... 27
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

i
ISC: PublicDatabase Manager ........................................................................................................................................................... 28
Project Manager ................................................................................................................................................................ 30
Computer world – PHREEQC Input ................................................................................................................................. 31
Implementing the model ............................................................................................................................................... 31
Phreeqc conventions ....................................................................................................................................................... 32
SOLUTION_MASTER_SPECIES ...................................................................................................................................... 34
SOLUTION_SPECIES ......................................................................................................................................................... 35
PHASES ................................................................................................................................................................................. 37
SOLUTION............................................................................................................................................................................ 38
Speciation model: Sea and ground water .................................................................................................................... 40
Problem Definition ........................................................................................................................................................... 40
Hands-on ............................................................................................................................................................................. 40
Explore the output............................................................................................................................................................ 45
Answer................................................................................................................................................................................... 48
Charge balancing and target saturation index .......................................................................................................... 49
Charge balance a solution............................................................................................................................................. 49
Adjust to a saturation index of a phase ................................................................................................................... 50
Session 2 – Modelling Reaction Paths ................................................................................................................................ 51
Learning outcomes................................................................................................................................................................ 53
Reaction path: Mixing two solutions .............................................................................................................................. 53
Learning outcome ............................................................................................................................................................ 53
Problem definition............................................................................................................................................................ 53
Conceptual model ............................................................................................................................................................ 53
Implementing the model ............................................................................................................................................... 54
Hands-on: Mixing two solutions ................................................................................................................................. 55
Reaction path: Sliding concentration path .................................................................................................................. 55
Learning outcomes .......................................................................................................................................................... 55
Speciation of Cd as a function of Cl concentration ............................................................................................ 56
Reaction path: Sliding activity path ................................................................................................................................ 66
Learning outcomes .......................................................................................................................................................... 66
Solubility of Fe(OH)3 as a function of pH ................................................................................................................ 66
Session 3 – Reactive transport modelling ......................................................................................................................... 87
Learning outcomes................................................................................................................................................................ 88
Reaction transport modelling ........................................................................................................................................... 88
Coupled Reactive Transport ......................................................................................................................................... 88
Examples of reactive transport model studies ...................................................................................................... 91
Mathematical Model .....................................................................................................................................................107
Diederik Jacques (djacques@sckcen.be)
SCK CEN/39560003

ii
ISC: PublicReactive transport modelling – The HPx model ......................................................................................................109
Reactive transport modelling (HP1) – Equilibrium mineral precipitation and dissolution .....................114
Dissolution and precipitation of gypsum and calcite .......................................................................................114
Transport of Cd ................................................................................................................................................................121
Effect of a higher Cl concentration on Cd mobility ...........................................................................................124
Session 4 – Equilibria at the solid-liquid interface .......................................................................................................129
Learning outcomes..............................................................................................................................................................130
Geochemical modelling – Cation Exchange ..............................................................................................................130
Mathematical world – Cation Exchange ................................................................................................................130
Implementing the model .............................................................................................................................................132
Speciation model - Equilibration of an exchanger with a solution .............................................................135
Reaction path – Exchange isotherm ........................................................................................................................140
Reaction path – Changing exchange composition during evaporation ...................................................145
Reactive Transport Modelling (HP1) – Cation Exchange ......................................................................................149
Transport and Cation Exchange – Single Pulse...................................................................................................149
Transport and Cation Exchange – Multiple Pulse ..............................................................................................155
Session 5 – Kinetics ..................................................................................................................................................................159
Learning Outcomes .............................................................................................................................................................160
Numerical world - Implementing kinetics..................................................................................................................160
Implementing the rate equation: RATES ...............................................................................................................160
Kinetic reactants: KINETICS RATES ...........................................................................................................................161
Examples of RATES and KINETICS .................................................................................................................................162
Reaction Path – Kinetic dissolution of a mineral ................................................................................................162
Reaction Path – Biodegradation using Monod kinetics ..................................................................................171
First-order kinetic degradation network: Degradation of PCE ..........................................................................176
Learning outcomes ........................................................................................................................................................176
Background .......................................................................................................................................................................176
Geochemical Modelling (HPg) – Modelling first-order degradation .........................................................179
Geochemical Modelling (HPg) – Modelling first-order degradation network .......................................185
Reactive Transport Modelling (HP1) – Leaching and biodegradation of PCE ........................................186
Reactive Transport Modelling (HP1) – Leaching of initial PCE source .......................................................190
Reactive Transport Modelling (HP1) – Leaching of initial PCE source – Source trapped in immobile
water region ......................................................................................................................................................................195
References....................................................................................................................................................................................201
Index ...............................................................................................................................................................................................205

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

iii
ISC: PublicDiederik Jacques (djacques@sckcen.be)
SCK CEN/39560003

iv
ISC: Public

List of Figures
Figure 1 – General flowchart for the modelling environmental- ................................................................................2
Figure 2 - Coupling of different processes - bold text indicates the minimum requirement to have acoupled reactive transport model ...........................................................................................................................................3
Figure 3 – Example of a continuum representation of a three-phase porous system. Chemicalcomponents will be distributed between these phases according to equilibrium (double arrows) andkinetic (single arrows) reactions. The system will react when mass enters or leave the system (blue arrows
at the left hand side) and might be buffered by a gas or solid phase. ....................................................................4
Figure 4 - Concept of a specific reactive transport code (left – conventional model code) and a genericreactive transport code (right – modelling environment) (after Meysman et al. (2003a). ................................5
Figure 5 – Spatial scope of reactive transport models for environmental applications: Watershedcatchment scale (km) (left), macro-scale / profiles (m – mm, center), and micro-scale /pore scale (µm).
From https://www.juliaperdrial.com/research-projects.html after Bao et al. (2017) ...........................................6
Figure 7 – HPx coupling HYDRUS software with PHREEQC. .........................................................................................7
Figure 8 – Processes solved by respectively the HYDRUS software (left) and PHREEQC (right) and theexchange of information of state variables between them. .........................................................................................8
Figure 9 – Schematic representation of different models based on thermodynamic and kinetic
geochemical modelling (from Zhu and Anderson (2002)) ......................................................................................... 10
Figure 10 Representation of a geochemical system. Left: Beaker model (Deutsch and Siegel 1997). ..... 17
Figure 11 – Scheme of possible speciation and reactions of Hg in a soil system (Leterme et al. 2014). . 18
Figure 12 – Schematic representation of the model to calculate aqueous speciation (left) and key
elements of the computer model......................................................................................................................................... 19
Figure 13 – Main form of HP Geochemistry ..................................................................................................................... 28
Figure 14 – Database Manager ............................................................................................................................................. 28
Figure 15 – Project Manager .................................................................................................................................................. 31
Figure 16 - Implementation of the model for aqueous speciation......................................................................... 32
Figure 17 – Extracts of SOLUTION_SPECIES from the phreeqc.dat database (PHREEQC 3.6.2) .............. 36
Figure 18 – Extracts of PHASES from the phreeqc.dat database (PHREEQC 3.6.2). Lines 903-930 showminerals; Lines 1138-1158 show gases. ............................................................................................................................. 37
Figure 19 – View of the Input Work Space ....................................................................................................................... 41
Figure 20 – Text editor in the Input Work Space ........................................................................................................... 42
Figure 21 – Illustration of the different insert options. The grey text between angle brackets should be
replaced by integers, float, text or strings. ....................................................................................................................... 42
Figure 22 - SOLUTION_MASTER_SPECIES data block in the phreeqc.dat thermodynamic database. Line
14 illustrates that the element name to be used in SOLUTION data block should be Na (first column,
element name), and not Na+ (second column, master species). Line 30 shows the master species
Alkalinity. ........................................................................................................................................................................................ 44
Figure 23 – Output work space with the phreeqc.out file. At the left hand-side, there is a navigation tree.
............................................................................................................................................................................................................ 45
Figure 24 – Phreeqc.out – Composition of sea water. .................................................................................................. 46
Figure 25 - Phreeqc.out – Description of the sea water solution. The red rectangle indicates the percent
error on the electrical balance (Q1); the green rectangle indicates the ionic strength (Q2). ....................... 46
Figure 26 – Phreeqc.out - Distribution of species for sea water. The red rectangles indicate the molality
and activity of two aqueous species (Q3). The blue rectangles indicate the log10 of the activity correction
factor (Q4). The green rectangles indicate the most abundant species of some of the components (Q5).
The yellow rectangles indicate the log10 activity of some species needed to calculate a mass action
constant (Q6). ............................................................................................................................................................................... 47
Figure 27 – Phreeqc.out: Saturation indices for sea water. The red rectangle indicates the saturation index
of aragonite and calcite (Q7). ................................................................................................................................................. 48
Figure 28 – Variants of titration reaction paths adding reactants (batch reactors in chemical engineering,
Levenspiel (1999). ....................................................................................................................................................................... 52
Figure 29 – Schematic representation of mixing two aqueous phases. ................................................................ 54
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

v
ISC: PublicFigure 30 - Implementation of the model for mixing two solutions. ..................................................................... 54
Figure 31 – Schematic representation of a sliding concentration path with fixed activity. ........................... 56
Figure 32 – Implementation of a sliding concentration path with fixed pH ....................................................... 57
Figure 33 – Output file obtained with Code 23 (for a given geochemical model) opened in
HPGeochemistry. ......................................................................................................................................................................... 60
Figure 34 – Output file workspace: Graph ribbon.......................................................................................................... 63
Figure 35 – Graph work space – General View ................................................................................................................ 64
Figure 36 - Speciation of Cd as a function of Cl with CdTotal = 1 mM and a pH = 3. ....................................... 65
Figure 37 – Schematic representation of a sliding activity path to calculate the solubility of a mineral. 67
Figure 38 – Implementation of a sliding activity path.................................................................................................. 67
Figure 39 - Form to define the inline variables ............................................................................................................... 70
Figure 40 - Form to define the values of the inline variables.................................................................................... 71
Figure 41 – Help tree for SELECTED_GRAPH .................................................................................................................. 77
Figure 42 – Solubility of Fe(OH)3 as a function of pH in a 0.1 M NaF background solution. ....................... 85
Figure 43 – Speciation of Fe(III) as a function of pH in a 0.1 M NaF background solution........................... 86
Figure 44 – Schematic representation of a 1D reactive transport problem involving equilibrium mineraldissolution and precipitation reactions. ...........................................................................................................................114
Figure 45- The HP1 - Print and Punch Controls dialogue window ...................................................................116
Figure 46 – The Default Database for Project Group dialogue window. .......................................................117
Figure 47 - Profiles of C concentrations (Example T-Em-01) ..................................................................................121
Figure 48 – Profiles of percentage of Cd in solution (Example T-Em-02) ..........................................................124
Figure 49 – Schematic representation of a 1D reactive transport problem involving equilibrium mineraldissolution and precipitation reactions with a time-varying upper solute concentration boundary
condition. .....................................................................................................................................................................................125
Figure 50 – Time series of Cd (top) and profiles of otavite (bottom) for transport without (left, T – SS –Em – 02) and with (right, T – SS – Em – 03) higher Cl concentrations. ................................................................127
Figure 51 – Time series of Cl concentrations for T – SS – Em – 03. ......................................................................127
Figure 52 – Examples of isotherms for simple exchangers (from Pabalan and Bertetti (2001)) ................132
Figure 53 – Implementation of the model for cation exchange processes. ......................................................132
Figure 54 – Illustration of the implicit definition of EXCHANGE (see text for explanation) ........................135
Figure 55 – Schematic representation of the model to calculate cation exchange isotherms (left) and key
elements of the computer model.......................................................................................................................................140
Figure 56 – Ca/Na exchange isotherms at three NaCl concentrations (example RE_ex 01) .......................144
Figure 57 – Schematic representation of the model to calculate cation exchange during evaporation ofwater (left) and key elements of the computer model.. ............................................................................................145
Figure 58 – Equivalent fractions of Na and Ca in the pore water (left axis) and ionic strength (right axis)of a soil in equilibrium with an exchange complex during a drying event. .......................................................148
Figure 59 – Equivalent fractions of Na and Ca on the exchanger (left axis) and ionic strength (right axis)of a soil in equilibrium with an exchange complex during a drying event. .......................................................149
Figure 60 – Schematic representation of a 1D reactive transport problem involving cation exchange
reactions. ......................................................................................................................................................................................150
Figure 61 – Breakthrough curves of Na, K, Cl and Ca at the bottom of the profile during steady-statetransport with cation exchange (example T-CEC-01) .................................................................................................153
Figure 62 – Breakthrough curves (left) and equivalent sorbed fraction (right) at 4 cm (top) and 8 cm
(bottom) for example T-CEC-01. .........................................................................................................................................155
Figure 63 – Schematic representation of a 1D reactive transport problem involving cation exchange
reactions with a time varying upper solute boundary condition...........................................................................156
Figure 64 – Time series of K at different observation nodes during transport with variable boundary
conditions (Example T-CEC-02). ..........................................................................................................................................158
Figure 65 – Breakthrough curves of Cl, Na, K, and Ca during transport with variable boundary conditions(Example T-CEC-02). ................................................................................................................................................................158
Figure 66 – Implementation of the model for kinetic processes. ..........................................................................160
Diederik Jacques (djacques@sckcen.be)
SCK CEN/39560003

vi
ISC: PublicFigure 67 – (Left) Dissolution rate of kaolinite as function of DG (Nagy et al. 1991); (Right) Dissolution
rate of K-Feldspar as a function of DG (Gautier et al. 1994). ...................................................................................163
Figure 68 – FT as a function of saturation sate for 2 values of n (m=1 in Eq. (27)) for dissolution (Ω < 1).
..........................................................................................................................................................................................................164
Figure 69 – Dissolution rate of quartz as a function of temperature (Dove 1995) .........................................165
Figure 70 – Examples of pH-dependency of steady-state dissolution rates of quartz (left) and albite (right)(Walther 1996). ...........................................................................................................................................................................166
Figure 71 – Schematic representation of the model to calculate kinetic dissolution of a mineral. .........167
Figure 72 – Quartz dissolution at pH 8.5 (left) and pH 10 (right). Si concentration (top), saturation index
of quartz (middle) and dissolution rate of quartz (bottom) as a function of time. ........................................170
Figure 73 – Defining different cases using the Global Variable form........................................................171
Figure 74 – Observed and modelled phenol concentrations during a methanogenic (bio)degradation
experiment (Bekins et al. 1998) ...........................................................................................................................................171
Figure 75 – The biodegradation rate as a function of substrate concentration defined by the Monod
kinetic rate equation (from Appelo and Postma (2005). Parameters for the equation are: kmax = -1.39 mg
phenol/L/day and k0.5 = 1.7 mg phenol / L. ...................................................................................................................172
Figure 76 – Schematic representation of the model to calculate biodegradation in the aqueous phase.
..........................................................................................................................................................................................................173
Figure 77 – Dialogue window to add Thermodynamic data – Aqueous Components .................................174
Figure 78 – Phenol concentration as a function of time ...........................................................................................176
Figure 79 - Perchloroethylene (PCE) degradation pathway (after Schaerlaekens et al. (1999)) ................177
Figure 80 - Degradation pathway of PCE using first-order rate constants. .......................................................177
Figure 81 - Profiles after 200d for the PCE biodegradation problem. (a) Comparison of HP1 results fordifferent values of ωs with the analytical solution for PCE, (b) – (d) comparison of HP1 simulations withωs = 0.1 with the analytical solution for 7 components ............................................................................................179
Figure 82 – PCE concentration as a function of time. ................................................................................................183
Figure 83 – PCE concentration as a function of time on a semi-log plot. ..........................................................183
Figure 84 - PCE and TCE concentrations as a function of time (stable TCE). ....................................................184
Figure 85 – PCE and TCE concentrations as a function of time. .............................................................................184
Figure 86 – DCE, VC and Eth concentrations as a function of time. .....................................................................186
Figure 87 – Schematic representation of a 1D reactive transport problem involving a kinetic reaction
network (boundary condition problem). .........................................................................................................................187
Figure 88 – Profiles (top) and time series (bottom) of some selected degradation products from PCE forsteady-state flow conditions. ...............................................................................................................................................189
Figure 89 – Outflow concentrations of PCE and its degradation products during steady-state flow
conditions. ...................................................................................................................................................................................190
Figure 90 – Schematic representation of a 1D reactive transport problem involving a kinetic reaction
network (Initial value problem). ..........................................................................................................................................191
Figure 91 – Profiles of aqueous PCE (left) and immobile PCE (right) during leaching of an initiallycontaminated soil profile. ......................................................................................................................................................194
Figure 92 – Outflow concentrations of PCE and its degradation products during leaching of an initiallycontaminated soil profile. ......................................................................................................................................................194
Figure 93 – Total amount of PCE in the aqueous phase (top left), immobile phase (top right) and summed
amount (bottom) in the source zone (0-50 cm), pristine soil (50-200 cm) and total soil profile duringleaching of an initially contaminated soil profile. ........................................................................................................195
Figure 94 – Schematic representation of a 1D reactive transport problem involving a kinetic reaction
network (Initial value problem, mobile-immobile solute transport). ...................................................................196
Figure 95 – Profiles of PCE in the mobile region (top left), in the immobile region (top right), the averaged
concentration over both regions (bottom left), and the PCE in the immobile phase in the immobile region(bottom right) .............................................................................................................................................................................198
Figure 96 – Concentration at the bottom node in the mobile (top) and immobile (bottom) region. ....199

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

vii
ISC: PublicFigure 97 – Cumulative PCE degradation in different soil compartments during leaching of PCE trappedin the immobile soil region (top plots and bottom left). The bottom right plot shows the cumulative
leaching of PCE from the soil profile.................................................................................................................................200

Diederik Jacques (djacques@sckcen.be)
SCK CEN/39560003

viii
ISC: Public

List of Tables
Table 1 – Text font and colouring of PHREEQC related input as used in the lecture notes and the graphicaluser interface. ............................................................................................................................................................................... 11
Table 2 - Conventions for syntax .......................................................................................................................................... 12
Table 3 – Species in a geochemical system with Ca, C, O and H ............................................................................. 17
Table 4 - Composition of sea and ground water (units mg/l, alkalinity defined as HCO3-)......................... 40
Table 5 – Equations expressing the pH dependency of the steady-state dissolution rate for someminerals.........................................................................................................................................................................................166
Table 6 - Definition of parameters and their values for the PCE biodegradation problem (from Case 1
and 2 in Sun et al. (2004)). Rate parameters are for a reference temperature of 20°C.................................177

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

ix
ISC: Public

List of Code Extracts
Code 1 - Code example ............................................................................................................................................................ 12
Code 2 – Example of conventions in PHREEQC / HPGeochemistry ........................................................................ 33
Code 3 – Example of conventions in PHREEQC / HPGeochemistry ........................................................................ 34
Code 4 – Syntax and example of SOLUTION_MASTER_SPECIES ......................................................................... 35
Code 5 – Syntax of SOLUTION_SPECIES ........................................................................................................................ 36
Code 6 – Syntax of PHASES .................................................................................................................................................... 37
Code 7 – Definition of a pure solution ............................................................................................................................... 38
Code 8 – Syntax of SOLUTION – Properties of the solution...................................................................................... 38
Code 9 – SOLUTION composition........................................................................................................................................ 38
Code 10 – SOLUTION definition with different concentration units ...................................................................... 39
Code 11 – Examples to charge balance a solution ........................................................................................................ 49
Code 12 – Syntax of a solution when using pH to charge balance ........................................................................ 49
Code 13 – Example of a SOLUTION definition in equilibrium with different phases ...................................... 50
Code 14 – Syntax of a solution in equilibrium with a solid ........................................................................................ 50
Code 15 – Mix solution 1 and solution 2 with a fraction of, respectively, 0.5 and 0.25 of the original
solutions. ........................................................................................................................................................................................ 54
Code 16 – Syntax of MIX.......................................................................................................................................................... 55
Code 17 – Definition of an EQUILIBRIUM_PHASES. .................................................................................................. 57
Code 18 – Syntax of EQUILIBRIUM_PHASES. ............................................................................................................. 57
Code 19 – Definition of a phase to fix pH of a system. ............................................................................................... 57
Code 20 – Definition of an REACTANT with reactants added incrementally in equal amounts. ................ 58
Code 21 – Syntax of REACTION. ......................................................................................................................................... 58
Code 22 – Definition of an REACTANT with an explicit list of added reactants. ............................................... 58
Code 23 – Definition of a SELECTED_OUTPUT keyword. .......................................................................................... 59
Code 24 – Examples of enhanced text options in gnuplot – Superscripts and subscripts. ........................... 63
Code 25 – Syntax of GLOBAL_VARIABLE ........................................................................................................................ 68
Code 26 - Example of the of inline variables defined with GLOBAL_VARIABLES ........................................... 69
Code 27 – Examples of defining a GLOBAL_VARIABLE with simple numerical expressions. ..................... 69
Code 28 – Syntax of COUNTER and COUNTERINCR ...................................................................................................... 71
Code 29 – Example of use of COUNTER and COUNTERINCR .................................................................................... 72
Code 30 – Syntax of INCLUDE$............................................................................................................................................ 72
Code 31 – Syntax of REPEAT ................................................................................................................................................. 73
Code 32 – Syntax of USER_PUNCH ...................................................................................................................................... 74
Code 33 – Syntax of SELECTED_GRAPH – XY Plot........................................................................................................ 75
Code 34 – Definition of three plots using SELECTED_OUTPUT .............................................................................. 76
Code 35 – Curve options for formatting a given curve. .............................................................................................. 76
Code 36 – Illustration of using GLOBAL_VARIABLE as an inline variable. ......................................................... 79
Code 37 – Illustration of using COUNTER as an inline variable................................................................................. 80
Code 38 – Illustration of inline expressions, COUNTERINCR and REPEAT .......................................................... 82
Code 39 – Syntax and example of EXCHANGE_MASTER_SPECIES.....................................................................133
Code 40 – Definition of half reactions for cation exchange reactions on the exchanger X using
EXCHANGE_SPECIES. ............................................................................................................................................................133
Code 41 – Implicit definition of an EXCHANGE. ............................................................................................................134
Code 42 – Explicit definition of an EXCHANGE. .............................................................................................................135
Code 43 – Illustration and syntax of SAVE and USE..................................................................................................136
Code 44 – Input file for Example Ce_ex01. .....................................................................................................................137
Code 45 – SELECTED_OUTPUT and PUNCH keyword. ..............................................................................................141
Code 46 – Examples of adding a dataset index in SELECTED_GRAPH to plot a specific dataset in a multidataset file ...................................................................................................................................................................................142
Diederik Jacques (djacques@sckcen.be)
SCK CEN/39560003

x
ISC: PublicCode 47 – Syntax and examples of –fig_multi as used in SELECTED_GRAPH to make time series or
profile plots during reactive transport simulations. ....................................................................................................154
Code 48 – Typical syntax of RATES. ..................................................................................................................................161
Code 49 – Basic syntax of the KINETICS block ..........................................................................................................162
Code 50 – Script for modelling quartz dissolution using Eq. (32). ........................................................................168
Code 51 – KINETIC block for Quartz dissolution .......................................................................................................168
Code 52 – Script for modelling Phenol biodegradation using the Monod equation ...................................174
Code 53 - KINETICS data block for Phenol degradation using Monod kinetics ..........................................175
Code 54- Script for modelling PCE degradation using first-order degradation kinetics .............................181
Code 55 – Syntax of –horizon in REACTIVE_TRANSPORT_OUTPUT. .............................................................192

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

xi
ISC: PublicDiederik Jacques (djacques@sckcen.be)
SCK CEN/39560003

xii
ISC: PublicPREFACE
Soil and groundwater pollution is an ever-increasing, worldwide problem. Tens of billions of dollars are
spent each year to remediate groundwater pollution, and to limit or prevent future contamination of thesubsurface. Most subsurface pollution problems stem from activities involving the unsaturated (vadose)zone between the soil surface and the groundwater table. The unsaturated zone hence provides the best
opportunities to limit or prevent groundwater pollution. Once contaminants enter groundwater,
pollution is essentially irreversible, or can be remediated only with extreme costs.
Numerical modeling is becoming an increasingly important tool for analyzing complex problems
involving water flow and contaminant transport in the unsaturated zone. Coupled multicomponent
reactive transport models for the so-called critical zone are essential tools to analyze nonlinear
interactions between hydrological, geochemical and (micro)biological processes and to link the criticalzone to other spheres of the earth system such as the hydrosphere, atmosphere, biosphere, andgeosphere. This course is designed to familiarize participants with the principles and mathematicalanalysis of coupled reactive transport models, and the application of state-of-the-art numerical codes
to site-specific subsurface flow and transport problems.

Course objectives
The objective of the course is to familiarize the participants with HPx software to implement reactive
transport problems. As geochemical model implementation forms an important component of reactive
transport modelling, the participant will also learn to define geochemical models using the PHREEQC
geochemical model. A lot of examples and hands-on will illustrate how to implement specificgeochemical and reactive transport models in an efficient way. Only a few theoretical lectures will begiven.

Course software
The graphical user interface HYDRUS-1D is used to implement the reactive transport model.
HPGeochemistry is used to implement the geochemical models.

Course developer and instructor

Dr. Diederik Jacques is group leader of the unit Engineered and Geosystem analysis from the Expert

group Waste & Disposal, Institute of Environment, Health and Safety of the Belgian Nuclear ResearchCentre (SCK•CEN) in Mol, Belgium. He received a B.S. and M.S. in Bio-engineering land and forest
management at the Catholique University of Leuven, Belgium, a Master of Statistic at the Catholique
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

xiii
ISC: PublicUniversity of Leuven, and a Ph.D. in soil physics at the Catholique University of Leuven. His expertise isin modeling water flow and solute transport in unsaturated porous media including characterizing spatial
variability and estimating parameters. He is working on different aspects of coupling unsaturated water
flow, solute transport and geochemical reaction including the development and testing of the coupledcode HP1, application to (long-term) solute transport in soils and interaction between different systems(clay – concrete, or soil – concrete). Other field of expertise, among others, is safety and performanceassessments of radioactive waste disposal sites, and long-term durability of cementitious materials.
Contact: djacques@sckcen.be.

Diederik Jacques (djacques@sckcen.be)
SCK CEN/39560003

xiv
ISC: Public

INTRODUCTION
Goals
•
•
•
•

Define the goals of the course
Understand the scope of geochemical and reactive transport models
Understand the general workflow from problem definition to model implementation
Learn about the course software

Goals of the course
Theoretical goals•••

Identify the underlying theoretical models of input definitions in the numerical model
Understand the major elements and assumptions in the underlying theoretical andmathematical models
Understand the major principles of the mathematical and numerical models

This course is not a course on geochemistry or flow and transport phenomena in porous media. It is
expected that course participants have some basic knowledge in these domains. During the course, wewill make reference to some theory and mathematical models and link them to specific options you
select or commands you give in the input for the numerical models. The course will also not give details
on numerical methods to solve the mathematical equations, although we will highlight at some points
the major techniques and assumptions of the numerical methods.

Practical goals••••

Get familiar with the graphical user interfaces (HPGeochemistry and HYDRUS-1D) to
implement geochemical and reactive transport modelling projects.
Be able to formulate a problem definition in the language for the computer code
Be able to extract the relevant simulation output
Interpret the simulation output in view of the problem definition

An important point here is get familiar – within this course, only a small part of all options will be covered;and this for different levels – possible physical problems that can implemented, available conceptualmodels, available theoretical/mathematical models, options and commands in the graphical userinterfaces or the input files (both at the pre-processing (input) as the post-processing site). However, we
hope that we can give a solid basis during the course for further exploration of the broad range ofpossibilities of the course software.

Reactive Transport Modelling
The Modelling Flowchart
Figure 1 shows a general flowchart of different steps when modelling a given system:•

Physical world – For an environmental-engineered system, one may seek for understanding
and predicting its current state, its evolution towards the current state or its evolution in the
future from the current state. The physical world enters the modelling flow chart via a clear
description and a problem definition

Introduction•••

Conceptual world – The complexity of the physical system is captured in choices on
representation of key properties (to handle for example physical and geochemical
heterogeneities) and key processes (flow, transport) and geochemical reactions.
Mathematical world – Based on assumptions on applicable theories and laws, a set ofmathematical equations are defined for calculating the processes. The chemical components
and their relations by a reaction network together with thermodynamic and other parameters
is a second element in the mathematical world. Finally, initial state variable sand properties(initial conditions) and entry/exit points of matters and energy (boundary conditions) need
to be specified.
Computer world - The computer world consists of the numerical solvers of themathematical world and an interface with the user to communicate parts of theconceptual/mathematical world to the numerical world, at minimum via input files
sometimes generated via a graphical user interface. Via output, the user obtains informationfrom the computer world for interpretation or comparison with the physical world.
Depending on the specific tool, part of the mathematical world is fixed within a source code,
whereas another part can be defined via the input files by the users.

Figure 1 – General flowchart for the modelling environmental-

Coupled Reactive Transport CodeA coupled reactive transport model integrates at minimum the numerical implementation of two sets ofprocess-based equations for:••Transport of aqueous species
o Focus in course on multi-phase continuum physics
Interaction between components (chemical elements) – thermodynamic equilibriumo Focus in course on low-temperature, low-pressure multi-component geochemistry
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

2
ISC: PublicIntroduction

Figure 2 - Coupling of different processes - bold text indicates the minimum requirement to have a coupled
reactive transport model

Pore-scale model and continuum representation
A key point in conceptual and mathematical models is the representation of heterogeneity. In thiscourse, the focus is on porous media (soil, aquifers, construction materials) which consists of aheterogeneous structure of void space (pores) and solid phases with sharp boundaries. A reactive
transport model representing the microstructure explicitly is often referred to as coupled reactive porescale model (e.g., Patel et al. (2014), Oostrom et al. (2014)). Alternatively, a continuum-based approach
is followed in which the heterogeneous structure is mathematically homogenized over a given volume
such that physical and chemical properties vary continuously in space. Rather than the existence of only
one phase at a given location in the simulation domain (as is the case in pore scale models), solids andfluids coexist simultaneously at each point in space.
Figure 3 shows a continuum representation of a three-phase porous system. Rather the representing the
distribution of solid, aqueous and gaseous phases explicitly, the phases are quantified in terms ofvolumetric content. Physical processes (e.g. mass entering or leaving the system, represented by the blue
arrows at the left hand side) may trigger a redistribution of chemical components between the different
phases including sorption on surfaces of the solid phases by a set of chemical equilibrium reactions. At
the same time, redistribution can be a continuous process with time when kinetic reactions are
happening. In some systems, the redistribution of species can also be controlled by a buffer reactions
via exchange with a solid or gas phase buffer (boxes at the right hand side).

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

3
ISC: PublicIntroduction

Figure 3 – Example of a continuum representation of a three-phase porous system. Chemical components will be
distributed between these phases according to equilibrium (double arrows) and kinetic (single arrows) reactions.
The system will react when mass enters or leave the system (blue arrows at the left hand side) and might bebuffered by a gas or solid phase.

Generic Coupled Reactive Transport Codes
Multi-process
Modern coupled reactive transport codes (Steefel et al. 2015, Steefel et al. 2005) offer much more
possibilities than the coupled processes of solute transport with geochemistry. They combine processes
that are relevant for understanding and assessing natural and engineered systems and their evolution.
Without being complete, following main processes are included in some/most modern coupled reactive
transport codes (Figure 2):•••

Physicalo Water flow – steady-state saturated water flow, variable-saturated flow (Richards
equation), multiphase flow, density-drive flow, non-isothermal flowo Heat transporto Transport of colloids, electro-migration (Nernst-Planck equation)o Mechanical processes
Chemistryo Geochemical equilibriumo Kinetic reaction networks
Biologicalo Microbiologically mediated reactions
o Biomass growtho Effects of roots and vegetation

Many of these codes allows for different feedback mechanisms between properties and parameters forthe different processes. Temperature has an effect of thermodynamic and rate parameters. Precipitationand dissolution processes of solid phases (and biomass) affect different flow and transport parameters.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

4
ISC: PublicIntroduction
Problem solving environment
Besides including different processes, modern coupled reactive transport models are also generic innature. Generality refers to the ratio between the part of the conceptual/mathematical model that isincluded in the source code and the part that users can adapt via (text) files (database – input in Figure1). Specific reactive transport codes (e.g. UNSATCHEM, Šimůnek and Suarez (1994)) includes a significant
part of the model in the source code (such as components, reaction equations and network) and users
can adapt only specific properties (such as initial amount) – the model formulation is fixed at compilation
time. Generic reactive transport codes allows users to define the reaction network, implement different
rate equations, choose between different types of models (transport models, activity correction models,surface complexation model, etc.) – the compiled code allows for solving different model formulations.

Figure 4 - Concept of a specific reactive transport code (left – conventional model code) and a generic reactive
transport code (right – modelling environment) (after Meysman et al. (2003a).

Meysman et al. (2003a) introduced the term problem solving environment (PSE) to specify a genericcoupled reactive transport model. Two important key words are:••Flexibility – the user can control the specific model formulation starting from a set (toolbox)of available model components (geochemical entities as species, surfaces, minerals andprocesses as equilibrium, kinetics, uniform or heterogeneous flow and transport, definitionof properties, etc.) without the need to write a line of code that requires compilation.
Extensibility – the user can extend the modelling environment with new components (suchas new species, reactions, kinetic laws, relations between porosity and transport properties).

Modern coupled reactive transport codes embed these two elements although the level and theimplementation of flexibility and extensibility may vary. Together with the hard-coded modelcomponents, flexibility and extensibility will determine the choice of a coupled reactive transport modelby the user.

Environmental and Engineered Applications
Scope of reactive transport models
Coupling geochemical models with flow and transport solvers allows for an explicit representation of
the spatial-temporal evolution of multiple geochemical species and solid phases. The capabilities of the
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

5
ISC: PublicIntroduction
transport solver mostly determines the spatial scope of the reactive transport models that ranges fromµm-mm to several hundreds of kilometres (Figure 5).
Models developed for simulating one-dimensional flow and transport processes have typically a spatial
range of mm to several tens of meters in horizontal, vertical or inclined direction. The possibility toinclude heterogeneity within the domain is evidently limited to one direction with homogeneous
properties and gradients in the other directions and bounded by impermeable boundaries in the
direction parallel to the flow and transport direction. Such reactive transport models are typically used
for laboratory applications both in the field for environmental as for engineering applications. In
environmental applications, this scale is often referred to as the field-scale especially for processes
related to near-surface environments (soils) with complex water redistribution processes.
Two- and three-dimensional forms of the water flow, heat transport and solute transport equations
broaden the applications to cases in which heterogeneity and gradients may exists in different scales.
Spatial range varies again from mm and can extend to hundreds to thousands of meters. Although not
limited to the two or three dimensional flow and transport solvers, one of the typical environmental
applications is the study of the critical zone, the zone from the outer extent of vegetation down to thelower limits of groundwater (Brantley et al. 2007). Catchment or watershed models extending up to
several tens to hundreds of kilometres includes additional processes such as subsurface-surface waterinteractions which coupled with geochemical solver opens several additional opportunities (Bao et al.
2017, Li et al. 2017, Prommer et al. 2003).

Figure 5 – Spatial scope of reactive transport models for environmental applications: Watershed- catchment scale(km) (left), macro-scale / profiles (m – mm, center), and micro-scale /pore scale (µm). Fromhttps://www.juliaperdrial.com/research-projects.html after Bao et al. (2017)

Applications
Many studies exists for different applications within environmental, engineered or combined
environmental-engineered systems such as in following non-exhaustive list:•••

Environmental applications
o Diagenesis (e.g. Meysman et al. (2003b)), weathering and soil development (e.g., Finke
and Hutson (2008), Keyvanshokouhi et al. (2016))
Engineered applicationso Ageing of construction materials such as concrete
Environmental-engineered applicationso Disposal of nuclear waste: Bildstein et al. (2019)

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

6
ISC: PublicIntroduction

Course software
HPx
This course will use the HPx reactive transport code (Jacques et al. 2018, Jacques et al. 2008a, b). HPx(Figure 6) couples the HYDRUS software (Simunek et al. 2013, Šimůnek et al. 2016, Simunek et al. 2008)(see https://www.pc-progress.com/en/default.aspx) with the geochemical solver PHREEQC (Parkhurst
and Appelo 1999, 2013a) (see https://www.usgs.gov/software/phreeqc-version-3).
HYDRUS software (HYDRUS1D, HYDRUS2D/3D) solves the water flow equation for variably saturated
porous media, heat transport, root water and solute uptake, advective-dispersive transport process inthe aqueous phase, and diffusive transport processes in the gaseous phase.
PHREEQC is a powerful geochemical code for aqueous speciation, mineral equilibria, multi-site cationexchange, complex surface complexation modelling, exchange with a gas-phase, solid solutions andkinetic reactions. It gives a large flexibility for modelling different geochemical systems and experimental

set-ups including titration, speciation plots, etc. Moreover, PHREEQC has the possibility to perform one-dimensional advective-dispersive transport under steady-state flow conditions and multi-dimensional
diffusive transport..

Figure 6 – HPx coupling HYDRUS software with PHREEQC.

As a result, HPx solves flow and transport with biogeochemical processes in an integrated way () allowing
for exchange of information on state variables between the two codes (Figure 7).
HPx includes a slightly modified version of PHREEQC 3.6.1 offering:
•

•
•
•

An alternative scripting language, MyBasic, which is a light weight BASIC interpreter withstructured blocks, object-oriented possibilities etc. (see https://github.com/paladin-t/my_basic)
including all standard PHREEQC basic functions
Use of inline variables and some linked HP specific keywords
Other keyword for plotting
Enhanced functionalities of scripting with CALCULATED_VALUES, especially related to reactive
transport problems

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

7
ISC: PublicIntroduction

Figure 7 – Processes solved by respectively the HYDRUS software (left) and PHREEQC (right) and the exchange of
information of state variables between them.

For one-dimensional transport, the current version is HPx 3.01.003, included in HYDRUS-1D 4.18.007.
This version is used during the course. Beside the latest version of HPx, the HYDRUS-1A 4.18.007
graphical user interface includes new dialogue form for pre-processing the geochemical input and postprocessing of geochemical output variables. We refer to the combination of these two dialogue forms
as HPGUI. The current version is 01.01.
For two and three dimensional transport (Šimůnek et al. 2013), an older version of HPx (version 2.02) isstill included in HYDRUS2D/3D graphical user interface. An update is foreseen end of 2020 or beginning2021.
The same PHREEQC library as in HPx is also used in combination with MT3D-USGS (Bedekar et al. 2016)
which is under development as MFTHP (Leterme et al. 2019) . MFTHP is a framework in which coupledflow and transport between unsaturated (soil) and saturated (aquifers) systems are coupled. Modflow(MF) is used for water flow in the aquifer system (and surface water system), MT3D-USGS (T) is used
for solute transport in the aquifer system (and surface water system), HYDRUS-1D (H) is used for flow
and transport in the soil system, and PHREEQC (P) is used for geochemical calculations. The HYDRUS
package for MODFLOW (Seo et al. 2007) is adapted and for solute transport, algorithms similar to thosedescribed in(Beegum et al. 2019) are used.
PHREEQC and its graphical user interfaces
PHREEQC is used for geochemical modelling during the course. The current version is 3.6.2(https://www.pc-progress.com/en/default.aspx). Although PHREEQC input can be made with standard
ASCI text editors and the program can be run from command prompt, we use a graphical user interface
that combines the pre- and post-processing for HPx from the HYDRUS graphical interface into
HPGeochemistry. HPGeochemistry embeds the same PHREEQC version as HPx (modified PHREEQC 3.6.1)
and is available as a separate software (HPGeochemistry,) installed during installation of HYDRUS1D (4.18.007). The PHREEQC input in HPGeochemistry is compatible with HPx input, and some aspects
are shared between HPGeochemistry and HPx (e.g. database management).
For completeness, we list here some alternative versions and graphical user interfaces for PHREEQC. The
most common graphical user interface for PHREEQC is PHREEQCi (see https://www.pcprogress.com/en/default.aspx). Alternative interfaces are the notepad++ interface from C.A.J. Appelo(see . https://www.hydrochemistry.eu/downl.html). For repetitive PHREEQC calculations, plotting andinversion, the PHREEPLOT version can be useful (Kinniburgh and Cooper 2009) (seehttp://www.phreeplot.org/). Recently, an alternative PHREEQC version with graphical user interface andadvanced project and data management options was launched (GibbsStudio, https://gibbsstudio.io/).
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

8
ISC: PublicIntroduction
Note that PHREEQC is also used in a number of other reactive transport models, mainly because of the
PHREEQC modules (iPHREEQC (Charlton and Parkhurst 2011) , and the reaction module PHREEQCRM(Parkhurst and Wissmeier 2015)). Getting familiar with PHREEQC opens the road to select from a numberof different code for coupled reactive transport modelling. Examples are:
•

•
•
•

YANTRA (Patel et al. 2014, Patel et al. 2018) a pore-scale modelling software that uses a lattice
Boltzmann approach for the transport processes (Navier-Stokes, diffusion, advection-diffusion).
iPHREEQC is used as one of the geochemical solvers. The spatial scale is typically small (100µm3).
Coupling with COMSOL (Nardi et al. 2014, Perko et al. 2015, Wissmeier and Barry 2011)
PHT3D (Prommer et al. 2003) coupling with MT3DMS
Coupling with OPENGEOSYS (He et al. 2015, Jang et al. 2018)

Further developments
Developments are ongoing for:
•
•
•
•
•

Embedding python
Alternative geochemical solvers (replacing PHREEQC), either by script functions or using
ORCHESTRA (Meeussen 2003); testing surrogate model approaches
Linux version for HP1
Further development of the graphical user interface HPGeochemistry and HPGUIImprovement of manual, help file, benchmarks and examples

Detailed course outline
During the course, examples are given for different types of models (Figure 8):
Speciation models: Given information on temperature, pressure, and composition (moles of
elements), a single calculation will give us the distribution of the elements over different phases,
the speciation, solid phase composition, saturation indices, etc. We can refer to it as closed, static,
batch or beaker-type of system.
Reaction path models: Starting from a give state of a system (determined by its temperature,pressure and composition), we analyse how the system changes when one or more of thesevariables changes. Composition or any other geochemical variable of interest is calculated as afunction of reaction progress (ξ) expressing e.g. the amount of a reactant that is added to thesystem, time (e.g. kinetic processes), temperature, pressure or some buffered values (e.g. sliding
activity or fugacity studies).
Reactive transport model (RTM): The distribution of elements between different phases (timedependent biogeochemical processes and reactions) is coupled to the movement (transport) of
solutes by advection and hydrodynamic dispersion processes. In essence, the equation for solutetransport (advection-dispersion equation) is coupled with equation to solve for thermodynamic
equilibrium or kinetic processes. Other processes can be integrated as well such as water/fluid
flow or heat transport.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

9
ISC: PublicIntroduction

Figure 8 – Schematic representation of different models based on thermodynamic and kinetic geochemical
modelling (from Zhu and Anderson (2002))

Session 1 – Aqueous speciation modelling
Based on information on the composition of ground or sea water, we calculate the speciation and thesaturation index of different minerals. A short introduction is given on geochemical systems and thebackground of thermodynamic equilibrium and its mathematical formulation. The graphical userinterface for HPGeochemistry is introduced that uses a modified version of PHREEQC 3.6.2 as thegeochemical solver. Aspects of project and database management are demonstrated.
Next, the PHREEQC keywords for aqueous speciation calculation (SOLUTION_MASTER_SPECIES,
SOLUTION_SPECIES, PHASES, SOLUTION) are introduced and applied in an exercise to calculate thespeciation of sea and ground water. The initial calculation of a solution composition is completed byshowing how one can charge balance a solution or fix the saturation index of a mineral.
Session 2 – Reaction path modelling
Several implementations for reaction path models are demonstrated. A first exercise shows how one canmix two solutions (MIX). In addition, two specific reaction paths are applied: a sliding concentration
reaction path (titration) applied to speciation of Cd when adding Cl, and a sliding activity reaction pathapplied to the solubility of Fe(OH)3 as a function of pH. We illustrate how structured input files can bemade, repetitive PHREEQC calculations are optimized and how HPGeochemistry can be used to visualize
output results.
Several PHREEQC keywords are introduced: geochemical entity (EQUILIBRIUM_PHASES), reaction pathmodelling (REACTION), output (SELECTED_OUTPUT, USER_PUNCH, SELECTED_GRAPH) and some
HPGeochemistry specific keywords for repetitive calculations (GLOBAL_VARIABLE, COUNTER,
COUNTER_INCR, REPEAT).
Session 3 – Reactive transport modelling
Several examples from literature illustrate the power of reactive transport models. After a short
introduction on the mathematical equations of coupled reactive transport, the reactive transport code
HPx is introduced. We implement three examples of a reactive transport via the HYDRUS-1D graphicaluser interface. The examples includes equilibrium mineral dissolution and precipitation in systems withcalcite, gypsum and/or otavite.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

10
ISC: PublicIntroduction
Session 4 – Equilibria at the solid-liquid interface
Several geochemical processes control sorption at a solid-liquid interface. Only cation exchange at asingle cation exchange site is introduced during the course using the PHREEQC keywords
EXCHANGE_MASTER_SPECIES, EXCHANGE_SPECIES, and EXCHANGE. First, an exercise to calculate thespeciation of an exchanger in equilibrium with a solution is used to illustrate the keywords USE and
SAVE. Two reaction path models including cation exchange are implemented: cation isotherms andevaporation of a soil solution. Finally, the process of cation exchange is integrated in a reaction transport
problem (major cations, one or more pulses).
We also introduce the concept of a multi-dataset file which is illustrated by calculating and plotting
isotherms at different ionic strengths. We also illustrate how a time series plot with different geochemicalvariables at a given observation point can be defined for the reactive transport problem.
Session 5 – Kinetics
For many biogeochemical processes, kinetic reaction are rather the rule than equilibrium reactions. We
introduce the power of PHREEQC for implementing user-defined rate equations via scripting using the
keyword RATES and including these rate equation in the geochemical model via KINETICS. These two
keywords are illustrated to model (i) kinetic dissolution of minerals using a rate equation based ontransition-state theory, (ii) kinetic biodegradation of phenol using a Monod type of rate equation, and(iii) a sequential-parallel first-order degradation network (PCE and daughter products). The latterexample is then included in a reactive transport problem for three conditions: (i) boundary value
problem, (ii) initial value problem, and (iii) initial value problem for a mobile-immobile region transport
problem.
In addition, we demonstrate the capabilities of the user interface for HPGeochemistry to define different
cases within when project simulation. In the initial value problem examples, the option to obtain summed
values for horizons or materials is used.

About the notes
Text font and colouring
In the notes, text related to input or output of the numerical codes does have typically a different fontand colouring. In general, PHREEQC input names are in Courier New and coloured according to theirfunction. The colours in the text are similar to the colour used in the graphical user interfaces. Table 1
gives an overview of the text font and colour conventions used throughout these lecture notes.
Table 1 – Text font and colouring of PHREEQC related input as used in the lecture notes and the graphical userinterface.

KEYWORD
Identifier
Text
Command
Script function
MyBasic standard
function
Inline Variable
NumberQuoted text
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

11
ISC: PublicIntroduction
Thermodynamic name
Text within angle brackets <: > with acolon to indicate the syntax of the
input (see further)

Syntax

Instructions for the hands-on training examples are indicated with a vertical bar:
Go to Geochemistry tab
Add SOLUTION 1001

Code examples
These notes contain several examples for PHREEQC related input. Typically, these codes are printed ona grey background.
When the code examples are discussed in the text, they have a caption number (Code) and line numbers
for reference during the discussion (Code 1).

SOLUTION 1
pH 7 charge
Na 1

Code 1 - Code example

When part of a hands-on training example, they typically contain no caption number nor line numbers:
Go to Geochemistry tab
Add SOLUTION 1001
SOLUTION 1
pH 7 charge
Na 1
Add SOLUTION 1002

Syntax examples
Throughout these lecture notes, several input commands for PHREEQC and Gnuplot are introduced.
Beside given specific examples, the syntax of these commands is given in which user-defined variables
and parameters are indicated by angle brackets <> that contain the type of the variable (e.g. integer,
floating point, …) and the name:<:int:Number>
A full list of conventions is given in Table 2.
Table 2 - Conventions for syntax

intInteger number

flt

Floating-point number

txt

A alphanumerical variable

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

12
ISC: PublicIntroduction
str

A alphanumerical variable between
double quotes “”[]Optional input|

Separate alternatives{}

Can be repeated more than one
time

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

13
ISC: PublicIntroduction

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

14
ISC: PublicSession 1 – Aqueous speciation calculations

Session 1

Learning outcomes
Understand the major principles of thermodynamic equilibrium modelling
Familiarize with the graphical user interface HPGeochemistry
Use PHREEQC for aqueous speciation problems

Geochemical System
Geochemical system
A geochemical system consists of different phases with different chemical species that are formed ofcomponents.
A phase consists of one or more chemical species within a separate physical entity and withhomogeneous properties. The most relevant phases for systems within the scope of this course are:•••Liquid phaseso A single aqueous phase (Leal et al. 2017) – indicated as (aq) - with the H2O as thesolvent and several other solutes, i.e. aqueous dissolved species.
o Some other non-aqueous liquid phases may exist
Gaseous phase indicated as (g) – only a single one in a system because of mixing capabilities(Leal et al. 2017)
Solid phases indicated as (s) – present as pure minerals or solid solutions. In a pure mineral,the phase consists of a single species, whereas solid solutions consists of two or more
species in variable amounts.

The species are composed of components; the latter being chemical elements or a linear combinationsof them. Components are also called primary of master species. For example, a system with components
O2-, H+, Ca+2 and C+4 can contain :•••aqueous species as H2O (liquid water), H+, OH-, Ca(OH)+,
gaseous species as H2O(g), CO2(g),solid phases as different CaCO3 polymorphs calcite, aragonite 1

As discussed in Leal et al. (2017), a unique name is required for each species by specifying its chemical
formula, its physical state, and the name of the phase. Throughout the text, we uses following
conventions:•••Aqueous species are identified just by their chemical formula – in case of possible confusion,
the subscript (aq) will be added (except (l) for H2O). As only one aqueous phase exists, thephysical state (aq) and the name of the phase (aqueous) are the same.
Gaseous species are identified by their chemical formula and the subscript (g). (g) indicatesthe physical state. Again, as only one gaseous phase exists, the physical state (g) and phasename (gaseous) are the same.
Solid species are indicated by the phase name (e.g. calcite, aragonite, quartz). In case thename is a chemical formula, the subscript (s) is added. We prefer to indicate solid phases by
their phase name rather than by their chemical formula as the same chemical formula is used
for polymorphs of a mineral (see CaCO3 example above).

1 In a geochemical system in equilibrium, only one CaCO

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

3 polymorph exists.

16
ISC: PublicAqueous speciation calculations

Figure 9 Representation of a geochemical system. Left: Beaker model (Deutsch and Siegel 1997).

Key question – Chemical state
When the temperature T, the pressure P and composition (the amount of master species, in moles) of ageochemical system are defined, the main interest is to know the chemical state of the system, i.e. themolar amounts of the chemical species in the system for the given temperature and pressure. The
relevance of knowing a chemical state is when one wants to know derived variables such as pH, amountof a given solid phase, saturation state, total amount of a master species in a given phase, mole fractions,concentrations, etc. Chemical state calculations needs to be done:•••

To know the equilibrium chemical state of a system
To know the evolution of the chemical state of a system when one or more reactions aretime-dependent (e.g. a phase transition from solid to aqueous by dissolution)
To know how the chemical state of a system changes when T, P or composition changes (e.g.,titration, during transport processes or evaporation)

ExamplesCa-C-O-H system
An overview of the species in a geochemical system with three states of matter (liquid, gas, solid) andmaster species Ca, C, O and H is given in Table 3. The species are selected from supcrt98 database which
is an updated version of the supcrt92.data database (Johnson et al. 1992). 2 There are two phases in thesolid state – calcite and aragonite – both with the same chemical formula of the species. However,species with the same chemical formula but in a different phase have distinct thermodynamic andchemical properties.
Table 3 – Species in a geochemical system with Ca, C, O and H

State
LiquidPhasesAqueous

Gas

GaseousSpecies
H2O(l)
OHH+
CO2(aq)HCO3H2CO32Ca2+
CaOH+
CaHCO3CaOH+
H2O(g)

2 An xml version of the supcrt98 database is provided in the complementary course material folder. This xml version

was obtained from https://reaktoro.org/thermodynamic-databases.html.
PHREEQC-HPx Course, Prague, September 3-4, 2020

SCK CEN/39560003

17
ISC: PublicSession 1
Solid

CalciteAragonite

CO2(g)CaCO3
CaCO3

Mercury in soils
This example illustrates how diverse a chemical element can occur in different phases in a multiphasesystem. When assessing e.g. the risk for groundwater pollution, it is important to understand thespeciation and quantify their distribution. Mercury levels, Hg, presence under three valence states andpossible many different species (schematic model shown in Figure 10), in soils are legacy of different
industrial processes over several centuries. The original form of the Hg pollution could originate as solid
cinnabar (HgS), residual non-aqueous liquid mercury (Hg0, not present as a continuous liquid phase, but
as a residual phase trIn apped in pore spaces) or as mercury chloride. In the aqueous phase, both Hg0
and HgII are present, the latter forms complexes with inorganic and organic (dissolved organic matter)
species. Methylmercury is an organic form of mercury that is more abundant under anoxic conditions.
Hg0 is volatile and occurs in the gaseous phase as well. Mercury binds to solid organic matter, especially
via thiol (SOH) functional groups.

Figure 10 – Scheme of possible speciation and reactions of Hg in a soil system (Leterme et al. 2014).

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

18
ISC: PublicAqueous speciation calculations

Aqueous Speciation
Problem Definition
Given the composition, temperature and pressure of an aqueous phase, calculate the equilibrium
speciation.

Conceptual Model
The problem involves following components:•••

Equilibrium reactions
One phases (aqueous)
Well-mixed

We model an aqueous speciation problem.
Figure 11 shows a schematic representation of the key elements to calculate aqueous speciation withthe most important PHREEQC keyword SOLUTION.

Figure 11 – Schematic representation of the model to calculate aqueous speciation (left) and key elements of thecomputer model.

Mathematical world – Aqueous speciation
Principles of thermodynamic equilibrium
Major principles
The goal of solving an geochemical equilibrium problem is to find the most stable state of a system
for a given set of pressure, temperature and compositional constraints.
The major principles of modelling (biogeo)chemical equilibrium processes are(i) conservation of mass, and(ii) minimum of free energy in the chemical system.
The first principle is given by the law of conservation of mass resulting in a set of mass/moles balance
equation for components. The second principle will be achieved by minimizing the free energies of thereacting species or by defining a set of mass law equations for all chemical reactions involved in the
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

19
ISC: PublicSession 1
system. The mass law equation of a chemical reaction states that the ratio of the product of activities of
the reaction products over the product of activities of the reaction reactants is constant.

Chemical reactions and equations
A chemical reaction is a process of a chemical change of the reactants, i.e., the substances initiallyinvolved, to products. The changes are typically the motion of electrons in the forming and breaking ofchemical bonds. Besides homogeneous reactions as acid-base, aqueous complexation and redox
reactions, precipitation/dissolution and sorption reactions on minerals or colloids are important
heterogeneous reactions in geochemical systems.
The symbolic representation of a chemical reaction is the chemical equation. A chemical equation should
follow the law of conservation of mass stating that the total mass of all products of a chemical reaction
is equal to the total mass of all reactants of that reaction. Therefore, a chemical equation is also calledmole-balance equation, especially in aqueous chemistry. In addition, a chemical reactions must also
conserve charge.

Mole-balance equation
A general mole-balance equation in a closed system involving N species at most is written as:No/ = ∑ν i Ai

(1)

i =1

where o/ refers to the null species, vi are the stoichiometric coefficients of the reaction, and Ai the involved
species. If the stoichiometric coefficient of a species is positive, the species is called a product, otherwise
it is called the reactant.
The same holds in a closed system in which NR simultaneous chemical reactions take place:

o/
=

N

( r 1,..., N )
ν A
=
∑
i =1

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

ir

ir

R

(2)

20
ISC: PublicAqueous speciation calculations

Conservation of mass

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

21
ISC: PublicSession 1

Thermodynamic equilibrium

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

22
ISC: PublicAqueous speciation calculations

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

23
ISC: PublicSession 1

Activity correction models for aqueous species

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

24
ISC: PublicAqueous speciation calculations

Saturation state and saturation index

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

25
ISC: PublicSession 1

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

26
ISC: PublicAqueous speciation calculations

Computer world - Interface – HPGeoChemistry
Learning Outcomes
Explore the HPGeochemistry graphical user interface
Set up the course project
Define the thermodynamic database
Define a simulation case and explore the input part

HPGeochemistry graphical user interface
Start HPGeochemistry to obtain the start screen (Figure 12).
The toolbar on the start screen is divided in three parts:•

Access to managerso••Project Manager to open/create/edit cases and projectsoThermodynamic database manager
Access to different workspaces and actionsoAccess to input workspaceoSave the projecto

Execute the projectoAccess to the phreeqc.out fileoAccess to the plotting workspaceoAccess to the output file workspace
Access to the preferences, settings and help filesooAccess to preferences and settings

Access to different help files
 HPx manual – instruction on the graphical user interfaces of HPx
specifically for geochemical input and output ; Overview of input
and output commands for the geochemical calculations (e.g. data
blocks of PHREEQC, specifically the commands not available instandard PHREEQC) and scripting ; How to instructions ;verification and example projects
 Phreeqc (USGS) – the official PHREEQC manual distributed withthe phreeqc program Phreeqc (T. Appelo) – a phreeqc help file created by C.A.J. Appelo
 Gnuplot – the help file of gnuplot

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

27
ISC: PublicSession 1

Figure 12 – Main form of HP Geochemistry

Database Manager
Manage databases
The thermodynamic database manager allows for defining shortcuts to different databases.
Thermodynamic databases are grouped in database group (e.g., a group can contain different versionsof a database).
Click Database Manager

buttonFigure 13 – Database Manager
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

28
ISC: PublicAqueous speciation calculations
DataBaseGroup:
•
•

•
•

To manage different database groups
A database group can contain different databases or different version of a database. A projectis always linked to a specific version of a database which improves traceability of yourgeochemical calculations
By default after HYDRUS installation, different versions of databases from different PHREEQC
versions are included – in Figure 13, this is only the latest version 3.6.1.
The side toolbar allows to:oCreate a new database groupo

Edit the description of the database group (not allowed for the database groups
installed by HYDRUS)oRemove the database group from the list (files are not deleted; not allowed for the
database groups installed by HYDRUS)

Items:•••

To add or remove databases
For the database groups available after HYDRUS installation, the databases from PHREEQC are
already added.
The side toolbar allows to:o

Add an existing databaseo

Edit the description of the database (not allowed for the database groups installedby HYDRUS)oRemove the database from the list (files are not deleted; not allowed for the databasegroups installed by HYDRUS)

Add Database Group and DatabaseIn this section, we illustrate how a database group and a database is added to the HPGeochemistrycodes. We add the database Thermoddem from the French National Geological Survey and add it to a
separate folder in the default folder for the databases.
To find the default folder for the databases and create a new folder:
Click the Settings and Preferences

buttonSelect Settings tab
Select and copy the text of the folder for the databases
Create the folder Thermoddem in that folder using a file explorer.
Next, we download the database from the internet.
Go to the website https://thermoddem.brgm.fr/
Click “Download Database”
Select ”Phreeqc”
Click “Download database”
Save and/or open the database in the Thermoddem folder of the default database folder.
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

29
ISC: PublicSession 1

Finally, add the database group Thermoddem and the current version of the database.
Click the Database Manager

buttonClick the Create a new Database Groupwindow for a new database group opens

button – The Meta Information Dialogue

Select Default Folder
Type the name of the database group in the Name textbox (Thermoddem)
Click the OK button
Click the Open Existing DatabasebuttonSelect the database file (if you save the file in the correct folder, this file will be there)
The Meta Information Dialogue window for a new database opens
Click OK button
We are now able to use this database in our projects.

Project Manager
Manage Projects
The Project Manager allows for defining simulation cases (add existing cases, copy an existing case or
create a case), accessing and editing meta-data to cases, remove cases from the project manager ordelete cases. Simulation cases can be grouped in projects.
Click Project Manager

buttonD. Jacques (djacques@sckcen.be)
SCK CEN/39560003

30
ISC: PublicAqueous speciation calculations
Figure 14 – Project Manager

Project:•••

To manage different projects
A project can contain different simulation cases
The side toolbar allows to:oCreate a new projecto

Edit the description of the selected projectoRemove project (files are not deleted)

Items:••

To add or remove cases
The side toolbar allows to:o

Add an existing simulationo

Copy an existing simulationoCreate a new simulationo

Edit the description of a simulationoRemove the simulation from the list (files are not deleted)oDelete the simulation from the list (files are deleted as well

Create the Course Project
Click Create new project

buttonGive Name Course HP 2020
Click Browse button
Select root folder for the project. In that folder, the different simulation cases are added as
separate folders.

Computer world – PHREEQC Input
Implementing the model
Figure 15 shows the PHREEQC keywords for implementing an aqueous speciation problem. The
keywords SOLUTION_MASTER_SPECIES, SOLUTION_SPECIES, PHASES and SOLUTION are
introduced in this section.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

31
ISC: PublicSession 1

Figure 15 - Implementation of the model for aqueous speciation

Phreeqc conventions
•

•

•
•
•
•
•

Simulations within an single input file are separated by the keyword END. A simulation consistsof a number of entities that are delimited with key words (e.g. SOLUTION). The order of the
keywords in a simulation is not important for the standard PHREEQC keywords, but can beimportant for some of the HPGeochemistry specific keywords (such as GLOBAL_VARIABLE, seelater). Most data blocks consists of a number of identifiers (e.g. –water) which order with a data
block is also not crucial. However, in most cases, the information given after the identifier isorder specific
In most cases, the input is case insensitive; the most important exception are the element names,
see later. Some definitions for the interaction with the output work space is also case sensitive(see later).
Comments start with #; the remaining of the line is neglected.
The semi-colon ; is interpreted as a line separator. Different input lines can be combined in one
physical line.
The backslash \ indicates line continuation. The input lines continues on the next physical line
@ indicates an inline value – the text will be replaced by a previously defined value (see later).
A hyphen - between integers indicates a range of integers.

Code 2 illustrates some of these conventions.:
•
•
•
•
•
•

Line 12 – Separation of two simulations
Line 1 and 2 – Commented lines
Line 4, 5, 8 , 9 – Comments within a line
Line 6 – Separation of input lines by ;Line 14-15 : Input line continuation by \Line 18 : Range of integers

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

32
ISC: PublicAqueous speciation calculations#Commented line
#First simulation
GLOBAL_VARIABLE CEC 0.25 "eq/kgw"
SOLUTION 1 #Begin of data block-pH 7 #identifier within a data block-pe 4 ; -temperature 25 #Separated lines
EXCHANGE 1 #Second data block
X @CEC #inline variable-equilibrate with solution 1
END #separate two simulations
SELECTED_OUTPUT 1-molalities Cd+2 CdCl+ CdCl2 CdCl3- CdOH+ \Cd(OH)2 Cd(OH)3- Cd(OH)4-2 CdOHCl
END
SOLUTION 2-6
Code 2 – Example of conventions in PHREEQC / HPGeochemistry

A input file can contain several simulations separated by the keyword END. For example, Code
3 contains 4 simulations. Each simulation contains a number of entities. During the simulations,first, for some entities (e.g. SOLUTION, EXCHANGE) an individual initial calculation will be
performed followed by a single (e.g. simulation 1) or multiple (e.g. simulation 3) batch reaction
simulations for all entities. If requested, calculations results are saved into memory (SAVE, seelater) or written to output files (see later). Entities can be defined:• Explicitly by the corresponding keyword (e.g. simulations 1 and 2); or• Implicitly by the keyword USE followed by the chemical entity, e.g. simulations 3 and 4).

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

33
ISC: PublicSession 1#Simulation 1
SOLUTION 1-pH 7
Na 0.1
Cl 1 charge
EQUILIBRIUM_PHASES 1
CalciteEND#Simulation 2
SOLUTION 2-pH 9 charge
Ca 0.1C 1 calcite 0
EXCHANGE 2
X 0.1-equilibrate with solution 2
END#Simulation 3
USE exchange 2
SOLUTION 3-pH 7 charge
Mg 0.2
REACTION 1
Ca 1
0.1 in 100 steps
END#Simulation 4
USE solution 1
USE reaction 1
EQUILIBRIUM_PHASES 4
Dolomite
END
Code 3 – Example of conventions in PHREEQC / HPGeochemistry

SOLUTION_MASTER_SPECIESSOLUTION_MASTER_SPECIES defines the correspondence between element names and aqueous
primary and secondary master species . The syntax and an example are given in Code 4. Line 5 and 6
define a primary master species, line 7 and 8 define a secondary master species.
•

•
•
•

•

Name: The name of the element. For redox sensitive elements, the valence state is added as well(e.g. S, S(6) and S(-2) in line 6-8. Element names are used to define the concentration in the
aqueous phase (SOLUTION) or to identify for which elements total concentrations are requiredin the output.
Master species: Formula of the primary (e.g. for Ca in line 5) or secondary (e.g. line 7-8)
master species. Used in the reactions of the solution species (SOLUTION_SPECIES)
Alkalinity: the alkalinity contribution of an element (or element valence state) which isdefined as the number of protons the elements consumes upon titration to pH 4.5.
Gwf|formula : The gram formula weight (gwf) is used when the input concentration is definedin mass/L or mass/kgw (e.g. g/L) to convert it to mole/kgw. In case of a formula, the formula isused to calculate the gram formula weight.
Gwf elem: The gram formula weight of the element (only for primary master species).

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

34
ISC: PublicAqueous speciation calculations
SOLUTION_MASTER_SPECIES<:txt:Name><:txt:masterspecies><:flt:alkalinity><:flt:gwf>|<:txt:formula>[<flt:gfw_el>]
SOLUTION_MASTER_SPECIESCa Ca+2 0 Ca 40.08
S SO4-2 0 SO4 32.064
S(6) SO4-2 0 SO4
S(-2) HS- 1 S
Code 4 – Syntax and example of SOLUTION_MASTER_SPECIES

More information on the keyword SOLUTION_MASTER_SPECIES is available on the online manual
page.

SOLUTION_SPECIESSOLUTION_SPECIES define the chemical reaction and its corresponding thermodynamic(thermodynamic constant, activity correction model and parameters) and supplementary (e.g. specificconductance, specific volume, …) information.
Code 5 shows some identifiers in SOLUTION_SPECIES; Figure 16 shows an extract from the phreeqc.dat
database with some examples:•••Line 3 : reaction: The association reaction. For each master species, the identity reaction
equation needs to be defined (e.g. Ca+2 = Ca+2). For secondary aqueous species, the secondary
species is the first species to the right of the equation.
Line 4: log_k: The log of the equilibrium constant at 25°C
Line 5: delta_h: Enthalpy of the reaction at 25°C – is used in the van ‘Hoff equation to calculate
the temperature dependency of the equilibrium constant: K T  ∆H  1 1 
− 
ln=



R  Tref T  K ref (3)<here KT and Kref are the equilibrium constants at the required temperature T [K], and thereference temperature (25°C), ∆H is the enthalpy of the reaction and R is the universal gas
constant.
•

Line 6: -analytic: indicates the six parameters of an analytical expression for the temperaturedependency (T is temperature in Kelvin):

Log ( K ) = A1 + A2T +•A3
A+ A4 log(T ) + 52 + A6T 2
T
T

(4)

Line 7: -gamma: indicates the activity correction model.
o If –gamma (and –llnl_gamma 3)is not present, the Davies equation is used to calculate
the activity:µlog ( γ ) =− Az 2 − 0.3µ  1+ µ(5)

uncharged species : log ( γ ) = 0.1µ3 Only relevant for the Lawrence Livermore National Laboratory aqueous model in the llnl.dat database

PHREEQC-HPx Course, Prague, September 3-4, 2020

SCK CEN/39560003

35
ISC: PublicSession 1
where γ is the activity correction coefficient, z is the charge, µ is the ionic strength and
A is a constant at a given temperature.
o

If –gamma is present, the equation from WATEQ is used:=log ( γ )− Az 2 µ
1 + Ba0 µ

(6)

+ bµ

where a0 is the Debye-Hückel a coefficient (<:flt:DHA>) and b the Debye-Hückel b
coefficient (<:flt:DHB>), and B is a constant ag at given temperature (other symbols as
defined before).
SOLUTION_SPECIES<:txt:reaction>log_k <:flt:log_k>delta_h <:flt:enthalpy> [<:txt:units>]-analytic<:flt:A1><:flt:A2><:flt:A3><:flt:A4> <:flt:A5> <:flt:A6>-gamma <:flt:DHA> <:flt:DHB>
Code 5 – Syntax of SOLUTION_SPECIES

Figure 16 – Extracts of SOLUTION_SPECIES from the phreeqc.dat database (PHREEQC 3.6.2)

More information on the keyword SOLUTION_SPECIES is available on the online manual page.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

36
ISC: PublicAqueous speciation calculations

PHASES
Phases defines the name of mineral phases and gas components together with the reaction,
thermodynamic information (equilibrium constant and temperature dependency) and supplementary
information (e.g. molar volumes).
Code 6 – Syntax of PHASES shows some identifiers in SOLUTION_SPECIES; Figure 18 shows an extractfrom the phreeqc.dat database with some examples:•••Line 3: name: The name of the mineral (alphanumeric, no spaces)Line 4 : reaction: Dissolution reaction from phase to aqueous species. The chemical formula ofthe phase should be at the first position with a stoichiometric coefficient of 1.
Lines 5-7 have the same meaning as for the aqueous species (see page 35).

Gas components are defined in the same way. To discriminate them from minerals, the names contain(g) in most databases. HPGUI uses also (g) to identify a phase as a gas component. Gas components canbe used to fix the fugacity in the system. Explicit definition of a gas phase is with the keyword
GAS_PHASE.
PHASES<:txt:name><:txt:reaction>log_k <:flt:log_k>delta_h <:flt:enthalpy> [<:txt:units>]-analytic<:flt:A1><:flt:A2><:flt:A3><:flt:A4> <:flt:A5> <:flt:A6>
Code 6 – Syntax of PHASES

Figure 17 – Extracts of PHASES from the phreeqc.dat database (PHREEQC 3.6.2). Lines 903-930 show minerals;
Lines 1138-1158 show gases.

More information on the keyword PHASES is available on the online manual page.
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

37
ISC: PublicSession 1

SOLUTION
SOLUTION – Pure solution
SOLUTION 1 pure water
temp 25
pH 7
pe 4-units mmol/kgw-density 1.0-water 1.0
Code 7 – Definition of a pure solution

SOLUTION <:int:number> <:txt:description>
temp <:flt:temperature degree>
pH <:flt:pH>
pe <:flt:pH>
units <:txt:unit>
density <:flt:density>-water <:flt:mass>
Code 8 – Syntax of SOLUTION – Properties of the solution
SOLUTION – Units and Composition
SOLUTION 1 pure water
units mmol/kgw
Na 1
K 1
Code 9 – SOLUTION composition
UNITS
3 groups(1) per liter /L(2) per kg solution /kgs
(3) per kg water: /kgw
In mol or in gram with prefix m(milli) or u(micro)
In /kgs: ppb, ppm, ppt

Composition
Element namemust correspond to the items in the first column
in SOLUTION_MASTER_SPECIES
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

38
ISC: PublicAqueous speciation calculations

SOLUTION 1 pure water
units mmol/kgw
Na mg/kgw
Code 10 – SOLUTION definition with different concentration units

Concentration
Note that other units of a single group can be used
These concentrations units are only used for the definition of the initial solution composition.
All output is in mol/kgw.

More information on the keyword SOLUTION is available on the online manual page.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

39
ISC: PublicSession 1

Speciation model: Sea and ground water
Problem Definition
Calculate the saturation index of calcite in sea water and ground water at 25°C. The composition of these
two types of water is given below (Table 4, from Appelo and Postma (2005), Ex. 5.16).
Table 4 - Composition of sea and ground water (units mg/l, alkalinity defined as HCO3-)
pH

Na+

K+

Mg2+

Ca2+

Cl-

Alkalinity

SO42-

Sea water

8.18

10000

423

1180

410

17700

180

2610

Ground water

7.1

14

3.3

7.1

80

34

253

1.8

The PHREEQC calculation gives also answers to following questions:
1.
2.
3.
4.
5.
6.
7.

Are the solution charged balanced?
What is the ionic strength of the two water types?
What are the molality and the activity of the aqueous complexes HCO3- and Ca2+?
What is the activity correction coefficient of HCO3- and Ca2+?
Which aqueous species of the different components is the most abundant?
What is the mass action constant for the reaction CO2 + H2O = HCO3- + H+?
What is the saturation index of calcite and aragonite?

Hands-on
The solution of this problem requires only the definition of two SOLUTION data blocks separated by the
keyword END. A first data block defines the composition of sea water; the second data block defines thecomposition of ground water.

Create the simulation caseOpen the Project Manager
Select the Course project
Click Create New CasebuttonGive Name – “Example CE_a 01” 4
Give Description (optional) – “Aqueous speciation of sea and ground water”
Click OK button
Select the simulation case
Click Select button

4 C: Geochemical problem; E: Equilibrium calculation, a: only aqueous specation

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

40
ISC: PublicAqueous speciation calculations

Select database and components
Go to the Database tab
Click Select button
Choose the database phreeqc.dat in the most recent Phreeqc database group (now phreeqc 3.6.2)
Click Select button
The pop-up shows the aqueous components that are in the selected database. You can selectthose that are relevant for the simulation case, but this is optional for pure geochemical
calculations. However, syntax colouring of components, species and minerals in the input
workspace is based on the selected components.
Select relevant components: Na, K, Mg, Ca, Cl, S and C (based on the solution composition given
in Table 4.
Click OK button

Input Work Space
Figure 18 shows the input work space for a geochemical project created with HPGeochemistry. Bydefault, there are three tab pages:•••

Database: To select the database. When no database is selected, this tab page is the only tab
page that is accessible. After selection of the database, one can select components. There is alsoa read-only print out of the database.
Additions to database: an input tab that is used typically to define new components, species orrate equations.
Geochemistry: an input tab to define the geochemical input and the output.

Figure 18 – View of the Input Work Space

The text editor has a few advanced text editors options, such as:•

Actions as undo, redo, find, goto, etc.
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

41
ISC: PublicSession 1
The text editor is based on Fast Colored Text Box and inherit many functions from it.
See https://github.com/PavelTorgashov/FastColoredTextBox for more details.
o Some of these are accessible via the text editor toolbox.
Syntax highlighting, including coloring of numbers, text between quotes, inline variables andcomponents, species and names of phases (given that components are selected in the databasetab page) (Figure 19)
Autocompletion
Bracketing matching for () and []Multiple line (un)commenting (# in the text editor toolbox)
Different insert option (see below and Figure 20) – the gray text between angle brackets should
be replaced by actual input; otherwise errors will be generated when executing the code
Folding of parts of input file by pairs of:
o

•

•
•
•
•
•

o
o

#region and #endregion#{ and #}Figure 19 – Text editor in the Input Work Space

The toolbar of the input work space has following functionalities:•••

Options to open, save and close additional input files. Note that opening additional input files
in the input work space does not imply that these input files are included in the simulation.
These files have to be explicitly included in one of the two default tabs (additions to database,geochemistry) using the keyword INCLUDE$.
Insert options: Double clicking an item on the tree will insert text in the text editor. The format
depends on the selection in the insert options (Figure 20)
Autocompletion: To select the autocompletion option, for example only keywords, all valid
phreeqc statements, species, phases etc.

Figure 20 – Illustration of the different insert options. The grey text between angle brackets should be replaced byintegers, float, text or strings.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

42
ISC: PublicAqueous speciation calculations
The tab pages right to the text editor gives access to different trees:••••

Outline: The outline of the input file based on the PHREEQC keywords
HP Keywords: List of valid keywords, identifiers and options organized in different groups. Also
information from the database is available for the components selected in the database tab.
CPS: List of Components, Phases and Species using information from the selected database andthe components selected in the database tab.
Scripting: List of valid script functions (MyBasic). Also information from the database is available
for the components selected in the database tab.

After double clicking some item in the lists, an appropriate text depending on the insert options will beinserted at the current position in the text editor.

Sea Water
Go to Geochemistry tab
Define solution 1 with the keyword SOLUTION 1 Sea Water (i.e., solution number 1 and solution
description Sea Water)
Define the pH with the identifier pH
Define the units of the solution with the identifier –unit. The unit is mg/l.
Note that the default unit in PHREEQC is mmol
Define the list of components in the solution together with the concentration:
The list of components consists of the element name of an aqueous component as defined under
the data block SOLUTION_MASTER_SPECIES. Components to be added are Na, K, Mg, Ca, and
Cl. Also Alkalinity and S(6) need to be defined (see below).
Na+ is not a valid input in the SOLUTION data block (not an element name); Na is the valid
input (see line 19 in Figure 21).

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

43
ISC: PublicSession 1

Figure 21 - SOLUTION_MASTER_SPECIES data block in the phreeqc.dat thermodynamic database. Line 14
illustrates that the element name to be used in SOLUTION data block should be Na (first column, element name),and not Na+ (second column, master species). Line 30 shows the master species Alkalinity.

Alkalinity is a master species (see line 35 in Figure 21). The formula which is used to convert mg(as defined in the input) to mol is Ca0.5(CO3)0.5 (line 35 in Figure 21). In the problem definition,
alkalinity was defined as HCO3 (Table 4). Therefore, a new formula in the input indicates how mg
is converted to mol for the component Alkalinity using the identifier as:
Alkalinity 180 as HCO3
The molar weight of Ca0.5(CO3)0.5 is 50.04 g/mole. The molar weight of HCO3 is 61.092. When
HCO3 is not defined in the PHREEQC input, the moles of alkalinity of the solution are calculatedusing the molar weight Ca0.5(CO3)0.5. This gives a higher initial molality for alkalinity
SO42- is defined by the secondary master species S(6) (line 36 in Figure 21). S(6) is defined as mg
of SO4 Table 4) which is identical to the formula used for conversion defined in thethermodynamic database (in Figure 21). There is no need to add the as identifier.
Put the keyword END.
When a PHREEQC input file is executed, PHREEQC first reads the input file until the keyword END
is encountered. Before reading further, PHREEQC performs the read instructions. In this case, thisis the composition, properties and concentration of aqueous complexes in the sea water.

Ground water
Define solution 2 with data block SOLUTION, solution number 2 and solution description Ground
Water.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

44
ISC: PublicAqueous speciation calculations
Define the composition of ground water in the same way as is described for sea water

Run the simulation case
Click the ExecutebuttonIf there are errors or possible errors in the input file, both the screen output and the output file willcontain a list of error or warnings. Some typical errors and warnings with SOLUTION are listed below.

•

Warning message: Unknown input, no keyword has been specified
Check spelling of keyword (SOLUTION)•

Warning message: Could not find element in database, XXX. Concentration is set to zero.
Check spelling of XXX (should be an element name, Figure 21) or check database(SOLUTION_MASTER_SPECIES)•

Error message: No element or master species given for concentration input
Check element list in SOLUTION data block. Probable the first letter of the element name is not acapital

Explore the output
The results of the PHREEQC simulation are summarized in an output file called phreeqc.out. The content
of the file can be explored in the output workspace accessible via
for:
•
•

(Figure 22). This workspace allows

Explore if there are any errors or warnings.
Navigation through the file via:
o The outline treeo The small toolbox.

Figure 22 – Output work space with the phreeqc.out file. At the left hand-side, there is a navigation tree.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

45
ISC: PublicSession 1
Phreeqc.out gives information on the two geochemical calculations (sea water and ground water). For
each simulation, the output reports the solution composition, a description of the solution, thedistribution of the aqueous complexes (the species), and the saturation indices. The simulation output
gives answers to questions as posed in the problem definition.

Solution composition
Figure 23 shows the solution composition of sea water. It gives the molality of the components[mol/kgw] and the amount of moles in the water phase. The input concentration (mg/l) is converted tomolality.

Figure 23 – Phreeqc.out – Composition of sea water.

Description of the solution – Question 1-2
The description of solution summarizes the main geochemical characteristics of the solution. It gives thepH, pe, specific conductance, temperature, electrical balance, ionic strength, amongst others (Figure 24).
This output block gives answers for Q1 and Q2. For sea water, the solution is charged balance (percent
error on electrical balance is 0.63) and the ionic strength is 6.230 × 10-1.

Figure 24 - Phreeqc.out – Description of the sea water solution. The red rectangle indicates the percent error on
the electrical balance (Q1); the green rectangle indicates the ionic strength (Q2).

Distribution of species (Q3-Q6)
Figure 25 shows the distribution of aqueous complexes in sea water. For each species, molality andactivity are given. The log10 of these quantities and the log10 of the activity correction factor (gamma) are
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

46
ISC: PublicAqueous speciation calculations
also reported. This output block gives the answer for Q3. The molality of HCO3- and Ca2+ in sea water are
1.997 × 10-3 and 9.378 × 10-3 mol/kgw, respectively. The activity of HCO3- and Ca2+ in sea water are 1.358
× 10-3 and 2.368 × 10-3 mol/kgw, respectively. The activity correction factor for these two species are
0.697 (= 10-0.168) and 0.252 (= 10-0.598), respectively (Q4). The activity correction factor of a two-charged
species is larger than that of a single-charged species. This is because the activity correction depends
on the squared of the charge of the species.

Figure 25 – Phreeqc.out - Distribution of species for sea water. The red rectangles indicate the molality and activityof two aqueous species (Q3). The blue rectangles indicate the log10 of the activity correction factor (Q4). The green
rectangles indicate the most abundant species of some of the components (Q5). The yellow rectangles indicate
the log10 activity of some species needed to calculate a mass action constant (Q6).

The aqueous species for each component are ordered from the most abundant to the least abundant(Q5). The uncomplexed (free) aqueous species is the most abundant for Na, K, Mg, Ca, Cl. HCO3- is themost abundant aqueous species for C(4), followed by aqueous complexes between HCO3- and Mg or
Na.
The ratio of the activity coefficients of species in a reaction gives the mass action constant. For thereaction CO2 + H2O = HCO3- + H+ (Q6), the mass action constant is given as: HCO3-   H + K=
[CO2 ][ H 2O](7)

log ( K=
) log(  HCO3- ) + log(  H + ) − log([CO2 ]) − log([ H 2O])

(8)

Or after log10 transformation:

Using the log activities in the section on the distribution of species (Figure 25), the value of log K is:log ( K ) =−2.867 − 8

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003(9)47
ISC: PublicSession 1

Saturation indices
PHREEQC reports also the saturation index of all minerals in the thermodynamic database which can beformed with the components in the solution (Figure 26). Besides the saturation index, also the log10 of
the ion activity product (log IAP) and the log10 of the mass action constant at the defined temperature(here 25°C) are given (log KT). Sea water is over saturated with respect to both aragonite (SI = 0.69)
and calcite (SI = 0.84).

Figure 26 – Phreeqc.out: Saturation indices for sea water. The red rectangle indicates thesaturation index of aragonite and calcite (Q7).

Answer
Q1-Q7
The answers on Q1 to Q7 for sea water were discussed above. For ground water, values for the different
questions are found in the PHREEQC output tab:Q1 Q2 Q3 Q4 Q5 Q6 Q7 -Percent error on electrical balance: 1.24
7.329 × 10-3
Molalities of HCO3- and Ca2+ are 4.050 × 10-3 and 1.919 × 10-3 mol/kgw. Activities of
HCO3- and Ca2+ are 3.712 × 10-3 and 1.353 × 10-3
Activity correction factors of HCO3- and Ca2+ are 0.916 and 0.704, respectively. Note thatthe activity correction factors in ground water are smaller than those in sea water ; this isbecause the ground water is more diluted with a lower ionic strength.
The most abundant aqueous species in ground water are the same as those in sea water.
The mass action constant is the same.
The saturation index of calcite and aragonite are -0.05 and -0.19, respectively.

Saturation indices
The saturation index of calcite in sea water is 0.84 (and thus the sea water is over saturated with respect
to calcite). The saturation index of ground water is -0.05 ; ground water is in equilibrium with calcite.

Interpretation
Because calcite is oversaturated in sea water, there is no thermodynamic equilibrium between calcite
and sea water. Precipitation of calcite in sea water is kinetically inhibited. One possible hypothesis is that
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

48
ISC: PublicAqueous speciation calculations
Mg, which has a five time higher concentration than Mg in ground water, inhibits calcite precipitation(Appelo and Postma, 2005) (Appelo and Postma 2005).

Charge balancing and target saturation index
Charge balance a solution
Using the identifier charge in the solution data block allows to adapt the concentration of an aqueouscomponent to obtain a solution that is charged balance. For example, one can calculate that pH of asolution when the concentration of the other components are known or the amount of a component to
add to reach a given pH.
Examples of charge balancing
Code 11 shows two example of charge balancing a solution:
SOLUTION 1: Adapt the pH of this solution to obtain charge balance. pH = 10.98
SOLUTION 2: Adapt Cl concentration to obtain charge balance at a pH of 7. [Cl] = 2 × 10-3 mol/kgw
Note that charge can be used only once in a SOLUTION data block.
SOLUTION 1
pH 7 charge
Na 3
C 2
END
SOLUTION 2
pH 7
Na 3
C 2
Cl 1
END

Code 11 – Examples to charge balance a solution

SOLUTION <:int:number> <:txt:description>
pH <:flt:pH> charge
Na <:flt:concentration>
C <:flt:concentration>

Code 12 – Syntax of a solution when using pH to charge balance

Hands on: Charge balance of a solution
Question 1: You are in a laboratory to make a solution. You add 1 mmol of NaHCO3 and 1 mmol of
Na2CO3 to a kg of water. Prior to the experiment, you want to know the pH of the solution.
[solution: pH 10.106]
Question 2: You want to start your experiment with the solution from question 1 but first change to pH
to 7. How much HCl or KOH do you have to add to get that target pH?
[solution: 1.348 mmol of HCl or KOH]
Create the case CE_a 02 and implement the simulation

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

49
ISC: PublicSession 1

Adjust to a saturation index of a phase
When defining a solution, in some cases one wants to adjust the concentration in equilibrium with phaseor to a specific saturation index:•••

If you can assume that a solution is in equilibrium with some mineral (e.g. calcite)
To estimate the concentration of a solution in equilibrium with a mineral
To adapt the solution in equilibrium with atmospheric conditions

Examples of adjusting to a saturation index
Code 13 illustrates how a solution is defined in equilibrium with calcite and a pCO2 of 10-3.5 atm byadjusting respectively the Ca and C concentration.

SOLUTION 1
pH 7 charge
Ca 1 calcite
C 1 CO2(g) -3.5
Code 13 – Example of a SOLUTION definition in equilibrium with different phases
Initial value must be different from zero.

SOLUTION <:int:number> <:txt:description>
pH <:flt:pH> charge
Ca <:flt:concentration> <:txt:mineral> <:flt: targetSI>
C <:flt:concentration> <:txt:mineral> <:flt: targetSI>
Code 14 – Syntax of a solution in equilibrium with a solid

Hands on: Phase equilibrium with a solution
Question: What is the pH, the Ca and the C concentration of a solution in equilibrium with calcite and apartial pressure of CO2 of (a) 10-2 atm and (b) 10-3.5 atm?
[solution: (a) pH 7.297, Ca 1.646 mmol/kg, C 3.622 mmol/kg(b) pH 8.279, Ca 0.493 mmol/kg, C 0.980 mmol/kg]
Create the case CE_a 03 and implement the simulation

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

50
ISC: PublicSession 2 – Modelling Reaction Paths

Session 2
A key element in the study of equilibrium geochemical systems is to simulate how a system change
when its composition, pressure or temperature changes. Tracing the geochemical state of a system as aconsequence of irreversible reactions or processes by modelling is referred to as reaction pathmodelling. For systems involving only equilibrium reactions, the evolution of states is expressed in termsof a kind of reaction progress variable. When kinetic reactions are relevant for the geochemical system,
time is the progress variable. This type of reaction paths is called kinetic reaction path and will be
discussed in Session 5 – Kinetics (page 159). Spatial coordinates / information are lacking in reaction
path models; modelling of geochemical systems where both time and spatial information are relevant isindicated with the term reactive transport models. Reactive transport models are introduced in Session
3 – Reactive transport modelling (page 87) and used in the subsequent sessions as well.
We use here the term reaction path models in a broad sense where it can refer to more theoreticalstudies to understand a geochemical system (e.g., calculating an isotherm) or to more real, physicalactions involving one or another mean of irreversible mass transfer (e.g. titration) – of course, a specificmodelled reaction path may serve both goals. As for speciation models, they combine the principle ofmass balance with thermodynamic equilibrium. Some examples of equilibrium reaction paths are:MixingTwo or more geochemical systems are mixed together. This reaction path is also
fundamental to some other reaction paths, e.g. titration.
Example: Reaction path: Mixing two solutions on page 53
Titration
A system in which reactants are gradually added to or removed from the system.
From a practical application, this type of reaction path is relevant in chemical engineering when
the term batch reactor is used (Levenspiel 1999). Different variants are possible (Figure 27).
Example: Reaction path: Sliding concentration path on page 55
Example: Reaction path – Exchange isotherm on page 140

Figure 27 – Variants of titration reaction paths adding reactants (batch reactors in chemical engineering,
Levenspiel (1999).

Evaporation – reaction path in which only water (H2O) is removed.
Example: Reaction path – Changing exchange composition during evaporation
on page 145
Buffering

The activity of a gas or an aqueous species is fixed. Several variants exist:
Fixed fugacity path – one or more gases are buffered
Sliding fugacity path – fugacity of one or more gases change
Fixed or sliding activity paths – activity of one or more aqueous species are bufferedor changed. The most typical form is a speciation plot as a function of pH.
Example: Reaction path: Sliding activity path on page 66

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

52
ISC: PublicModelling reaction paths
Flush
Tracing the evolution of a solid phase when unreacted fluid is added to the system.
After each step, the solid phase is separated from the aqueous phase and reacted with a newvolume of the aqueous phase.
Flow-throughTracing the evolution of an aqueous phase. After each time step, the solid
phase is separated from the aqueous phase, and further reactions in the aqueous phase are
simulated.
Many theoretical and practical applications are a combination of two or more of these basic reaction
paths (or combined with kinetic reaction paths).

Learning outcomes
Use PHREEQC for simulating irreversible reactions
Understand the concept of inline variables in HPGeochemistry
Use HPGeochemistry specific keywords for repetitive geochemical calculations and plotting
Implement different basic reaction path simulation (mix, sliding concentration, sliding activity)

Reaction path: Mixing two solutions
Learning outcome
Use the PHREEQC keyword MIX

Problem definition
Two solutions in equilibrium with calcite, but in a different atmospheric condition, are mixed. How does
the saturation index of calcite change?

Conceptual model
The problem involves following components:
•
•
•
•
•

Open system
Equilibrium reactions
One phase (aqueous)
Well-mixed
Physical action of mixing two solutions

Figure 30 shows a schematic representation of mixing of two solutions.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

53
ISC: PublicSession 2

Figure 28 – Schematic representation of mixing two aqueous phases.

Implementing the model
Figure 29 shows the PHREEQC keywords for implementing this problem. The keyword MIX is introducedin this section.

Figure 29 - Implementation of the model for mixing two solutions.

MIX
MIX indicates which solutions will be mixed, and the fraction of each solution that is taken in the new
solution. In Code 15, 50 % of solution 1 is mixed with 25% of solution 2. The new solution is calculated
by adding these two solutions. If both solutions consists of 1 kg of water, the new solution has a total
amount of 0.75 kg water. If Na in solutions 1 and 2 is 0.5 and 0.1 moles, respectively, the total amountof Na in the new solution is 0.275 moles. Code 16 shows the syntax of the MIX data block.
MIX 1
1 0.5
2 0.25
Code 15 – Mix solution 1 and solution 2 with a fraction of, respectively, 0.5 and 0.25 of the original solutions.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

54
ISC: PublicModelling reaction paths
MIX <:int:number> <:txt:description><:int:solution number> <:flt:fraction><:int:solution number> <:flt:fraction>
Code 16 – Syntax of MIX
More information on the keyword MIX is available on the online manual page.
In the graphical user interface of HPGeochemistry, MIX is found at:
HPKeywords -> Physical Action -> MIX

Hands-on: Mixing two solutions
Solution 1 is equilibrated in an atmosphere with a partial pressure of CO2 of 10-2 atm, solution 2 in an
atmosphere with a partial pressure of CO2 of 10-3.5 atm.
Mix these solution in following fractions:•••

25% of solution 1 with 75% of solution 2
50% of solution 1 with 50% of solution 2
75% of solution 1 with 25% of solution 2
Create the case CE_a 04 and implement the simulation[Solution: SI of calcite increase from -0.3, -0.22 to -0.11]

The new solution is under saturated with respect to calcite, although we started from two solutions
saturated with calcite. Although mixing is a linear process (just add masses of water and moles of
elements), the activity and speciation calculation is nonlinear. It may lead to dissolution of calcite.

Reaction path: Sliding concentration path
When a reactant is added to a system, the chemical entities (the aqueous phase and the assemblages of
exchange sites, surfaces, and phases) will change. The keyword REACTION allows for adding a numberof components to the system in a step-wise manner.

Learning outcomes
Use REACTION to describe titration experiments
Use EQUILIBRIUM_PHASES to fix pH, impose partial pressure of gases and impose equilibration
with pure minerals
Use SELECTED_OUTPUT to control the output of the geochemical calculations
Implement a sliding concentration reaction path
Construct a sliding activity plot using post-processing in HPGeochemistry.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

55
ISC: PublicSession 2

Speciation of Cd as a function of Cl concentration
Problem definition
Cd2+ forms aqueous complexes with OH- and Cl-. The chemical form of an element determines its
mobility, transfer, fate and toxicity. In this particular example, we will invest the speciation of Cd2+ at a
specific pH as a function of Cl concentration.
Consider a system of 1 kg of water and a concentration of 1 mM of Cd. We will add in total 1 mol of HCl
to the system, and buffer the solution up to a pH of 3 by adding NaOH.

Conceptual model
The problem involves following comopnents:
•
•
•
•
•
•

Open system
Equilibrium reactions
One phase (aqueous)
Well-mixed
Adding reactant – only elements
Buffer pH

We model a Sliding concentration path with fixed activity.
Figure 30 shows the mathematical representation of the problem. There is one phase with equilibrium
reactions and activity corrections defined in a database (phreeqc.dat). We define a buffer (solid buffer)
to fix the pH at the requested value. We incrementally increase the Cl concentration by adding HCl to
the system.

Figure 30 – Schematic representation of a sliding concentration path with fixed activity.

Implementing the model
Figure 31 shows the required input and output keywords to implement a model that calculates
speciation as a function of concentration at a fixed pH. Following new keywords are required:
EQUILIBRIUM_PHASES and REACTION. In addition, the keyword SELECTED_OUTPUT is introduced tocontrol the output.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

56
ISC: PublicModelling reaction paths

Figure 31 – Implementation of a sliding concentration path with fixed pH

EQUILIBRIUM_PHASESThe keyword EQUILIBRIUM_PHASES defines an assemblages of pure equilibrium phases.
EQUILIBRIUM_PHASES 1
Calcite 0 10
Calcite 0 Gypsum 1.0
Fix_H -5 HCl 10
Code 17 – Definition of an EQUILIBRIUM_PHASES.
Code 17 shows an example of the EQUILIBRIUM_PHASES keyword:•••Line 2 gives the standard use together with the defaults (see Line 2 in Code 18 for syntax): Nameof mineral, target saturation index (zero is forcing equilibrium between phase and solution) andamount (in moles of mineral in your system).
Line 3 attains the saturation index of the mineral (calcite) by adding or removing the elementsof a given formula (here defined by the chemical formula of gypsum)Line 4: This option is also useful to force a given pH at your system (line 4) - it is required todefine Fix_H in the phases data block (see Code 19).

EQUILIBRIUM_PHASES <:int:number><:txt:description><:txt:phase name> <:flt:SI> <:flt:amount (moles)>Code 18 – Syntax of EQUILIBRIUM_PHASES.
PHASESFix_H; H+ = H+; log_k 0
Code 19 – Definition of a phase to fix pH of a system.

More information on the keyword EQUILIBRIUM_PHASES is available on the online manual page.
In the graphical user interface of HPGeochemistry, EQUILIBRIUM_PHASES is found at
HP Keywords -> Geochemical entity ->EQUILIBRIUM_PHASES
To include the phase name and saturation index, use the subnode (only when components are defined):
Si and AmountA template to define the Fix_H phase (Code 19) is found at:
HP Keywords ->Database -> PHASES -> Templates -> Fix the pH
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

57
ISC: PublicSession 2

REACTION
The keyword REACTION defines an irreversible reaction that transfer element to or from the aqueous
solution.
REACTION 1
NaCl
2
Calcite -0.001
1 moles in 100 steps
Code 20 – Definition of an REACTANT with reactants added incrementally in equal amounts.
Code 20 shows an example of the REACTANT keyword, the syntax is given in Code 21:•••

Line 2 and 3 define the reactant to be added (NaCl) or removed (Calcite) from the system.
Reactants are defined as a chemical formula (NaCl) or as a phase name (Calcite).
The numerical value after the reactant indicates the relative stoichiometry compared to the otherreactants. A negative number indicates a removal of the elements.
Line 4 defines the amount of the reaction and the number of steps.

REACTION <:int:number> <:txt:description><:txt:reactant> <:flt:relative stoichiometry><:float:amount> moles in <:int:steps> steps
Code 21 – Syntax of REACTION.
Instead of defining a reactant amount with constant increment, a list of amounts can be defined as well.
Line 4 in Code 22 defines four steps. The total amount of moles added to the system depends on thestate of the INCREMENTAL_REACTIONS keyword. If false (INCREMENTAL_REACTIONS false) whichis the state by default, PHREEQC will add the amount for each step to the starting solution, meaning
that, in this particular example, we obtain results for systems in which 0.25, 0.45, 0.95 and 1 moles of
NaCl (and corresponding amounts of calcite being removed) are added. If true
(INCREMENTAL_REACTIONS true), the next step is added to the solution obtained in the previous step
such that the results are for systems in which 0.25, 0.7, 1.65 and 2.65 moles of NaCl are added (andcorresponding amounts of calcite being removed).
REACTION 1
NaCl
2
Calcite -0.001
0.25 0.45 0.95 1 moles
Code 22 – Definition of an REACTANT with an explicit list of added reactants.
Note that INCREMENTAL_REACTIONS have also a computational effect when reactants added
incrementally as in Code 20. However, the total amount added is always 1 for the example in Code 20.
When INCREMENTAL_REACTIONS is false, each step calculates the system as the amount at a certain
step is added in the starting solution. When INCREMENTAL_REACTIONS is true, the incrementallyamount is added to the resulting solution from the previous step. 5
Note also that one can use millimoles or micromoles instead of moles.
More information on the keyword REACTION is available on the online manual page.
More information on the keyword INCREMENTAL_REACTIONS is available on the online manual page.

5

INCREMENTAL_REACTIONS has a similar effect on the time steps in a KINETIC calculations.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

58
ISC: PublicModelling reaction paths
In the graphical user interface of HPGeochemistry, REACTION is found at
HP Keywords -> Physical Actions -> REACTION
Linear steps (e.g. Code 20) or an explicit list (e.g. Code 22) are found at the subitems, respectively
Linear steps
andLists of steps
Within each of these subitems, the graphical user interface gives support for including reactants basedon phase names (only when components are defined):
Reactant - phases
And the syntax for the amounts of the reactants (for different units):
units
INCREMENTAL_REACTIONS is found at
HP Keywords -> Physical Actions -> INCREMENTAL_REACTION
SELECTED_OUTPUT
The keyword SELECTED_OUTPUT specifies the geochemical variables that will be written to a output
file that consists of columns (for variables) and rows (for different calculations). Such data files allows fora more easy data manipulation in spreadsheets and/or plotting programs. SELECTED_OUTPUT has
many identifiers to specify the type of geochemical data to be exported, ranging from information on
the calculation (e.g., -state, -simulation), general geochemical information of the system (e.g. –pH,
-pe, -charge_balance), concentrations of components (-total), molalities or activities of species(-molalities), -activities)), and information on the different chemical entities (e.g., equilibrium_phases), -gases) – see the online manual page.
Code 23 shows an example of the SELECTED_OUTPUT keyword:••Line 1: a number (integer) is associated with the SELECTED_OUTPUT. As such, different output
files can be defined each identified with that number.
o Note that when the number is 1, a number of identifiers will be true by default. For
numbers different from 1, the default will print no data to the output file (all identifiers
are false by default). For example, Line 3 is not needed in the SELECTED_OUTPUT 1 as
the pH is written by default to the output file. However, for a SELECTED_OUTPUT n (ndifferent from 1), the pH will only be written to the output file when Line 3 is present.
o The number of the SELECTED_OUTPUT links the specific output file (name defined inLine 2) to two other keywords with the same number: USER_PUNCH which will write
additional data columns to the output file based on a user-defined script, and
SELECTED_GRAPH will produce plots based on the data in the outputfile (the latter isonly available in HPGeochemistry and HPx, not in other versions of PHREEQC) (seelater).
Line 2: -file) defines the name of the output file. Implicitly, it defines also the folder, as
HPGeoChemistry writes the file relative to the case folder. Thus, in Line 2, the file is written to
the case folder. When –file ..\results.out is used, the file is written to the parent folderof the case. Note that, if the file already exists, it will be overwritten.
SELECTED_OUTPUT 1-file results.out-pH true–totals Na K
-molalities CO2 HCO3- CO3-2–equilibrium_phases Calcite
Code 23 – Definition of a SELECTED_OUTPUT keyword.
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

59
ISC: PublicSession 2

Figure 32 – Output file obtained with Code 23 (for a given geochemical model) opened in HPGeochemistry.

There is a large flexibility for writing these output files:
•
•

•

•

A single input file and set of calculations can have different output files identified by the
numbers.
Beside the main geochemical variables that can be defined with the identifiers in
SELECTED_OUTPUT, other output values can be written as well by using the USER_PUNCH
keyword. The USER_PUNCH keyword will be used in other exercises.
During the calculations, writing to the output file can be switch off and on:o Within the keyword PRINT, the identifier –selected_output) will switch off (if false)or on (if true) writing to all defined selected_output files.
o Within the keyword SELECTED_OUTPUT, the identifier –active) will switch off (if false)or on (if true) writing to this specific_output file.
o Within the keyword SELECTED_OUTPUT, the identifier –punch) will switch off (if false)or on (if true) writing of the output defined with PUNCH having the same number as thisspecific output file.
If a SELECTED_OUTPUT with the same number is defined for a second time in a single input
file:o If only –active and/or –punch are defined: the status of writing is changed.
o If any other identifier is defined: the data file is closed, and the existing
SELECTED_OUTPUT definition is replaced by the new one. Note that if you have thesame file name, the existing one is overwritten.

These options allows for smart writing to output files that can be directly used for e.g. plotting.
More information on the keyword SELECTED_OUTPUT is available on the online manual page.
In the graphical user interface of HPGeochemistry, SELECTED_OUTPUT is found at
HP Keywords -> Output -> SELECTED_OUTPUT
Some of the identifiers with their options are also available.

Hands-onCreate the simulation caseOpen the Project Manager
Select the Course project
Click Create New CasebuttonGive Name – “Example RE_a 01” 6
Give Description (optional) – “Speciation of Cd as a function of Cl at fixed pH”
Click OK button

6 C: geochemical problem; E: Equilibrium calculation, a: only aqueous speciation

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

60
ISC: PublicModelling reaction paths
Select the simulation case
Click Select button

Select database and components
Go to the Database Panel
Click Components button
Select relevant components: Cd, Cl, Na
Click OK button

Create input
The input file consists of the (i) solution definition (SOLUTION), (ii) the reaction definition(EQUILIBRIUM_PHASES, REACTION), and (iii) the output definition (SELECTED_OUTPUT). In addition,a new phase needs to be defined to fix the pH at a value of 3.
Addition to Database

Go to Additions to Thermodynamic Database Panel
Type the content of Code 19
or
Select Help tree HPKeywords ->Database ->Phases ->Templates ->Fix the pH

Solution Definition
Go to Geochemistry Panel
Define SOLUTION 1 with the Cd concentration and pH – charge balance by Cl

Reaction Definition
Define EQUILIBRIUM_PHASES to fix the pH at a value of 3 by adapting NaOH (buffer)
Define REACTION
Add HCl with a stoichiometric coefficient
Define the addition of 1 mol in total in a number of steps (e.g. 100)
Note that we add here only the components H and Cl to the system, without changing the volumeof the system. This is equivalent to simulating a batch reactor with a constant volume as was
illustrated in the most left sketch in Figure 27.
Controlling output and plotting
SELECTED_OUTPUT is used to select relevant output results. To keep the input file clear, we also makean additional input file.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

61
ISC: PublicSession 2
Add SELECTED_OUTPUT 1000
Choose file name (e.g. CdSpeciation.out)
Set pH true
Add –totals for Cl
Add -molalities for the Cd species (select 9 at max)
Run (

)

OutputPhreeqc.out contains one initial solution calculation and 100 reaction steps in the batch simulation. For
each reaction step, the output file describes what geochemical entities are used (solution 1, pure phaseassemblage 1 and reaction 1), how much reactant has been added, and the geochemical speciation.
It is more useful to look at the results in the file defined in the SELECTED_OUTPUT data block(CdSpeciation.out). This is possible via the HPGeochemistry graphical user interface by selecting theoutput file workspace.
Click Output file workspace (Click Open data file () button) button

Select the output file (CdSpeciation.out) – by default, the project folder opens.

Datafile workspace
The data file opens in a grid view assuming the data file structure used by Gnuplot:•••

First line contains headings
Next rows contain the data
Two datasets are separated by a double empty line

File Structure Ribbon
The user can define different data structures using the File Structure ribbon (via). In the course, wewill only work with the standard structure. The File Structure ribbon can also be used to select text view
instead of grid view.
Graph Ribbon
Graphs can be created from a data file using the Graphs ribbon (viaribbon:••••). Figure 33 shows the Graph

Create graph:o To select the type of graph (XY or YX graph)o To define the location of the graph; either as a floating window or as an item in theGraph work space (see later)o Button to create the graph
Select X: Selection of the X-variable (in a XY graph)
Select Y: Selection of Y-variables (in a XY graph)
Select DataSet: only active in a multi dataset to select the dataset

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

62
ISC: PublicModelling reaction paths

Figure 33 – Output file workspace: Graph ribbon

We create an XY plot in which the speciation of Cd (only Cd+2, and the three Cd-Cl complexes) is plotted
as a function of Cl concentration.
Select Cl in Select X
Select m_Cd+2, mCd_Cl+, mCdCl2, mCdCl3- in Select Y
Click Create Fig button
Note that X and Y selections can be done by selecting the columns in the grid view as well.
The Graph Settings dialogue window opens with following tabs:•••

General: Graphs are ordered at two levels in the Graphs work space (Group and Type). In
addition, a name can be given.
Axis: Some settings for the X and Y axis. Note that, by default, the axis labels are constructedfrom the column headings in the data files. These labels can be redefined here. Enhanced text
options should follow the gnuplot format and can be changed by gnuplot commands (see later).
Curves: Some settings for the curves, such as the labels for the legend. Again, enhanced text
options of gnuplot can be used.

Code 24 shows example of sub- or superscript enhanced text options in gnuplot. Subscripts are denotedby _ ; superscripts by ^. Curly brackets {} are used for sub- or superscripts of multiple characers. Using@ will align sub- and superscripts.
Gnuplot enhanced text
Cd_{2+}
CdCl^CdCl_3^CdCl@_3^-

Results
Cd2+
CdClCdCl3CdCl-3

Code 24 – Examples of enhanced text options in gnuplot – Superscripts and subscripts.

Go to General
Type Cd-Cl speciation in textbox Name
Go to Axis
Go to X
Type “Cl [mol/L]” in textbox Title of Axis
Go to YType “Molalities [mol/L]” in textbox Title of Axis
Type 1E-6 in Minimum Value
Type 1E-3 in Maximum Value
Select Log scale
Go to Curves
Remove in the column “Labels” m_ from each entry
Click OK
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

63
ISC: PublicSession 2
The figure is created and listed in the items of the Graph Workspace.
Click Graph work space () button

Graph workspace
Figure 34 shows the Graph work space (after selection of the figure Speciation Cd-Cl). The tree at theleft side lists all figures defined within the simulation case (e.g., figures made from output file or figures
defined in the input file). The toolbar above allows for:•View multiple selected graphs•

Navigate through the list

Delete the figure (not for all figures)•

The right side shows one or more graphs (depending on the selection). Each graph have a toolbar to:••

Edit the graph with gnuplot commands and some limited changes via the graphical userinterface
Copy to the clipboard

Export the graph via the activated Graph Export ribbon.

•

Figure 34 – Graph work space – General View

Click the Edit Graph () button or Go to the Edit Graph tab

For a graph made from an output file, there are two tabs:••

Commands: a text editor to enter gnuplot commands (see gnuplot help file)
Curves: To change properties of the individual curves

Add “set key inside bottom” in the Commands tab (position of the legend)
Change the labels of the curves in the Curves tab (use enhanced text options for the species)

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

64
ISC: PublicModelling reaction paths

Figure 35 - Speciation of Cd as a function of Cl with CdTotal = 1 mM and a pH = 3.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

65
ISC: PublicSession 2

Reaction path: Sliding activity path
Many geochemical entities and component exhibit large changes in their speciation or solubility of thefull pH range; up to many order of magnitude. Information on how speciation and solubility of different
elements are a function of pH is crucial information for transfer and fate of the element in theenvironment.
Solubility of Fe(OH)3 shows a so-called amphoteric behaviour with high solubility at low and high pH
and low solubility at intermediate pH values due to the speciation of ferric iron (Fe(3)). In addition, thebackground concentration may play a role as anions can form aqueous complexes with Fe3+ as well.
In this exercise, we will calculate the solubility and speciation of Fe in a system with Fe(OH)3 and abackground concentration of NaF. We will introduce a number of keywords (also a few specific for
HPGeochemistry) that enables to have a structured input file for calculating such so-called slidingactivity plots; these plots forms the basis for geochemical interpretation.
The approach is to a given extent similar to the one implemented in PhreePlot (Kinniburgh and Cooper,2009) which is a PHREEQC wrapper enabling a more straightforward way of making plots.
HPGeochemistry works with an element called inline variable that is similar to the concept of tags in
PhreePlot.
Using INCLUDE$, GLOBAL_VARIABLE and different forms of inline variables allows for a structured
input file for making this kind of fundamental geochemical plots.

Learning outcomes
Understand the concept of inline variables and apply them under the form of a
GLOBAL_VARIABLE, COUNTER and inline expression
Use INCLUDE$ and GLOBAL_VARIABLE to make a structured input file
Implement a sliding activity calculations
Produce one-dimensional sliding activity plots (SELECTED_GRAPH)

Solubility of Fe(OH)3 as a function of pH
Problem Definition
Calculate the speciation of Fe+3 and the solubility of Fe(OH)3 as a function of pH in a solution containing
0.1 mol Fe and a NaF background concentration of 10-1 M.

Conceptual model
The problem involves following components:••••

Equilibrium reactions
Two phases (aqueous, solids)
Well-mixed
Buffer pH at different values

We model a Sliding activity path with varying pH (activity of H+)..

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

66
ISC: PublicModelling reaction paths
Figure 36 shows a number of independent systems in which equilibrium is calculated within the aqueous
phase and between the aqueous and solid phase. A property is fixed at a given value by an external(solid) buffer (in the exercise, the pH will be fixed at different values in the different systems).

Figure 36 – Schematic representation of a sliding activity path to calculate the solubility of a mineral.

Implementing the model
Figure 37 shows the input and output keywords that can be used to implement a model that calculates
solubility as a function of pH. Following new keywords are discussed: INCLUDE$, GLOBAL_VARIABLE,
COUNTER and COUNTER_INCR. The latter 3 are keywords only valid in HPGeochemistry and are used tomake repetitive statements more easily; in this particular case, an equilibrium calculation of Fe(OH)3 at
each pH. Note that the same problem can be implemented with PHREEQC keywords only as well.
In addition, the keyword SELECTED_PUNCH is introduced to control the output and SELECTED_GRAPH
to produce an plot. Also SELECTED_GRAPH is a keyword specific for HPGeochemistry; it replace the
keyword USER_GRAPH from the classical phreeqc; USER_GRAPH does not work in HPGeochemistry.

Figure 37 – Implementation of a sliding activity path

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

67
ISC: PublicSession 2
Inline variables
An inline variable is a part of text in the input file delineated by the symbol @ followed by a text string.
During execution time, the inline variable is replaced by a corresponding numerical or alphanumericvalue that is defined in the input file or calculated during execution. Different types of inline variables
exist:-

Defined as a GLOBAL_VARIABLE with a fixed value
Defined in a DATA BLOCK
Defined as a COUNTER with a value that is changed after its use
Defined as an inline expression

Note that these keywords are not available in standard PHREEQC. Examples of some of these types are
given during the course.
GLOBAL_VARIABLE
A GLOBAL_VARIABELE defines the name (case-sensitive) and value of a variable that can be used as an
inline variable or within a script. The global variable can be either a numerical value or an alphanumericvalue. In case of an alphanumeric value, the name should end with $. Numerical global variables can be
defined either as a number or as a set of simple arithmetic expressions.
Code 25 shows the syntax of global variable:•••

Case sensitive name (with $ to indicate an alphanumerical value, Line 2)
The numerical value (or string “” for a alphanumerical value, Line 2)(optional) a description of the global variable
GLOBAL_VARIABLE<:txt:name><:flt:value>[<:str:description>
GLOBAL_VARIABLE<:txt:name>$<:str:text>[<:str:description>]
Code 25 – Syntax of GLOBAL_VARIABLE

To illustrate the use these inline variable, we implement some examples in a temporary case.
Open the Project Manager
Select the Course project
Click Create New CasebuttonGive Name – “tmp”
Because the input file is read and executed in a linear way, inline variables are only available when they
are previously defined. A rather good place to define them is in the top of the Additions todatabase Panel. We define a solution in which the pH, Cl and concentrations are defined by inline
variables.
Add Line 1-3 from Code 26 in the Additions to the database Panel
Add Line 5-8 from Code 26 in the Geochemistry Panel
Run
Check the output file – you can see that the inline variables @pH, @IS and @CdConc are replacedby their current values.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

68
ISC: PublicModelling reaction paths
GLOBAL_VARIABLE pH 7
GLOBAL_VARIABLE IS 1E-3
GLOBAL_VARIABLE CdConc 1E-6
SOLUTION 1-pH @pH charge
Cl @IS
Cd @CdConc

Code 26 - Example of the of inline variables defined with GLOBAL_VARIABLES

Inline variables are case-sensitive.
For example, using @Cdconc in Code 26 (Line 8) will give following error:
ERROR: Inline value for Cdconc is not defined.
@Cdconc
To use a set of simple numerical expressions 7, name_variable should be proceeded by the equal symbol(=). Single (Code 27, Line 2) or multiple expressions (Code 27, Line 3) can be used with & as theexpression separator. For the latter, the global variable is set to the last expression.

GLOBAL_VARIABLE radius 0.5 "radius of an object,in m "GLOBAL_VARIABLE area=3.1415*@radius^2 "area of theobject, in m2"
GLOBAL_VARIABLE area=pi=3.1415 & pi*@radius^2"area of the object, in m2"
Code 27 – Examples of defining a GLOBAL_VARIABLE with simple numerical expressions.
Define @CdConc ten times smaller than @IS
Check the output file.
These GLOBAL_VARIABLE keywords can be used anywhere in the input file, but, as mentioned above,
they can only be used after they have been defined. When you redefine a GLOBAL_VARIABLE, the oldvalue is replaced.
In the graphical user interface of HPGeochemistry, GLOBAL_VARIABLE is found at
HP Keywords -> Scripting -> GLOBAL_VARIABLE
HPGeoChemistry provides a graphical user form for defining inline variables via the Global Definition
panel. These are introduced in the input file in the beginning of the input using the above mentioned
GLOBAL_VARIABLE keyword. The form for defining the inline variables has some advantages:••••

The input for the values of the inline variables is via a nice tab and table format.
There is some limited assistance for writing the numerical expressions.
The inline variables defined via these forms are available in the help trees and the
autocompletion option in the input space.
Runs with different parameter values can be defined.

In the graphical user interface of HPGeochemistry, global variables defined via the Global Defintion Panelare available via

7 MuParserX is used as the mathematical parser (see https://beltoforion.de/en/muparserx/). For more complex

recalculations, CALCULATE_VALUES with –state definition can be used (see help file)
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

69
ISC: PublicSession 2
CPS -> Global Defintions -> Global Variables with @CPS -> Global Defintions -> Global Variables
The first option will include the name of the global variable with a @ symbol (inline variable); the second
will just include the name (e.g. to be used in the script function get_var(“”), see help file).
Select the Global Definition panel in the input workspace
Select the Global Variable item in the list at the right hand side
The form has two functionalities:“Edit” (Figure 38) to construct and define a input structure for defining the values of the inline variables.
The input structure consists of one or more tabs having each one or more input group. An input groups
consists of four columns and the user can define one or more records. Tabs, groups and records can beinserted, deleted or moved with the three icons strips.
The input structure is visible as a tree (Table). At the right hand side, one or two tabs are available when
you selected a first level node or a second or third level node, respectively. The properties tab allows
for including information as a tool tip and/or, a short description for the tab (first level node), group(second level node) or column (third level node).
During this course we will only have examples with the standard table (one tab and one group) and add
different records if needed. No extra information will be added.

Figure 38 - Form to define the inline variables

Select the second level node Variables
Go to the table and fill in the Name and Value of the three variables (pH, IS, CdConc). Optionally,
you can add the unit and description :Go to the Values tab to confirm the new table structure.
Go to one of the input tabs (e.g. GeoChemistry)
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

70
ISC: PublicModelling reaction paths
In the CPS-tree, under the node Global Variables, the defined global variables are in the list.
When you type @ in the text editor, the autocompletion will list these variables as well.
Go back to the Global Definition tab, Global Variables Item
Go to the Values tab“Values” (Figure 39) to define the actual values that will be used during the simulation.
Actual numbers can be defined here, however it is not possible to add new records in this form (this is
only possible in the Edit tab).
The tool strip in the left area allows to define different set of variables that will use the same input file.

Figure 39 - Form to define the values of the inline variables

Add two case where IS =1E-3 and 1E-2
Adapt the value of IS in the second case to 1E-2

HPGeochemistry will use the inline variables for the selected case. Select the first case, run thesimulation and check the output. Repeat for the second case.

COUNTER and COUNTERINCR
COUNTER is a special global variable: after it is used as an inline variable (or in an inline expression), thevalue is changed. The default change is adding 1, but this can be redefined using COUNTERINCR.
Code 28 shows the syntax of these two keywords. It is similar as GLOBAL_VARIABLE; also simple
numerical expressions can be used.

COUNTER <:txt:name> <:flt:starting value>
COUNTERINCR <:txt:name> <:flt: value>
Code 28 – Syntax of COUNTER and COUNTERINCR

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

71
ISC: PublicSession 2
When the counter is used as an inline variable (@name), its current value will be inserted and thenchanged by adding value; other options are possible by adding a prefix. For example in Code 29,
In line 4, Conc will not change the value due to the prefix =In line 5, Conc will change by dividing the value by 10 due to the prefix /In line 6, Conc is reset to its starting value due to the prefix &•••COUNTER Conc 0.1
COUNTERINCR Conc 10
SOLUTION 1
Na @=Conc
F @/Conc
COUNTER &Conc
Code 29 – Example of use of COUNTER and COUNTERINCR

Following prefixes are allowed:
nothing

counter is used, and then incremental value is added&

counter is used and reset to initial value+

counter is used, and then incremental value is added-

counter is used, and then incremental value is subtracted*

counter is used, and then multiplied by incremental value/

counter is used, and then divided by incremental value=

counter is used, value is not changed

In the graphical user interface of HPGeochemistry, COUNTER and COUNTERINCR is found at
HP Keywords -> Scripting -> COUNTERHP Keywords -> Scripting -> COUNTERINCR
INCLUDE$INCLUDE$ indicates the insertion point of an external file in the input file. The inserted file is interpreted
as part of the input file and insertion points may be within a keyword block, a simulation block or
between simulation blocks (adding additional blocks). The insertion is done dynamically at running time,
so one can write these input files and insert them at a later point during the execution. This functionalityof PHREEQC gives the possibility of advanced scripting aspects including generating significant input atrunning times.
The syntax is simple: the keyword followed by the name or the full path to the file (Code 30).
HPGeochemistry looks for the file in the simulation case folder. If your file is located there, justthe file name and extension is sufficient. If the file is located somewhere else, then a full path isrequired, except if alternative search folders are defined with INCLUDEPATH (see later in thecourse).

INCLUDE$ <:txt:Filename>

Code 30 – Syntax of INCLUDE$

More information on the keyword INCLUDE$ is available on the online manual page.
In the graphical user interface of HPGeochemistry, INCLUDE$ is found at
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

72
ISC: PublicModelling reaction paths
HP Keywords -> Scripting -> INCLUDE$
The sub nodes open a file exploring starting at the folder of the simulation case (Selecting File –From Case) or from the last accessed folder(Select File – Only Name and Select File –Full path, where the former will transfer only the filename in the current position of the input panel
and the latter the full path).

REPEATREPEAT will create a file with repetitive statements. The major goal of that file is to include that file at agiven point in the input (using INCLUDE$).
Code 31 shows the syntax for the repeat include statement. :file_in is the file name of a file that contains a number of statements that need to be executeda number of times.
file_out is the filename of the file that will be created
number is the number of times the file_in should be executed.

REPEAT include<:int:number><:str:file_in><:str:file_out>

Code 31 – Syntax of REPEAT

The file_out will contain number number of lines, each line having following statement:

INCLUDE$ file_in
The aim is that file_out is then included at a specific position in the input to have a repetitive
execution of the statements in file_in.
In the graphical user interface of HPGeochemistry, REPEAT is found at
HP Keywords -> Scripting -> REPEAT
USER_PUNCH
The keyword USER_PUNCH allows for converting the standard Phreeqc output variables into userdefined quantities and to write them to the file specified by SELECTED_OUTPUT. It also writes variables
to the output file that are not available in the SELECTED_OUTPUT. The conversion of output variables
is done with a script included in the USER_PUNCH block 8.
HPGeochemistry allows for two scripting language:•

The standard scripting language in PHREEQC, a kind of old-school, non-structured BASIC withnumbered Basic statements

The same scripting possibility is also used in a number of other keywords including USER_PRINT,CALCULATE_VALUES and RATES.
8PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

73
ISC: PublicSession 2•

An alternative Basic-type of language, MyBasic, developed by R. Wang (seehttps://github.com/paladin-t/my_basic), which is a more modern variant of Basic allowing forstructured programming with support of object-oriented program possibilities.

Both language have the same functions for geochemistry (although MyBasic requires that many
functions ave parentheses, e.g. mu() instead of mu), but MyBasic offers more flexibility for structuring
programs and advanced scripting (including templates). The Mybasic language is invoked by theidentifier –mybasic.
Code 32 shows the syntax of USER_PUNCH:••••

Line 1: the number links the user_punch to the corresponding SELECTED_OUTPUT (for the
output file) and SELECTED_GRAPH (variables that can be plotted).
Line 2: invokes the alternative scripting languageLine 4 and Line 7: indicate the start and end of the scripting blockLine 6: the punch() 9 command indicates a variable that has to be written to the output file.

USER_PUNCH <:int:number>-headings <:text:heading>-mybasic-start<:txt:mybasic instruction>
punch(<:txt:variable name>)-end
Code 32 – Syntax of USER_PUNCH
More information on the keyword USER_PUNCH is available on the online manual page.
In the graphical user interface of HPGeochemistry, USER_PUNCH is found at
HP Keywords -> Output -> USER_PUNCH
Some script functions (MyBasic form- are given in Scripting.
SELECTED_GRAPH
The keyword SELECTED_GRAPH 10 defines one or more plots that are made based on output data thatis created with SELECTED_OUTPUT and USER_PUNCH.
Code 33 shows the syntax of SELECTED_GRAPH. The definition of a single figure is shown here:
•
•

•
•

•

Line 1: The SELECTED_GRAPH definition is linked to a specific SELECTED_OUTPUT and
USER_PUNCH via number.
Line 3: Different graphs can be defined within a single SELECTED_GRAPH. The different figures
are delineated by –fig_nr. When the same number is used a second time, additional
information is added to the existing figure.
Line 4: defines the name of the graphLine 5 and 6: The current implementation supports X-Y plots in which Y-variables are plottedversus X-variable(s) at the first or secondary Y-axis or Y-X plots in which X-variables are plottedversus Y-variable(s) at the first of secondary X-axis. The default is the X-Y plot. If we consider asingle X-variable, then we define a single –fig_x and one or more –fig_y. Code 34 illustrates
how these variables are defined (see below).
Line 7 and 8: defines the labels of respectively the X and Y axis

9 In the standard scripting language, the command is
10

punch.
SELECTED_GRAPH replaces the keyword USER_GRAPH in the standard Phreeqc – the latter cannot be usedin HPGeoChemistry or HPx.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

74
ISC: PublicModelling reaction paths•

Line 9 and 11: delineate a region for Gnuplot commands. This allows for further customizing thelayout of the graph.

SELECTED_GRAPH <:int:number>-fig_nr <:int:nr>-fig_name "Name"-fig_x <:txt:var>-fig_y {<:txt:var>}-fig_xt "X label"-fig_yt "Y label"-fig_start#gnuplot commands-fig_end
Code 33 – Syntax of SELECTED_GRAPH – XY Plot
Code 34 illustrates the definition of 3 plots based on a single output file defined with
SELECTED_OUTPUT. As can be seen for –fig_x and –fig_y, a single variable requires the specificationof the type of variable using the symbol &, and the name of the variable. The type of variable is required
because it determines the prefix that is used in the output files for a given type of variable. For example,
the molality and (log)activity of species CO2 (-molalities and –activities in SELECTED_OUTPUT)
are printed in the output file with the prefix m_ (m_CO2) and la_ (la_CO2), respectively.
•••Line 12, 17, 22: general variables have the type &none (e.g., pH, pe, charge balance, water,
temperature etc.). When the type of the variable is not defined (Line 22), the &none type isassumed.
Line 13, 18: here, we specifically ask for the total concentration of Ca (Line 13) via &tot and forthe saturation index of CO2(g) (Line 18) via &siLine 23: here we ask for the molality (&mol) of three species. Several Y-variables can be definedwithin a single –fig_y statement (here 3). The variable type &mol is used for all Y-variables; if
a variable type is not given in a –fig_y statement, the previously defined variable type is used 11.

11 If you want to change from e.g. &mol to a general variable, &none should be defined. For example, &mol CO2

HCO3- &none temp.

PHREEQC-HPx Course, Prague, September 3-4, 2020

SCK CEN/39560003

75
ISC: PublicSession 2
SELECTED_OUTPUT 1-file res01.out-pH true-totals Ca-saturation_indices CO2(g)-molalities CO2 HCO3- CO3-2
SELECTED_GRAPH 1-fig_nr 1-fig_name "Solubility"-fig_x &none pH-fig_y &tot Ca-fig_nr 2-fig_name "CO2"-fig_x &none pH-fig_y&si CO2(g)-fig_nr 3-fig_name "C Speciation"-fig_x pH-fig_y &mol CO2 HCO3- CO3-2
Code 34 – Definition of three plots using SELECTED_OUTPUT
The most convenient way of defining the variables of a plot is via the help tree HP Keywords (Figure40):
HP Keywords -> Output -> SELECTED_GRAPH
To include the statements required for a XY plot (Code 33), use:
HP Keywords -> Output -> SELECTED_GRAPH -> Templates -> XY Plot
To include variables with the specification of the type, use (if components are selected, list of species
etc. are available as well):
HP Keywords -> Output -> SELECTED_GRAPH -> Variables
Curve options
After each curve (i.e. a Y-variable for a X-Y plot), specific options for the curve can be added between
angle brackets <> . Code 35 shows the syntax for the most important curve options for formatting a
curve:•••••

k: the label that appears in the legend
vk: bool to indicate if the label is present in the legend
style: to indicate if the line is plotted as a line, points, a combination of lines and points or as
dots
ls: an integer referring to a specific line definition (defined with gnuplot commands)
dt: an integer referring to a specific dash pattern (defined with gnuplot commands)< k <:str:key> vk <:bool:key visible> style<:txt:line style> ls <:int:number of line style>dt <:int:number of dash type> >Code 35 – Curve options for formatting a given curve.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

76
ISC: PublicModelling reaction paths

Figure 40 – Help tree for SELECTED_GRAPH

To include curve options, use:
HP Keywords -> Output -> SELECTED_GRAPH -> Curve options
SELECTED_GRAPH has more possibilities to combine different data sources, specify layout etc. By asmart writing of results to selected output files, a large variety of graphs can be generated. Some aspects
will be illustrated further in the course; the manual contains also some other examples.

Hands-onCreate the simulation case, database and components
Open the Project Manager
Select the Course project
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

77
ISC: PublicSession 2
Click Create New CasebuttonGive Name – “Example RE_a 02”
Give Description (optional) – “Solubility of Fe(OH)3 as a function of pH”
Click OK button
Select the simulation case
Click Select button
Go to the Database Panel
Click Components button
Select relevant components: F, Fe, Na
Click OK button
Go to the Additions to Thermodynamic database Panel
Add the phase Fix_H under PHASES.

Step 1: Calculate solubility of Fe(OH)3 at a pH of 2
To construct the input file in such a way that it can be easily reused for similar problems, we will include
different other input files into the main input file. Each of these input file cover part of the problem.
Input files are easily being included using the INCLUDE$ statement followed by the (full) file name. For
this first part, we define two additional input files:••

SolutionDefinition.in : to define the solution composition
ReactionDefinition.in: to define the reaction
Click Create new input file
Click Save As File
SolutionDefinition.in

button to open a new tab with a text editor(toolbar input work space) to save the file in project directory as

Define solution at a pH of 1.9 and with the background NaF concentration and some initial Fe(3)(0.1 mol/kgw)
Click Save FileClick Create new input file
Click Save As Filebuttonto save file as ReactionDefinition.in

Add USE solution 1
Add EQUILIBRIUM_PHASES 1 to fix the pH at 2, to have oxidizing conditions equal toatmospheric conditions and to have equilibrium with Fe(OH)3(a)Choose NaOH as the alternative formula for Fix_H
Add also –force_equality after Fix_PH
Click Save File
Go to Geochemistry Panel
Include the SolutionDefinition.in file
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

78
ISC: PublicModelling reaction paths
With the Help tree HPKeywords -> Scripting -> Include$ -> Select File –From Case, a file browser dialogue windows allows for selecting the correct file
Include END
Include the ReactionDefinition.in fileRun the case ()

Step 2: Calculate solubility of Fe(OH)3 at pH 2, 3 and 4
Using GLOBAL_VARIABLE and COUNTER, solubility and speciation calculations will be done at different
pH values. These variables are then used as inline variables in the input file.
Code 36 illustrates how GLOBAL_VARIABLE is used as an inline variable. The GLOBAL_VARIABLE
pHStart is defined in Line 1 and used in Line 6 of the solution keyword. The phreeqc.out shows that thefile SolutionDefinition.in is included (line 40) and that @pHStart is replaced with a value of 2 (line 44).
Tab Geochemistry
GLOBAL_VARIABLE pHStart 2"lower value of pH"

INCLUDE$ SolutionDefinition.in
Tab SolutionDefinition.in
SOLUTION 1
pH@pHStart
units
mol/kgw
Fe(3)
1e-1
Na
1e-1
F
1e-1
Excerpt from phreeqc.out

Code 36 – Illustration of using GLOBAL_VARIABLE as an inline variable.

Code 37 illustrates how COUNTER works. Line 3 defines the pHIterator counter with an initial value
defined by the GLOBAL_VARIABLE pHStart. By default, every time pHIterator is used as in inline
variable, its value will increment by 1. pHIterator defines the target pH of the solution, i.e. thesaturation index of the phase Fix_H+ (see line 13). We include three times ReactionDefinition.in to havea calculation at a pH 2 (pHIterator = pHStart), pH 3 (pHIterator = pHStart + 1) and pH 4
(pHIterator = pHStart + 2), as illustrated in phreeqc.out on line 147, 266 and 385, respectively.

Tab Geochemistry
GLOBAL_VARIABLE pHStart 2"lower value of pH"

COUNTER pHIterator @pHStart
INCLUDE$ SolutionDefinition.in
END
INCLUDE$ ReactionDefinition.in
INCLUDE$ ReactionDefinition.in
INCLUDE$ ReactionDefinition.in
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

79
ISC: PublicSession 2
Tab ReactionDefinition.in
USE solution 1
EQUILIBRIUM_PHASES 1
Fix_H+ -@pHIterator NaOH 10
END
Excerpt from phreeqc.out

Code 37 – Illustration of using COUNTER as an inline variable.

We define global variables via the Global Definitions Panel.
Go to Global Definition Panel
Select Item Global Variables
Functionality Edit
Select Variables in the tree.
In the Record tab, give the name and the value

Validate the record by selecting the functionality Values.
Select tab SolutionDefintion.in
Express pH as an inline variable (@pHStart) – see Code 36.
Select Geochemistry Panel
Add COUNTER pHIterator
Select Tab ReactionDefinition.in
Express saturation index of Fix_H+ as inline variable (@pHIterator)
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

80
ISC: PublicModelling reaction paths
You need to have the – sign before @pHIterator to keep the negative saturation index.
Fix_H -@pHIterator NaOH 10
Add END at the end
Select Geochemistry Panel
Copy INCLUDE$ ReactionDefinition.in two times
Run the case ()

Solubility and speciation are available in the phreeqc.out file.

Step 3: Loop over pH 2 to 12
To produce a nice graph, we want to have calculations at a large number of pH values between 2 and12, let say 100. We will overwrite the default incremental step of 1 by (12-2)/100 using COUNTERINCR
and the inline expression option. The REPEAT statement helps to define the INCLUDE$
ReactionDefinition.in 100 times.
Code 38 illustrates these different aspects. Line 9 defines the alternative value for COUNTER
pHIterator by using the first name and then the value. Note that we could use the statement:
COUNTERINCR pHIterator 0.1
Alternatively, we introduce the inline expression operator @<> with the expression between the angle
brackets. Note that other inline variables can be used in that expression. Based on two additional
GLOBAL_VARIABLES (line 2-3), the incremental step is calculated.
REPEAT in line 11 defines to create a file (slidingpH.in) by repeating @<@pHSteps+1> (=101) times
the INCLUDE$ statement for ReactionDefinition.in. By including slidingpH.in in our input file,
ReactionDefinition.in will be executed 101 times.

Geochemistry
GLOBAL_VARIABLE pHStart 2"lower value of pH"#define pH range and resolution
GLOBAL_VARIABLE pHStart 2"lower value of pH"
GLOBAL_VARIABLE pHEnd 12"upper value of pH"
GLOBAL_VARIABLE pHSteps 100 "resolution"#END USER INPUT#make sliding pH loop
COUNTER pHIterator @pHStart
COUNTERINCR
pHIterator=@<(@pHEnd@pHStart)/@pHSteps>
REPEAT
include"ReactionDefinition.in""slidingpH.in" @<@pHSteps+1>
Excerpt from phreeqc.out

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

81
ISC: PublicSession 2
Excerpt from slidingpH.in

Code 38 – Illustration of inline expressions, COUNTERINCR and REPEAT

Instead of including the commands in the Geochemistry Panel, we use the Global Definition
panel
Go to Global Definitions panel
Go to Global Variables item
Select Edit
Add two GLOBAL_VARIABLES to define upper pH value and the number of steps
Validate by Values button
Go to Transformations item
Add a variable using
Define pHIncrement
Add formula in bottom panel(@pHEnd - @pHStart) / @pHSteps
Add description
Go to Geochemistry Panel
Add COUNTERINCR
COUNTERINCR pHIterator @pHIncrement
Add REPEAT to include the ReactionDefinition.in files
Delete the INCLUDE$ ReactionDefinition.in statements
Add INCLUDE$ slidingpH.in
Run the case ()

Explore phreeqc.out to see the sequence of calculations.
Step 4 : Controlling output and plotting
SELECTED_OUTPUT, USER_PUNCH and SELECTED_GRAPH are used to select relevant output results.

To keep the input file clear, we also make an additional input file.
Click Create new input file
Click Save As Fileto save file in project directory as Plot.in

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

button to open a new tab with a text editor

82
ISC: PublicModelling reaction paths
Add SELECTED_OUTPUT 1
Choose file nameSet pH trueF-

Add -molalities for the hydrolysis species of Fe3+ and the complexes between Fe3+ and

Add USER_PUNCH 1
In this example, USER_PUNCH is used to recalculate the total concentration to log10 concentrations. The–mybasic script language is used.
Add –headings logFe
Add -mybasic
Add –start
Add the command punch(log10(TOT("Fe")))Add –end
We create two graphs. The first one is the log-transformed total Fe concentration as a function of pH.
The second are the molalities of the different Fe species.
Add SELECTED_GRAPH 1
Add a first graph via Help Tree HPKeywords -> Output -> SELECTED_GRAPH ->TEMPLATES -> XY plot
Adapt the identifiers to create the graph (consider lines 8-14 as extra):–fig_nr 1-fig_name "Solubility Fe(OH)3"-fig_x pH-fig_y logFe < vk false ls 1 >-fig_xt "pH"-fig_yt "log_{10} Fe_T (mol/kgw)"–fig_start
set xrange [@pHStart:@pHEnd]
set yrange [-8:0]
set ytics 2
set title "Solubility of Fe(OH)_3"
set style line 1 lc "blue" lw 3
–fig_endLine 4 refers to the variables that are punched in USER_PUNCH as defined in –headings
Gnuplot interprets the headings as case-sensitive names
The text between the angle brackets indicate line options for the variable logFe; vk false
means that this line is not shown in the legend of the graph; ls indicates the plotting styleof the line.

For the key, we use enhanced text formats in gnuplot{} delineates a bracketed clause to which the enhanced text options apply_ indicates subscript
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

83
ISC: PublicSession 2^indicates superscript
More information on enhanced text format in GNUPLOT is available in the topic “Enhanced
text format” in the GNUPLOT manual or online. Note that some symbols may interfere with
the HPGeochemistry and HPx file reading (e.g. “ or @).
Between line 8 (-fig_start) and line 14 (-fig_end), gnuplot commands are given to
improve the quality of the plot 12. Explanation of these commands are given in the gnuplot

help file or in examples on the gnuplot website - demos. Note that inline variables can beused; in this particular example, the x-axes boundaries are linked to the start and end valueof the pH 13. Secondly, the plotting style of the line defined in line 4 is defined in line 13 asa blue line with a thickness of 3.
(Extra) Add the second graph

Adapt the identifiers to create the graph – the help tree can be used to insert the lines 4-11
–fig_nr 2-fig_name "Speciation Fe+3"-fig_x pH-fig_y &mol Fe+3 < k "Fe^{+3}" ls 1 dt 1 >-fig_y &mol FeOH+2 < k "FeOH^{+2}" ls 2 dt 2>-fig_y &mol Fe(OH)2+ < k "Fe(OH)@_2^{+}" ls 3 dt
1 >-fig_y &mol Fe(OH)3 < k "Fe(OH)_3" ls 4 dt 1 >-fig_y &mol Fe(OH)4- < k "Fe(OH)@_4^{-}" ls 5 dt
1>
-fig_y &mol
FeF3 < k "FeF_3" ls 1 dt 2>-fig_y &mol
FeF2+ < k "FeF@_2^+" ls 2 dt 2>-fig_y &mol
FeF+2 < k "FeF^{+2}" ls 3 dt 3>-fig_xt "pH"-fig_yt "Molality (mol/kgw)"–fig_start
set xrange [@pHStart:@pHEnd]
set logscale y
set ytics format "10^{%T}"
set yrange[1E-8:1E-1]
set title "Speciation Fe^{+3}"
set key inside#define line styles (see gnuplot)
set for [i = 1:5] style line i lw 3
–fig_endLines 4-11 defines the different lines on the plot, &mol indicates that we want to plot
molalities. The indicator k between the angle brackets indicates the legend for a specificline. We use 5 different ls (1-5) that differ in colour. The lines of the first 5 variables areplotted with dash type dt 1 (full line) ; the next three lines with dasht type dt 2.
Line 17 defines that the y-scale is logarithm
Line 18 defines the format of the labels on the y-axis, here as 10X.
In line 24, the style of the lines are defined with a width of 3.

12 Note that HPGeochemistry also allows for post-processing the graph using the same set of gnuplot commands.
Thus, rather than defining these commands in the input (prior to plot generation), they could be also defined in amore interactive way after running the simulation case and generating of the plot.
13 Note that the inline variables are not available when post-processing the figure.
D. Jacques (djacques@sckcen.be)84
SCK CEN/39560003
ISC: PublicModelling reaction paths
Select Geochemistry Panel
Before the include statement for slidingpH.in:
Add INCLUDE$ plot.in
Add PRINT
Add –reset false
The statements PRINT and –reset false suppress printing all the details in the phreeqc.out file.

Output
The resulting graphs are shown in Figure 41 and Figure 42.

Figure 41 – Solubility of Fe(OH)3 as a function of pH in a 0.1 M NaF background solution.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

85
ISC: PublicSession 2

Figure 42 – Speciation of Fe(III) as a function of pH in a 0.1 M NaF background solution.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

86
ISC: PublicSession 3 – Reactive transport modelling

Session 3

Learning outcomes
Understand the relevance and power of reactive transport model
Implement simple reactive transport models via the HYDRUS-1D graphical user interface
Explore time series and profiles generated with a reactive transport model
Set-up a reactive transport problem with equilibrium dissolution and precipitation reactions

Reaction transport modelling
Coupled Reactive Transport

Leterme et al. (2014) and Batlle-Aguilar et al. (2011)

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

88
ISC: PublicReactive transport modelling

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

89
ISC: PublicSession 3

Steefel et al. (2005)

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

90
ISC: PublicReactive transport modelling

Examples of reactive transport model studies

Jansen (2003)

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

91
ISC: PublicSession 3

Seuntjens et al. (2002)
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

92
ISC: PublicReactive transport modelling

Jacques et al. (2008b)
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

93
ISC: PublicSession 3

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

94
ISC: PublicReactive transport modelling
Jia et al. (2019)

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

95
ISC: PublicSession 3

Jacques et al. (2018)

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

96
ISC: PublicReactive transport modelling
Makselon et al. (2017)

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

97
ISC: PublicSession 3

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

98
ISC: PublicReactive transport modelling

Vandenhove et al. (2007), Vandenhove et al. (2011), Payne et al. (1996)

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

99
ISC: PublicSession 3

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

100
ISC: PublicReactive transport modelling

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

101
ISC: PublicSession 3

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

102
ISC: PublicReactive transport modelling

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

103
ISC: PublicSession 3

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

104
ISC: PublicReactive transport modelling

Mayer et al. (2015)

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

105
ISC: PublicSession 3
Xie

etD. Jacques (djacques@sckcen.be)
SCK CEN/39560003

al.

(2015)

106
ISC: PublicReactive transport modelling

Mathematical Model

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

107
ISC: PublicSession 3

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

108
ISC: PublicReactive transport modelling

Reactive transport modelling – The HPx model

Jacques et al. (2006), Jacques et al. (2008a), Jacques et al. (2018)

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

109
ISC: PublicSession 3

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

110
ISC: PublicReactive transport modelling

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

111
ISC: PublicSession 3

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

112
ISC: PublicReactive transport modelling

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

113
ISC: PublicSession 3

Reactive transport modelling (HP1) – Equilibrium mineral
precipitation and dissolution
Dissolution and precipitation of gypsum and calcite
Problem definition
Sulphate-free water is infiltrated in a 50-cm long uniform soil column under steady-state saturated flow
conditions. The reactive minerals in the soil column are calcite (CaCO3) and gypsum (CaSO4.2H2O), both-2with a concentration of 2.176 x 10 mmol / kg soil.

The physical properties of the soil column are a porosity of 0.35, a saturated hydraulic conductivity of 10
cm / day, bulk density is 1.8 g/cm³ and a dispersivity of 1 cm.
The input solution contains 1 mM CaCl2 and is in equilibrium with atmospheric partial pressure of oxygenand carbon dioxide. The soil solution is in equilibrium with the reactive minerals and with atmosphericpartial pressure of oxygen. The latter is to have only the oxidized components of S and C.
Calculate the dissolution fronts of calcite and gypsum over a period of 2.5 days.

Conceptual model
The problem involves following components:
•
•
•
•
•
•
•

Equilibrium reactions
Two phases (aqueous, solids)
Well-mixed
Homogeneous domain
1D flow
1D advective-dispersive solute transport
Constant boundary condition

Figure 43 shows a schematic representation of the conceptual model.

Figure 43 – Schematic representation of a 1D reactive transport problem involving equilibrium mineral dissolution
and precipitation reactions.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

114
ISC: PublicReactive transport modelling

Hands-on
Step 1 Create project
Project Manager
Button: "New"Name: "T-Em-01"
Description: "Dissolution of calcite and gypsum in soil profile"
Button: "OK"Step 2: Define the problem and main settings
Main Processes
Heading: Dissolution of calcite and gypsum in soil profile Check: "Water Flow"
Check: "Solute Transport" Select: "HP1 (PHREEQC)"
For a new project, the box HP1 version 2020 is checked. When checked:•

The most recent version of HP1 is used (version 3.01.003)•

The new GUI for defining the problem is used

For back-compatibility with previous projects, an older version of HP1 is used when unchecked,
although most previous projects should run with these new versions.
Button: "Next"
Geometry Information
Depth of the Soil Profile: 50 cm
Button: "next"
Time Information
Final Time: 2.5 days Maximum Time Step: 0.05
Button: "next"
Print Information
Unselect: T-Level information Select: Print at Regular Time Interval Time Interval: 0.025
Print Times: Number of Print times: 5 Button: "next"
HP1 – Print and Punch Controls
In most cases, the defaults in the HP1 – Print and punch Controls dialogue window are valid. The input
is divided in four panels:•

Punch Time and Locations: If controlled by HYDRUS, output times and locations arecontrolled by the graphical user interface of HYDRUS. Output times for observation points orsummed values are defined in the Print Information dialogue window, output locations forobservation points in the Soil Profile – Graphical Editor. Output times for profiles were also
defined in the Print Information dialogue window. When you want to use the HYDRUS1D
graphical user interface to look at profiles or time series, the two options (Observ. Nodes
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

115
ISC: PublicSession 3•••

Printed to Different Files, Mobile and Immobile Cells in different
Files) should be checked. If not, no plots will be generated. The checkbox Make Gnuplot
Templates is deprecated – when checked, Gnuplot templates are written and saved in the
project directory. Alternatively, one can select output times independently from the definitions
by defining the nodes (punch cell) and the frequency (punch frequency).
Print Options: To indicate when and for which cells a complete chemical description is written to
the phreeqc.out file. As this printing takes some time, it is disabled by default.
Selected Ouput: This option will give an output file with selected output for geochemical
calculations prior to the transport simulations.
PHREEQC dump: This option is to select when a dump file is created, i.e. a complete chemical
state that can be used to restart a simulation. We refer to the PHREEQC manual for more
information on the dump file.

Figure 44- The HP1 - Print and Punch Controls dialogue window

Button: "Next"
Step 3: Water flow
Water Flow – Iteration Criteria
Button: "Next"
Water Flow – Soil Hydraulic Model
Button: "Next"
Water Flow – Soil Hydraulic Parameters
Qs: 0.35
Ks: 10
Button: "Next"
Water Flow – Boundary Conditions
Upper Boundary Condition: Constant Pressure Head
Lower Boundary Condition: Constant Pressure Head

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

116
ISC: PublicReactive transport modelling
Step 4: Solute transport
The solute transport problem is defined by the next dialogue windows.
Solute Transport – General Information
Stability Criteria: 0.25
Button: "Next"

Solute Transport – HP Input
The HP Input dialogue window for transport is similar to the input of HPGeochemistry.
Database
When you open the HP Input dialogue window for the first project within a Project Group (ProjectManager in the HYDRUS 1D graphical user interface), the Default Database for Project Group dialogue
window (Figure 45) appears to define the default database for new projects in the Project Group. There
are different options:•••

When “Use Default” and “Use General Default” are checked (Figure 45), the general
default database is used for new projects, i.e. phreeqc.dat from the database group Phreeqc
3.6.2.
When “Use Default” is checked, but “Use General Default” is unchecked, a databaseselected by the user via the “Select” button (right column) will be used. The “Select” button
opens the Database manager.
When both are unchecked, a new project opens without a default database. The Databasemanager appears to define the database.

The database can be changed via the Database Panel.
Note that default databases for Project Groups can be changed via Preferences and Settings () -> Settings -> Databases -> “Set” button after “Default databases for HP1 projects(Defined in the HP1 GUI)”.

Figure 45 – The Default Database for Project Group dialogue window.

Button “Close” – we use the general default database for new projects during the course.
Reactive Transport Panel
This tab is used to define standard (aqueous components) and other more advanced options.
Select the item Selection of Components
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

117
ISC: PublicSession 3
Based on the database, lists of aqueous components (first column) and gases (second column) are given.
Note that gases are based on the phases defined in PHASES in the database; when the phase namecontains (g), it is considered as a gas.
The user select the components that need to be transported. In this example, we have a CaCO3-CaSO4
problem and a CaCl2 solution. Thus, the components are Ca, Cl, C(4), S(6). Note that two additional
components are always transported:••

Total_O: this (hypothetical) component contains all O that is not in H2O.
Total_H: this (hypothetical) component contains all H that is not in H2O
Selected Ca, Cl, C(4) and S(6)Note that components can be added that are not in the database. There are two options:••

You define the required database information via the Add information button in the
database Panel. An additional column will appear from which you can select the additional
components for the reactive transport.
You define the new components in the column “Components from Text Input”. You needto provide the required thermodynamic databases in the input text files (e.g., in the tab“Additions to database”).

Geochemical model
The additional tabs are used to prepare the information related to the geochemical model (PHREEQC
input):
•
•

•
•

•

Global definitions: same functionality as in HPGeoChemistry.
Additions to the thermodynamic database: used to define additional species,
reactions and rates (SOLUTION_MASTER_SPECIES, SOLUTION_SPECIES, PHASES,EXCHANGE_MASTER_SPECIES,
EXCHANGE_SPECIES,
SURFACE_MASTER_SPECIES,
SURFACE_SPECIES, RATES)
Definition of initial solutions: used to define the initial and boundary solution
compositions (SOLUTION)
Geochemical
model: used to define the geochemical model (EXCHANGE,
EQUILIBRIUM_PHASES, SURFACE_SPECIES, SOLID_SOLUTIONS, KINETICS) and somenumerical information if needed,
Additional output: used to define the output (SELECTED_OUTPUT, USER_PUNCH)

Only the most relevant or typical PHREEQC keywords are listed; however, these and other keywords canbe used in the different text boxes. It is recommended that the user is then familiar with the phreeqc.in
file that is generated from these four text boxes and additional information from the graphical userinterface.
In our particular example, we only need to define two solutions (1 boundary solution, 1 initial solution).
However, we will also ask for the concentration of the components in SELECTED_OUTPUT although thisinformation is also given in the standard output of HYDRUS.
The link between the solutions defined in the geochemical model (PHREEQC) and the initial andboundary conditions to be defined for the solute transport problem is by the solution number. It is
recommended not to use a number smaller than the number of nodes in the transport problem to avoid
any possible problem of conflict. For an one-dimensional problem defined in HYDRUS1D, the maximum
number of nodes is 1000. We recommend therefore to use always solution numbers larger than 1000.
In the exercise (or as a general rule of thumb), following reasoning is used:
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

118
ISC: PublicReactive transport modelling••••Numbers between 1001-2000: Initial condition of mobile cells
Numbers between 2001-3000: Initial condition of immobile cells
Numbers between 3001-4000: Upper boundary condition
Numbers between 4001-5000: Lower boundary condition
Select Solution Definitions Panel
Add SOLUTION 1001 : A solution in equilibrium with gypsum and calcite
SOLUTION 1001
EQUILIBRIUM_PHASES 1001
gypsumcalciteO2(g) -0.68
SAVE solution 1001
END
Add SOLUTION 3001 : upper boundary solution (right side) with CaCl2
SOLUTION 3001
pH 7 charge
Cl 2
Ca 1
O(0) 1 O2(g) -0.68
C(4) 1 CO2(g) -3.5

Define for each node the geochemical model. Note that the initial amount of a mineral must be definedas mol/1000 cm³ soil.
There are two options to link geochemical entities (EXCHANGE, EQUILIBRIUM_PHASES, SURFACE,
SOLID_SOLUTION; KINETICS) to nodes:••

Link the entity to the specific node number explicitly (e.g. EQUILIBRIUM_PHASES 1 or for arange EQUILIBRIUM_PHASES 51-101)
Link the entity to the nods by linking it to a material number using the identifier -material.
The spatial distribution of a material is defined in the Profile information form.

To avoid conflict with node numbers, it is recommended to use values larger than the number of nodes.
The same guidelines as for the solutions can be used (numbers between 1001-2000 for mobile cells,
between 2001-3000 for immobile cells)
Select Geochemistry Panel
EQUILIBRIUM_PHASES 1001
-material 1
gypsum 0 3.9E-5
calcite 0 3.9E-5
O2(g) -0.68
Select Output Panel
SELECTED_OUTPUT-totals Ca Mg Cl S C-equilibrium_phases gypsum calcite
Button OK or Next on the main tool strip

Solute Transport – Solute Transport Parameters
Disp: 1 cm
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

119
ISC: PublicSession 3
Button: "Next"
Solute Transport – Boundary conditions
Upper Boundary Condition

Bound. Cond. 3001

Step 5 : Geometry and initial conditions
Click OK
Click Yes
Soil Profile – Graphical Editor
Menu: Conditions -> Initial Conditions -> Pressure Head Button: "Edit Condition"
Select All Top Value : 0
Menu: Conditions -> Observation Points Button: "Insert"
Insert 5 observation nodes, one every 10 cm
Menu: File -> Save Data
Menu: File –> Exit
Soil Profile – Summary
Button: "Next"

Run Application
OutputThe standard HYDRUS-output can be viewed through the Post-processing unit. HP1 Output allows forviewing additional geochemical information. The specific dialogue window for HP1 output opens via
HP1 output and is similar to the output workspace of HP1Geochemistry.
In the phreeqc.out work space, there is a part indicated by “Overview model”. You can check how the
information from PHREEQC and the GUI are linked to the transport model.
Check the section “Overview model” in the phreeqc.out
The graph workspace lists the graphs that are created based on the information in the HYDRUS GUI, i.e.
time series at selected nodes, and profiles at selected times for the variables defined in the HP1 print
information dialogue window and with SELECTED_OUTPUT (and USER_PUNCH) definitions.
Explore the graphs in the output

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

120
ISC: PublicReactive transport modelling

Figure 46 - Profiles of C concentrations (Example T-Em-01)

Hereafter are some examples how to adapt the figures using the Gnuplot commands:
The numbers on the x-axis for the profiles of C (Figure 46) are obtained via
Select Profiles -> Total Aqueous Concentration -> C
Select Edit graphInsert “set xtics 1E-5 format "%g"”See “format specifiers” in gnuplot help file. Alternatively, use:“set xtics 1E-5 format "%0.0t x 10^{^%T}"”

Transport of Cd
Problem Definition
The transport of a heavy metal through the soil column is investigated. Under high pH conditions, Cd
precipitates as CdCO3 or otavite. However, due to changing geochemical conditions in the soil column,
solubility of Cd is changing. This example studies Cd mobility through the soil column.
The physical and geochemical set up as in Part I is used as the basis. The infiltrating water is contaminatedwith a small amount of Cd (1 x 10-6 M Cd). Otavite is added to the geochemical model with a small initialamount (1 x 10-10 mol/1000 cm³ of soil).
Calculate the percentage of Cd in the aqueous phase.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

121
ISC: PublicSession 3
Hint: The basic statement SYS("Cd") gives the total moles of Cd in the system, whereas the statement
TOT("Cd") gives the total concentration of Cd in the aqueous phase. The total amount of water (in kg)
is obtained by the basic statement TOT("water").

Geochemical system••••

Will otavite precipitate in the gypsum-calcite system for the given Cd contamination?
Will otavite precipitate in the calcite system after gypsum is completely dissolved for the given
Cd contamination?
Calculate the aqueous concentration of Cd in equilibrium with otavite in both systems. What is
the reason for the difference?
Will otavite precipitate in the incoming solution?

HPGeochemistry can be used to solve these questions.

Conceptual model
The conceptual model is similar to that of the previous model (see page 114 and Figure 43).

Hands-on
Step 1 – Create projectOpen project T-Em-01

Project Manager
Button: "Copy"
Name: "T-Em-02"
Description: "Dissolution of calcite and gypsum in soil profile – Cd Transport"
Button: "OK"Main Processes
Heading: "Dissolution of calcite and gypsum in soil profile – Cd Transport"
Button: “OK”

Step 2 – Adapt geochemical model
Database and Reactive Transport – Selection of Components
Select: Solute Transport – HP1 Input Data
Select: Reactive Transport Tab – Selection of Components item
Select: Cd in the Components from Database

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

122
ISC: PublicReactive transport modelling
Geochemical model and Output
Solutions
Select Solution Definitions Panel
Add Cd in solution 3001
Geochemical model
Select GeoChemistry Panel
Add otavite in EQUILIBRIUM_PHASES
OutputSelect Output Panel
Add otavite to SELECTED_OUTPUT –equilibrium_phases
Add USER_PUNCH data block to calculate the percentage of Cd in solution
For each of the variables, it is possible to define an axis label (plot_titles) and a unit
(plot_units)
USER_PUNCH-headings CdPer-plot_titles "Percentage Cd in solution"-plot_units "%"-mybasic-startif (sys("Cd") > 0) thenperCd = 100 * tot("Cd") * tot ("water") / sys("Cd")
else
perCd = 0
endif
punch(percd)-end
Button OK or Next on the main tool strip

Run Application
OutputLook at the profiles of Cd, Otavite and percentage of Cd in solution (Figure 47).
The small amount of octavite initially present in the soil column dissolves completely (100% of Cd is insolution at day 0) (see Starting questions). When only calcite is present (i.e., after complete gypsum
dissolution), otavite precipitates and almost 90% of Cd is in the solid phase. After calcite is dissolved, thesolubility of otavite increases and all Cd is mobilized again.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

123
ISC: PublicSession 3

Figure 47 – Profiles of percentage of Cd in solution (Example T-Em-02)

Effect of a higher Cl concentration on Cd mobility
Problem Definition
Aqueous components which form strong complexes with Cd will enhance the mobility of Cd.
In this example, the composition of the inflowing water is changed after 1 day to a solution
with a higher CaCl2 concentration (1 x10-2 M).
Evaluate the effect of the higher Cl concentration on Cd mobility.

Conceptual model
The conceptual model is similar to that of the previous model (see page 114 and Figure 43), except thata time-varying upper boundary condition for the solution concentration is needed (Figure 48).

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

124
ISC: PublicReactive transport modelling

Figure 48 – Schematic representation of a 1D reactive transport problem involving equilibrium mineral dissolution
and precipitation reactions with a time-varying upper solute concentration boundary condition.

Hands-on
Step 1 – Create ProjectOpen project T-Em-02

Project Manager
Button: "Copy"
Name: "T-Em-03"
Description: "Dissolution of calcite and gypsum in soil profile – Cd Transport; effect of Cl"
Button: "OK"Main Processes
Heading: "Dissolution of calcite and gypsum in soil profile – Cd Transport; effect of Cl"
Button: “OK”

Step 2 – Time information
Time Information
Check: “Time Variable Boundary Conditions”
Number of Time Variable Boundary conditions: 2
Button: “OK”

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

125
ISC: PublicSession 3
Time Variable Boundary Conditions
Time: 1, 2.5
cTop: 3001, 3002
cBot:4001, 4001
Button: “OK”
Step 3 – Geochemical model and output
Select Solution Definitions Panel
Add solution 3002: Boundary solution with higher CaCl Concentration
Add solution 4001: Bottom boundary solution – pure water
Button OK or Next on the main tool strip
Run Application
OutputTime series of Cl for the case of low and high CaCl2 and profile data for otavite are shown in Figure 49.
At the first observation depth, breakthrough curves of Cd are quite similar. The calcite and otavite
dissolution fronts already reached this depths after 1 day, and the high Cl concentration have only asmall effect on it: the increase in Cd concentration starts faster after 1 day.
However, deeper in the soil, a release of Cd occurs already after 1.5 days at 20 cm depth in case of highCl concentrations. Note that, in case of low Cl concentration, the Cd peak does not pass the depth 20
cm after 2.5 days. In case of the high Cl concentration, the otavite dissolution front does not follow thecalcite dissolution front as is the case for the low Cl concentration front, but the Cl concentration front(Figure 50). Due to the high Cl concentrations, otavite is not formed anymore after 1.5 days.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

126
ISC: PublicReactive transport modelling

Figure 49 – Time series of Cd (top) and profiles of otavite (bottom) for transport without (left, T – SS – Em – 02)
and with (right, T – SS – Em – 03) higher Cl concentrations.

Figure 50 – Time series of Cl concentrations for T – SS – Em – 03.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

127
ISC: PublicSession 3

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

128
ISC: PublicSession 4 – Equilibria at the solid-liquidinterfaceSession 4
Many solid phases in geochemical system can sorb chemicals. This is generally called sorption but it
may refer to various processes such as adsorption (adherence of a chemical to the surface of the solid),
absorption (uptake of the chemical inside the solid), or exchange (replacement of one chemical byanother on the solid surface). There are empirical equations that describe the sorption of elements to asurface, e.g. linear isotherms, Freundlich isotherms or Langmuir isotherms. These equations, available in
HYDRUS, but they can also be implemented in PHREEQC as well (e.g. example 19 in Phreeqc manual,(Parkhurst and Appelo 2013b), see also online version). Their parameters are conditionally on thephysical-geochemical conditions and do not take into account any competitive adsorption processes.
The two major groups of models that account for physical-geochemical conditions and competitive
adsorption processes are:Ion exchange models – a charge-neutral exchange of ions between the aqueous phase andexchange sites at the surface of solid phases in which the surface has a fixed charge (e.g. fixedcharge on clay minerals). The site capacity of a surface is thus always occupied with species of the
opposite charge such that the surface is charge neutral.
Surface complexation models – these models describes the speciation of surface functiongroups by accounting for chemical bound (chemical affinity) and electrostatic effects that depends
on the surface charge. Surface complexation models differ in the way the division is made between
charged surface and bulk solutions and the charge/potential relations that are used. The four mostcommon models are the constant-capacitance model (CCM), the diffuse double layer model(DDL), the triple layer model (TLM) and the CD-MUSIC model (Goldberg et al. 2007). These models
and simpler variants (e.g. neglecting the electrostatic effects) can be implemented in PHREEQC.
In this course, we treat only single-site cation exchange models.

Learning outcomes
Basic understanding of cation exchange processes
Implementation of exchange processes using EXCHANGE
Calculating reaction paths involving cation exchange processes and cation exchange isotherms
Including cation exchange in a reactive transport model

Geochemical modelling – Cation Exchange
Mathematical world – Cation Exchange
Exchange equationsz B Az A + z A BX zB ←→ z B AX z A + z A B zB

(10)( B ) ( AX )K =( A ) ( BX )

(11)

zB

zA

zB

zA

s

z A zBzAzB

The activity of the exchange species are typically expressed as molar or equivalent fractions relative toeither the total number of exchange sites or the total number of exchangeable cations, according to oneof the following conventions
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

130
ISC: Public

Equilibria at the solid-liquid interface
Gaines-Thomas convention – activity of an exchange species equals it equivalent fraction relative
to total number of exchange sites;
Vanselow convention – activity of an exchange species equals it molar fraction relative to thetotal number of exchange sites;
Gapon convention – activity of an exchange species equals the fraction (molar or equivalent)
relative to the total number of exchangeable cations.
For the Gaines-Thomas and the Vanselow convention, the exchange equation Eq. (10) is written in termsof mole per exchanging cation:(12)→ 1 zA AX z A + 1 zb BzB
1 zA A zA + 1 zB BX zB ←The equivalent fraction βe used in the Gaines-Thomas convention is calculated as:

=
β e,I

z I nI=∑ z J nJ

ne, I
ne, I
=
∑ zJ ne,J CEC

(13)

J A=
, B ,...
J A, B ,...
=

where I, A, B denotes exchange species, z expresses charge, n expresses the concentration on theexchanger (mol/kgs), ne expresses equivalents on the exchanger (eq/kgs), and CEC is the cation exchange
capacity (eq/kgs).
The molar fraction βm used in the Vanselow convention is calculated as:

β m,I
=

nI
=
∑ nJ

ne I z I
nI=∑ ne J zJ TEC(14), B ,...
J A=
J A, B ,...
=

Where TEC is the total exchangeable cations (mol/kgs), and the other symbols as defined above.
For the Gapon convention, the exchange equation Eq. (10) is written in terms of a single exchange site(per mole of exchange site):→ A1 zA X + 1 zb BzB
1 zA A zA + B1 zB X ←

(15)

Equivalent and molar fractions are identical in this case.
For homovalent exchange, the three convention produces similar results. However, the results differs forheterovalent exchange, see section 6.3 in Appelo and Postma (2005).

Exchange isotherm
An exchange isotherm relates the equivalent fraction of a cation in solution to the equivalent fraction of
the cation on the exchanger. In a system with two cations A and B, the equivalent fraction of B in solution
EB and the equivalent version of B on the exchanger, βe,B are defined as:=EBz B mBz B mB=zA mA + zBmB zA mA,0

(16)

=
β e,B

z BnB
z BnB
=
zA nA + zBnB CEC(17)

where m expresses the molality (mol/kgw), n expresses the concentration on the exchanger (mol/kgs), z
the charge, and CEC is the cation exchange capacity (eq/kg) and mA,0 is the initial concentration of cation
A (mol/kgw). The second equality in Eq. (16) is, when starting from a solution and exchanger in the Aform, A is replaced by B. The denominator zAmA,0 then represents the solution normality (which remains
constant during the ion exchange process).
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

131
ISC: PublicSession 4
Ion exchange isotherms illustrates the preference or selectivity of an exchanger for one ion over another.
Figure 51 shows three types of selectivity for ion A over an ion B: Curve (i) shows no selectivity for ion A
over B or equal preference for ion A and B. Curve (II) shows an exchanger with preference for ion A overion B. Curve (iii) is the case for preference of ion B over ion A.

Figure 51 – Examples of isotherms for simple exchangers (from Pabalan and Bertetti (2001))

Implementing the model
Figure 52 shows the PHREEQC keywords for implementing this problem. The keywords
EXCHANGE_MASTER_SPECIES, EXCHANGE_SPECIES and EXCHANGE are introduced in this section.

Figure 52 – Implementation of the model for cation exchange processes.

Database – EXCHANGE_MASTER_SPECIES and EXCHANGE_SPECIESEXCHANGE_MASTER_SPECIESThe keywords EXCHANGE_MASTER_SPECIES and EXCHANGE_SPECIES define the exchange equations
for the thermodynamic geochemical calculations.
EXCHANGE_MASTER_SPECIES is similar to SOLUTION_MASTER_SPECIES but is used to define one ormore exchangers. The syntax and an example are given in Code 39.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

132
ISC: Public

Equilibria at the solid-liquid interface
EXCHANGE_MASTER_SPECIES<:txt:Name> <:txt:master species>
EXCHANGE_MASTER_SPECIESX XCode 39 – Syntax and example of EXCHANGE_MASTER_SPECIES
More information on the keyword EXCHANGE_MASTER_SPECIES is available on the online manual
page.
EXCHANGE_SPECIESEXCHANGE_SPECIES defines the thermodynamics of the exchange reactions. An example is given in
Code 40 using the Gaines and Thomas convention:•Line 2 – 3: For each exchanger, an identity for the master exchange species is required

EXCHANGE_SPECIES
X = Xlog_k 0;
Na+ + X- = NaX
log_k 0;K+ + X- = KX
log_k 0.7
Ca+2 + 2 X- = CaX2
log_k 0.8
Code 40 – Definition of half reactions for cation exchange reactions on the exchanger X using
EXCHANGE_SPECIES.

The reaction equations for cation exchange (Eq. (10)) are defined as half reactions (line 5, 8, 9 in Code40) and the log K value is relative to a reference half reaction with a log K of 0.0. In Code 40, the reference
half reaction is for NaX (line 5-6) which is also the case in the PHREEQC v3.6.2 databases.
The log K of the other half reactions (for K in line 8 and for Ca in line 11) are relative to the reference
half reaction. Consider the exchange reaction Na/K written according to the Gaines-Thomas convention(Eq. (12)):

LogK Na K = −0.7

(18)

( NaX ) ( K + )  ( NaX )   ( KX ) =K Na K =( KX ) ( Na + )  ( Na + )( X - )   ( K + )( X - ) 

(19)

Na + +KX ←→ NaX+K +This mass action equation is

−1

Thus

Log=
K Na K LogK Na − LogK K(20)With a LogKNa/K of -0.7 and LogKNa of 0.0 (reference half reaction), the LogKK for the half reaction (Line8) is 0.7 (Line 9).
Similarly, for the Na/Ca reaction:

Na + +0.5CaX 2 ←
→ NaX+0.5Ca 2+
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

LogK Na Ca = −0.4

(21)
133
ISC: PublicSession 4This mass action equation is( NaX ) ( Ca + )  ( NaX )   ( CaX 2 ) 
K Na K ==
0.5
2( CaX 2 ) ( Na + )  ( Na + )( X - )   ( Ca 2+ )( X - ) 
0.5

−0.5

(22)

Thus
(23)

Log=
K Na Ca LogK Na − 0.5LogK Ca

With a LogKNa/Ca of -0.4 and LogKNa of 0.0 (reference half reaction), the LogKK for the half reaction (Line11) is 0.8 (Line 12).
The databases included with PHREEQC (v3.6.2) uses the Gaines-Thomas convention for exchange.
Appelo and Parkhurst (2002) gives examples on how the other conventions can be implemented in
PHREEQC – see online version.
For each of the half reactions, activity coefficients and temperature dependency of LogK can be definedwith identifiers –gamma, -davies, -delta_h, -analytical_expression.
More information on the keyword EXCHANGE_SPECIES is available on the online manual page.
To model sigmoidal exchange isotherms, either the active fraction model (defined with –gamma, see theonline manual page and following example on Appelo’s webpage, or with the Rothmund-Kornfeld
model (Bond 1995) (see following example on Appelo’s webpage or Appelo and Parkhurst (2002) see
online version).

Input – EXCHANGEEXCHANGE defines the exchanger properties used during the geochemical calculations.
EXCHANGE 1
X 0.1-equilibrate with solution 1
Code 41 – Implicit definition of an EXCHANGE.
Code 41 shows an example of the EXCHANGE keyword:••Line 2 defines the amount of exchanger X. The amount is defined in moles in your system. The
resulting concentrations of the exchange species are reported in mol/kgw (e.g. via –molalities in SELECTED_OUTPUT). This line can be repeated for different exchangers todefine a multi-site exchange assemblage.
In this definition, the composition of the exchanger is calculated as it is in equilibrium withsolution 1 (Line 3). This calculation will be performed before any batch reaction will be
performed and thus will not change the solution composition.
This is further illustrated in Figure 53. The input is shown at the top: a solution and an exchanger
are defined; the latter being initially calculated in equilibrium with the solution. The
phreeqc.out file contains an initial solution calculation, an initial exchanger calculation and abatch reaction calculations (figures in the middle). If you look in detail to the batch reaction
calculation in phreeqc.out (not shown), you see that both the solution and exchanger
composition are the same as in the two initial calculations and is thus an unnecessary calculation.
The batch reaction calculation can be suppressed by uncommenting line 10: there will be nobatch reaction because there is no solution in the geochemical model (see bottom figure).
Input

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

134
ISC: Public

Equilibria at the solid-liquid interface

OutputWhen Line 10 is uncommented

Figure 53 – Illustration of the implicit definition of EXCHANGE (see text for explanation)

As an alternative, the initial composition of the exchanger can be defined explicitly, as shown in Code42.
EXCHANGE 1
NaX 0.01
CaX2 0.06
Code 42 – Explicit definition of an EXCHANGE.

An exchanger can be linked to equilibrium phases or kinetic reactants.
More information on the keyword EXCHANGE is available on the online manual page.

Speciation model - Equilibration of an exchanger with a solution
Learning outcomesCalculate the cation exchanger composition in equilibrium with a solution
Understand the USE and SAVE keywords

Illustration of SAVE and USE
With the SAVE and USE keywords, geochemical entities (SOLUTION, EQUILIBRIUM_PHASES,EXCHANGE, etc.) are saved after a batch reaction calculation or are used in a specific batch reaction.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

135
ISC: PublicSession 4
After a (batch) reaction calculation when equilibrium is attained, the changed composition of thegeochemical entities are not saved 14 automatically. If one wants to use the calculated composition inanother simulations during the execution (i.e., in another data block that is delineated with END), thegeochemical entity has to be saved using SAVE (Code 43).
More information on the keyword SAVE is available on the online manual page.
Initial geochemical entities or saved geochemical entities after a reaction calculations can be used insubsequent simulations using USE. USE can also specify other entities such as REACTION. Moreinformation on the keyword USE is available on the online manual page.

SAVE solution 1
SAVE <:txt:entity> <:int:number>
USE solution 1
USE <:txt:entity> <:int:number>
Code 43 – Illustration and syntax of SAVE and USE.
Example CE_ex01 illustrates USE and SAVE when defining the composition of an exchanger. The input
file is listed in Code 44.

14 With save, we mean within the memory of the simulation. Depending on the specifications in the input file, thecomposition are written to the phreeqc.out file or to selected output files.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

136
ISC: Public

Equilibria at the solid-liquid interface#Simulation 1
SOLUTION 1
Na 0.1C 1 charge
EQUILIBRIUM_PHASES 1
calciteEND#Simulation 2
EXCHANGE 1
X 0.1-equilibrate with solution 1
END#Simulation 3
SOLUTION 2
Na 0.1C 1 charge
EQUILIBRIUM_PHASES
calciteSAVE solution 2
END

2#Simulation 4
EXCHANGE 2
X 0.1-equilibrate with solution 2
END#Simulation 5
USE exchange 1
SOLUTION 3-pH 7 charge
Mg 0.1SAVE exchange 1
END#Simulation 6
USE exchange 1
USE solution 2
USE equilibrium_phases 2
END

Code 44 – Input file for Example Ce_ex01.

Simulation 1
Simulation 1 performs an initial solution calculation and a batch reaction calculation with SOLUTION 1
and EQUILIBRIUM_PHASES 1 which is clear from outline of the phreeqc.out file:

The results of the calculations are printed in the output file which shown that both the solution and theequilibrium phase assemblage have changed:

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

137
ISC: PublicSession 4

However, the chemical entities SOLUTION 1 and EQUILIBRIUM_PHASES 1 are not changed in thememory after the batch reaction calculation.
Simulation 2
Simulation 2 equilibrates EXCHANGE 1 with SOLUTION 1. Because the results of the reaction calculations
from simulation 1 are not saved to the memory, the initial calculation of SOLUTION 1 is used toequilibrate the exchanger. Consequently, the exchanger is completely covered by Na:

Simulation 3
Simulation 3 is exactly the same as simulation 1 (using number 2 instead of 1), except that the solution
after the batch reaction calculation is saved in memory as solution 2.
Simulation 4
Simulation 4 equilibrates the exchanger with solution 2. Contrary to simulation 2, it is now a solution
containing the cations Ca and Na as solution 2 is the solution after the batch reactions from simulation
3. The exchanger contains now both Ca and Na:

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

138
ISC: Public

Equilibria at the solid-liquid interface

Simulation 5
Simulation 5 defines a new solution (SOLUTION 3 containing Mg) and performs first an initial
calculations. Then, a batch reaction calculation is performed with SOLUTION 3 and EXCHANGE 1.
EXCHANGE 1 was defined in simulation 2 and equilibrated with SOLUTION 1 from simulation 1 andcontains only Na. After the batch reaction calculation, the exchanger contains Na and Mg. Because thehigh selectivity for Mg compared to Na, almost all Mg is on the exchanger:

The resulting exchange composition is saved in memory as EXCHANGE 1.
Simulation 6
In simulation 6, entities are all defined in an implicit way: EXCHANGE 1 (the Mg(-Na) exchanger obtained
after simulation 5, SOLUTION 2 after simulation 3 and EQUILIBRIUM_PHASES 2 (the initial definition).
Note that this is also reflected in the phreeqc.out file:

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

139
ISC: PublicSession 4

Reaction path – Exchange isotherm
Learning outcomesCalculate a cation exchange isotherm – one element is replaced by another
Create a multi-dataset file for plotting with Gnuplot
Plot results for different datasets in a multi-dataset fileProblem Definition
Calculate the exchange isotherm for three different background concentrations (10-3M, 10-2M, 10-1M) in
1 litre of water with an exchanger (capacity 1 mmol/l) initially completely covered with Na.

Conceptual model
The problem involves following components:
•
•
•
•
•

Equilibrium reactions
One phases (aqueous)
Surface-liquid interface – Cation exchange process
Well-mixed
Adding and removing of aqueous components

We model a cation exchange isotherm. .
Figure 54 shows a schematic representation of the key elements to calculate a cation exchange isotherm.
In essence, no solid phase is present, only a geochemical entity to represent the cation exchange
complex (EXCHANGE). The idea is to remove in a step wise manner all moles in the system of a given
component (Na) and replace it charged-balanced by another component (Ca). To calculate the isotherm,
the keyword REACTION is used.

Figure 54 – Schematic representation of the model to calculate cation exchange isotherms (left) and key elementsof the computer model.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

140
ISC: Public

Equilibria at the solid-liquid interface
The exercise asks to calculate three isotherms and to plot them in a single graph. The calculation needs
to be repeated for three times (GLOBAL_VARIABLE, INCLUDE$). There are several ways to obtaingraphs of such repetitive calculations in HPGeochemistry. In this particular example, we use the keyword
PUNCH to create a so-called multi-dataset file.

Multi-Dataset File
Gnuplot uses text files as source for the data to make a plot. The first line typically contains column
headers and the following lines data records. Such file format is created when SELECTED_OUTPUT isused to write geochemical results.
Blanc lines have a specific meaning in Gnuplot. A single blanc line indicates a discontinuity in a curve –no line will join points separated by a blanc line. A double blanc line indicates a different dataset that isplotted as a separate curve.
To create a file following the Gnuplot conventions, following statements are used:••

Printing to an output file can be supressed by the indicator –active(line 7 in Code 45).
If false, new geochemical calculations are not printed in the output file. Printing will restart when
–active is reset to true. In this way, calculations not needed for plotting (e.g. initial
calculations) can be excluded from the plot.
The keyword PUNCH 15 helps to control a multi-dataset file. Line 8 defines a label for a specificdataset that will be used in plots for the file created with SELECTED_OUTPUT 1(labels areinserted sequentially). Line 9 inserts two empty lines (a dataset separator in GNUPLOT) and alabel in the file created with SELECTED_OUTPUT 1.

SELECTED_OUTPUT 1-file results.out-pH true–totals Na K
-molalities CO2 HCO3- CO3-2–equilibrium_phases Calcite–active false
PUNCH 1 index "Dataset"
PUNCH 1 gnuplot "Mixing line"
Code 45 – SELECTED_OUTPUT and PUNCH keyword.
The curve options in SELECTED_GRAPH, -fig_y for a XY plot for example are used to specify whichdataset in a multi-dataset file has to be included in a plot. More specifically, <i(<txt:Option>|<int:Index>)> is used for this purpose. Examples are given in Code 46 where Caconcentrations are plotted as a function of pH assuming a multi-dataset file.
••••

Line 5 – when no index is given, only the first dataset is plotted.
Line 6 – the curve option i all will plot all datasets in the file as separate curves
Line 7 – only the curve corresponding to the dataset with index 2 will be plotted. This is the third
dataset in the file as the index number is zero-based (index 0 is the first dataset).
Line 8 – different datasets, specified by the zero-based index are plotted. The list of indices starts
with the letter l and consists of comma-separated indices (0,2 in the example). A range a indicesis defined with a hyphen (6-8 in the example). Spaces are not allowed.

15 Punch is only available in HPGeochemistry and HPx. It is not available in standard PHREEQC.

PHREEQC-HPx Course, Prague, September 3-4, 2020

SCK CEN/39560003

141
ISC: PublicSession 4
SELECTED_GRAPH 1-fig_nr 1-fig_name "Name"-fig_x pH-fig_y Ca-fig_y Ca <i all>-fig_y Ca <i 2>-fig_y Ca <i ll0,2,6-10>-fig_xt "Label x"-fig_yt "Label y"
Code 46 – Examples of adding a dataset index in SELECTED_GRAPH to plot a specific dataset in a multi-dataset
file

Hands-on: Cation exchange isotherms
Create the simulation case, database and components
Open the Project Manager
Select the Course project
Click Create New CasebuttonGive Name – “Example RE_ex 01”
Give Description (optional) – “Exchange isotherm”
Click OK button
Select the simulation case
Click Select button
Go to the Database tab
Click Components button
Select relevant components: Ca, Na
Click OK button

Step 1: Calculate the cation isotherm at a specific isotherm
We create an extra input file (ReactionDefinition.in) that calculates the isotherm for a specificconcentration Na. The file will be included three times, each time with a different Na concentration.
Therefore, we define a global variable for the Na concentration. We first equilibrate the solution with an
exchanger, and then remove Na and add Ca (!keep charge balance) until all Na is replaced by Ca.
The results are written in a multi-dataset output file. That file is controlled such that the initial calculations
are not included, and that the three isotherms are separated by a two empty lines.
Starting conditions for exchange isotherm
Create a new input tab and save it as ReactionDefintion.in
Define a solution with the Na concentration (counterbalance by Cl) – use a global variable(@ConcNa) as the concentration
Define an exchanger and equilibrate with the solution. Use a global variable (@CEC) to define thecapacity of the exchanger.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

142
ISC: Public

Equilibria at the solid-liquid interface
These initial composition calculation are not needed for the plotting of the isotherms. We do notsave them in the output file. Writing to the output file is suppressed by:
SELECTED_OUTPUT 1-active false
Because the exchanger is defined in equilibrium with the solution, a batch calculation between
the solution and exchanger is not required. The reaction calculation can be excluded by adding:
Include END

USE solution noneCalculate the isotherm
We first set the instructions for the multi-dataset file
Activate the selected output file 1
Now, the instructions for the isotherm are given:
Include the USE statement to specify the solution and exchanger to be used in the calculation
USE solution 1
USE exchange 1
Use REACTION to remove all Na and replace it by Ca – for charge balance, for each Na removed,
only 0.5 Ca are removed in 100 steps.
The total amount of Na is the initial amount on the exchanger and the amount in solution. An
inline expression can be used to calculate the total of Na to be removed
Include END@<@CEC+@ConcNa> in 100 steps

Step 2: Repeat the calculation for different Na concentrations
Repeat the calculation three times
Go to GeoChemistry Panel
Define the global variable CEC
Insert three times
A (re)definition of the global variable ConcNa, each time with a different value
Delineate another dataset:
PUNCH 1 index “@ConcNa M”
PUNCH 1 gnuplot “@ConcNa M”
The first line (index) is only for the first block – as it corresponds with the first dataset, it
is not needed to include two empty lines. The PUNCH statement assigns a label for plotting
to the first dataset in the file..
The second line (gnuplot) is used for the two subsequent blocks. A “gnuplot” delineation(headings and two empty lines) is written to the file defined by SELECTED_OUTPUT 1. Alabel is also assigned to the dataset using the global variable.
Include ReactionDefinition.in
Include END
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

143
ISC: PublicSession 4
Define the selected output
In the first block, include the definitions for the output
SELECTED_OUTPUT with the relevant definitions (file name)
USER_PUNCH – calculate the equivalent fractions of Ca in
Solution: use tot(“”)
Exchange: use mol(“”)
At the end, include the definition for a graph
SELECTED_GRAPH with the relevant definitions
For -fig_y
Use the variable that expresses the equivalent fraction of Ca on the exchanger
Use curve options < i all >i – indicates the index
all – indicates that all datasets in the multi-dataset will be printed.

OutputFigure 55 shows the resulting graph with the three isotherms. When the solution is less diluted (e.g. 0.1
versus 0.001 M), selectivity for Ca decreases. In a more concentrated solution, monovalent cations are
preferred on the exchanger compared to bivalent cations. This means that in sea water, the exchanger
is dominated by Na whereas in fresh water, the exchanger will be dominated by Ca.

Figure 55 – Ca/Na exchange isotherms at three NaCl concentrations (example RE_ex 01)

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

144
ISC: Public

Equilibria at the solid-liquid interface

Reaction path – Changing exchange composition during evaporation
Problem Definition
During soil evaporation, the soil pore water is concentrating. Consequently, the selectivity for cations on
the exchanger will change in case of heterovalent cation exchange.
Consider a saturated soil sample (porosity of 0.4) in which the pore water (0.1 M NaCl, 0.01 M CaCl2 – astrongly concentrated solution is chosen to see the effect on the exchanger) in equilibrium with a cationexchange complex. The capacity of the exchanger is 5.5 meq / 100 g soil.
Calculate the equivalent fraction of Ca and Na in solution and on the exchanger when the soil sample isdrying from full saturation to the residual water content of 0.035. Assume further a particle density of
2.65 g/cm³.

Conceptual model
The problem involves following components:
•
•
•
•
•

Equilibrium reactions
Two phases (aqueous and inert solid phase)
Surface-liquid interface – Cation exchange process
Well-mixed
Adding and removing of water

We model evaporation. .
Figure 56 shows a schematic representation of the key elements to calculate a cation exchange during
evaporation. Here, the solid phase is present but it is considered to be inert (e.g. no dissolution andprecipitation). Exchange processes at the solid-solution interface are included (EXCHANGE) and water isremoved in a step wise manner using the keyword REACTION

Figure 56 – Schematic representation of the model to calculate cation exchange during evaporation of water (left)
and key elements of the computer model..

Hands-on: Exchange during evaporation
Create the simulation case, database and components
Open the Project Manager
Select the Course project
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

145
ISC: PublicSession 4
Click Create New CasebuttonGive Name – “Example RE_ex 02”
Give Description (optional) – “Exchange composition during soil evaporation”
Click OK button
Select the simulation case
Click Select button
Go to the Database tab
Click Select button
Select relevant components: Ca, Na
Click OK button

Step 1: Calculate the exchanger capacity
The input capacity of the exchanger (X) in EXCHANGE is in mol in your system. If we consider a systemof 1 bulk dm3, the capacity needs to be expressed in mol/bulk dm³. Note that also the volume of waterin SOLUTION should then be expressed in L / bulk dm3.
Go to Global Definitions Panel
Select the Global Variables item
In Edit – Create three global variables: CEC (in meq/100g), Porosity and ParticleDensity(in g/cm³)
Validate by Values
Select the Transformations item
Create a new variable – CECmol_dm3
Add the equation to transform CEC in CECmol_dm3
solid = 1 - @Porosity
cec_g = @CEC / 100
solid * @ParticleDensity * cec_g
Step 2: Define the geochemical problem
Go to Geochemistry Panel
Define the solution with SOLUTION – use the global variable Porosity to define the amount ofwater
To exclude the initial calculation from the selected output file (and graph), put an END after the solution
definition. When you define SELECTED_OUTPUT after this end (as in this example), the initial calculationof the solution will not be written to the output file. You need to have USE solution in the next
simulation.
Put END
Put USE solution 1

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

146
ISC: Public

Equilibria at the solid-liquid interface
Define the exchanger with EXCHANGE in equilibrium with the solution – use the global variable
CECmol_dm3 to define the capacity of the exchanger
During the evaporation process evaporation, water (H2O) is removed from the system whereas all othercomponents remain in the system. The evaporation process is modelled using REACTION by removing
the reactant H2O from the system. The total volume of water that is removed from the system (1 dm³
bulk soil) is the difference between full saturation and the residual water content. Note that the amountin REACTION is expressed in moles (see page 58) – the volume of water needs to be converted in molesof water.
Define the Reaction – use a inline variable (@<>) to convert from volume to moles. Calculate theevaporation process in 100 steps.

Step 3: Output and plots
We create two graphs – one with the equivalent fractions on the exchanger and one with the equivalent
fractions in the pore water as a function of water content. On both graphs, we plot the ionic strength on
the secondary Y-axis.
Define SELECTED_OUTPUT and request at least the water content (-water) and the ionicstrength (-ionic_strength)
Define USER_PUNCH and calculate equivalent fractions of Ca and Na in solution and on theexchanger
The script function mol(“”) returns the concentration of a species (mol/kgw). To calculate theequivalent fraction of an exchange species, the numerator should have the same unit as thedenumerator. If you use the global variable CECmol_dm3, you either have to convert the numerator
to mol/dm³ (mol(“”)*tot(“water”)) or you convert the denumerator to mol/kgw(get_var(“”)/tot(“water”))
Define SELECTED_GRAPH and define two graphs.
The secondary Y-axis is defined by –fig_sy followed by the definition of the curve(s). The
secondary Y-axis label is defined by –fig_syt.
Enhanced text format can be used for the key of the equivalent fractions on the exchanger:-fig_y BCa <k "{/Symbol:Italic b}_{e,Ca}">-fig_y BNa <k "{/Symbol:Italic b}_{e,Na}">/ indicate fonts specifier – here Symbol is defined to define Greek letters:Italic will have an italic font (:Bold a bold font)
In each of the two graphs, we will put the key inside at the bottom left corner. Gnuplot commands
for all plots defined in a single SELECTED_GRAPH can be defined between the identifiers –startand –end. Note that the gnuplot commands between –fig_start and –fig_end are applied
for a specific plot only.
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

147
ISC: PublicSession 4-startset key bottom inside left
set style line 1 lc "red" lw 3
set style line 2 lc "blue" lw 3
set style line 3 lc "green" lw 2
-end
Save and Run

Output
The two resulting graphs are shown in Figure 57 and Figure 58. In line with the exchange isotherms
calculated with the previous example (see Figure 55), the selectivity of Ca decreases during theevaporation process (i.e., decreasing water content resulting in increasing ionic strength of the solution).
The equivalent fractions of Na on the exchanger increases whereas they decrease in the pore water.

Figure 57 – Equivalent fractions of Na and Ca in the pore water (left axis) and ionic strength (right axis) of a soil inequilibrium with an exchange complex during a drying event.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

148
ISC: Public

Equilibria at the solid-liquid interface

Figure 58 – Equivalent fractions of Na and Ca on the exchanger (left axis) and ionic strength (right axis) of a soil inequilibrium with an exchange complex during a drying event.

Reactive Transport Modelling (HP1) – Cation Exchange
Transport and Cation Exchange – Single Pulse
Learning outcomes
Implement cation exchange processes in a reactive transport problem
Create time series plots at given nodes of geochemical variables using SELECTED_GRAPH

Problem Definition
This example is adapted from Example 11 of the PHREEQC manual (Parkhurst and Appelo, 1999, see
online documentation). We will simulate the chemical composition of the effluent from an 8-cm column
containing a cation exchanger. The column initially contains a Na-K-NO3 solution in equilibrium with thecation exchanger. The column is flushed with three pore volumes of a CaCl2 solution. Ca, K and Na areat all times in equilibrium with the exchanger. The simulation is run for one day; the fluid flux density isequal to 24 cm/d (0.00027777 cm/s).
The column is discretized into 40 finite elements (i.e., 41 nodes). The example assumes that the same
solution is initially associated with each node. Also, we use the same exchanger composition for all
nodes.
The initial solution is Na-K-NO3 solution is made by using 1 x 10-3 M NaNO3 and 2 x 10-4 KNO3 M. The
inflowing CaCl2 solution has a concentration of 6 x 10-4 M. Both solutions were prepared under oxidizing
conditions (in equilibrium with the partial pressure of oxygen in the atmosphere). The amount of
exchange sites (X) is 1.1 meq / dm³ soil.
In this example, only the outflow concentrations of Cl, Ca, Na, and K are of interest.
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

149
ISC: PublicSession 4

Conceptual model
The problem involves following components:
•
•
•
•
•
•
•
•

Equilibrium reactions
One phases (aqueous)
Surface-liquid interface – Cation exchange process
Well-mixed
Homogeneous domain
1D flow
1D advective-dispersive solute transport
Constant boundary condition

Figure 59 shows a schematic representation of the conceptual model.

Figure 59 – Schematic representation of a 1D reactive transport problem involving cation exchange reactions.

Hands-on
Step 1 – Define main settings and water flow
Project Manager
Button: "New"Name: "T-CEC-01"
Description: " Transport and Cation Exchange, single pulse "Button: "OK"Main Processes
Heading: Transport and Cation Exchange, single pulse
Check: “Water Flow”
Check: "Solute Transport"
Select: "HP1 (PHREEQC)" – “HP version 2020”
Button: "Next"
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

150
ISC: Public

Equilibria at the solid-liquid interface
Geometry Information
Depth of the Soil Profile: 8 cm
Button: "next"
Time Information
Time Units: Seconds
Final Time: 86400
Initial Time Step: 180
Minimum Time Step: 180
Maximum Time Step: 360
Button: "Next"
Print Information
Number of Print times: 12
Button: “Select Print Times”
Button: "Next"
HP1 – Print and Punch Controls
Button: “Next”
Water Flow – Iteration Criteria
Lower Time Step Multiplication Factor: 1
Button: “Next”
Water Flow – Soil Hydraulic Model
Button: “Next”
Water Flow – Soil Hydraulic Parameters
Catalog of Soil Properties: Loam
Qs: 1
Ks: 0.0002777
Button: “Next”
Water Flow – Boundary Conditions
Upper Boundary Condition: Constant Pressure Head
Lower Boundary Condition: Constant Pressure Head
Button: “Next”
Step 2 – Define Solute transport and geochemical model
Solute Transport – General Information
Button: “Next”
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

151
ISC: PublicSession 4
Solute Transport – HP1 Input Data
Go to Database tab
Select components: Na, K, Ca, Cl, N(5)Go to Solution Definitions tab
Define SOLUTION 1001 - K-Na-N(5) solution: Use pH to charge balance the solution, and adapt
the concentration of O(0) to be in equilibrium with the partial pressure of oxygen from theatmosphere.
Define SOLUTION 3001 - use pH to charge balance the solution, and adapt the concentration ofO(0) to be in equilibrium with the partial pressure of oxygen from the atmosphere.
Go to GeoChemistry tab
Define EXCHANGE 1001 – define it for material 1, and equilibrate with the initial solution.
Solute Transport – Transport Parameters
Bulk Density: 1.5
Disp: 0.2 cm
Button: “Next”

Solute Transport – Boundary Conditions
Upper Boundary Condition: Concentration Flux
Add the solution composition number for the upper boundary condition
Lower Boundary Condition: Zero Gradient
Button: “Next”
Step 3: Geometry and initial conditions
Soil Profile – Graphical Editor
Menu: Conditions -> Profile Discretization
Number (for sidebar): 41
Menu: Conditions -> Initial Conditions -> Pressure Head Button: "Edit Condition"
Select All Top Value : 0
Menu: Conditions -> Observation Points Button: "Insert"
Insert a node in the middle and at the bottom
Menu: File -> Save Data
Menu: File –> Exit
Soil Profile – Summary
Button: “Next”
Run ProjectNote: This exercise will produce following warnings: "Master species N(3) is present in solution n but is
not transported.". The same warning occurs for N(0). N(3) and N(0) are two secondary master species
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

152
ISC: Public

Equilibria at the solid-liquid interface
from the primary master species N. Only the secondary master species N(5) was defined as a component
to be transported (Solute Transport – HP1 Components). HP1, however, checks if all components which
are present during the geochemical calculations, are defined in the transport model. If not, a warning
message is generated. In our example, the concentrations of the components N(0) and N(3) are very low
under the prevailing oxidizing conditions. Therefore, they can be neglected from transport problem. If
one wants to avoid these warnings, one have to either include N(0) and N(3) as components to betransported or define an alternative primary master species representing nitrate (such as Nit-) using
SOLUTION_MASTER_SPECIES and SOLUTION_SPECIES.
OutputFigure 60 shows the breakthrough curves of the different ions (components) at the bottom of the profile.
The figure is obtained via:
Post-Processing – Observation Points – All Concentrations – Observation node 2.
Right click – Chart Designer
The figure is saved as “BTC – Bottom.vct” in the project folder.
The concentrations of node 41 (the last node) are plotted against time. Chloride is a conservative soluteand arrives in the effluent at about one pore volume. The sodium initially present in the column
exchanges with the incoming calcium and is eluted as long as the exchanger contains sodium. The
midpoint of the breakthrough curve for sodium occurs at about 1.5 pore volumes. Because potassium
exchanges more strongly than sodium (larger log K in the exchange reaction), potassium is released
after sodium. Finally, when all of the potassium has been released, the concentration of calcium increases
to a steady-state value equal to the concentration of the applied solution.

Figure 60 – Breakthrough curves of Na, K, Cl and Ca at the bottom of the profile during steady-state transportwith cation exchange (example T-CEC-01)

Output – Creating time series plot
Using the SELECTED_GRAPH keyword in the input files, plots similar to Figure 60 can be made forcomponents and other geochemical variables. Instead of the default figure type “XY”, the identifier –fig_type time is used. The time variable is then automatically assigned to the X-axis. Variables to be
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

153
ISC: PublicSession 4plotted to the Y-axis are defined with –fig_curve (and –fig_scurve) similarly to the –fig_ydefinitions. The label is defined as –fig_curvet.
To indicate which nodes you want to plot, the identifier –fig_multi is used – see Code 47:••

Line 2 will make a separate plot for each of the nodes (default)
Line 3 will make a separate plot for the observation node 1, 4, 5 and 6.

-fig_multi <:txt:multi>-fig_multi all-fig_multi l1,4-6
Code 47 – Syntax and examples of –fig_multi as used in SELECTED_GRAPH to make time series or profile
1
2
3

plots during reactive transport simulations.

In the next instructions, these graphs are made, for the two nodes, for the concentrations of Na, K, Ca
and Cl and for the equivalent sorbed fractions of Na, K and Ca.
Open Solute Transport – HP1 Input Data
Go to Output Panel
Add SELECTED_OUTPUT and for the total concentrations of Na, K, Cl and Ca
Add USER_PUNCH and calculate the equivalent fractions
Add SELECTED_GRAPH
Add template for time series
HP Keywords -> Output -> Selected_Graph -> Templates -> Time series
Change:-fig_nr to 1
-fig_name to “Concentrations”-fig_multi to all-fig_curve to &tot Cl Na K Ca
-fig_curvet to “Concentration (mol/kgw)”
Add gnuplot commands: set yrange [0:0.0014]
Add a second time series template and adjust the identifiers for plotting the equivalent fractions(variables defined in USER_PUNCH)
Add Gnuplot commands for all graphs-startset key inside bottom
set xrange [0:86400]
set xtics 14400-end
The resulting plots are shown in Figure 61.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

154
ISC: Public

Equilibria at the solid-liquid interface

Figure 61 – Breakthrough curves (left) and equivalent sorbed fraction (right) at 4 cm (top) and 8 cm (bottom) for
example T-CEC-01.

Transport and Cation Exchange – Multiple Pulse
Problem Definition
We use the same set-up is in the previous example, but time variable concentrations are applied at thesoil surface.
Following sequence of pulses are applied at the top boundary:0 – 8 hr: 6 x 10-4 M CaCl2
8 – 18 hr: 5 x 10-6 M CaCl2, 1 x 10-3 M NaNO3, and 2 x 10-4 M KNO3
18 – 38 hr: 6 x 10-4 M CaCl2
38 – 60 hr: 5 x 10-6 M CaCl2, 1 x 10-3 M NaNO3, and 8 x 10-4 M KNO3

Conceptual model
The model is similar to the previous example, except a time varying upper solute boundary condition(Figure 62).

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

155
ISC: PublicSession 4

Figure 62 – Schematic representation of a 1D reactive transport problem involving cation exchange reactions witha time varying upper solute boundary condition.

Hands-on
Step 1 – Define main settings and water flow
Project Manager
Select and copy project “T-CEC-01”
Name: "T-CEC-02"
Description: " Transport and Cation Exchange, multiple pulses "Button: "OK"Main Processes
Heading: Transport and Cation Exchange, multiple pulses
Button: "OK"Time Information
Time Units: Change to hours
Final Time: 60
Check: Time-Variable Boundary Conditions
Number of records: 4
Button: "Next"
Button: “OK”
Step 2 – Define Solute transport and geochemical model
Solute Transport – HP1 Input Data
Go to Geochemistry Tab
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

156
ISC: Public

Equilibria at the solid-liquid interface
Add SOLUTION 3002 (8-18 hr) and 3003 (38-60 hours)
Add SOLUTION 4001 (pure water)
Go to Output Tab
Change the xrange and xtics to [0:60] and 15 (SELECTED_GRAPH)
Button: “OK”
Variable Boundary Conditions
Fill in the time, and the solution composition number for the top boundary (3001/3002/3003) andthe bottom boundary condition (4001).
Step 3: Geometry and initial conditions
Soil Profile – Graphical Editor
Menu: Conditions -> Profile Discretization
Menu: Conditions -> Observation Points Button: "Insert"
Insert nodes at 2, 4, 6, and 8 cm
Menu: File -> Save Data
Menu: File –> Exit
Soil Profile – Summary
Button: “Next”
Run ProjectOutputExplore the output. As examples, Figure 63 gives the K concentration at different depths in the profile.
Figure 64 shows the outflow concentration. The first pulse is identical to the single pulse project. Thenadditional solute pulses of different solution compositions will restart the cation exchange processdepending on the incoming solution composition.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

157
ISC: PublicSession 4

Figure 63 – Time series of K at different observation nodes during transport with variable boundary conditions(Example T-CEC-02).

Figure 64 – Breakthrough curves of Cl, Na, K, and Ca during transport with variable boundary conditions (Example
T-CEC-02).

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

158
ISC: PublicSession 5 – Kinetics

Session 5
Many processes in porous media are not in instantaneous equilibrium; it takes some time. Such kinetic
processes occur for all geochemical processes - degradation, complexation, sorption,
dissolution­precipitation, oxidation-reduction. Radionuclide decay is another example of a kinetic
process. PHREEQC has an large flexibility in defining rate equation and allows thus for the definition ofrelatively complex reaction networks including e.g. colloid retention and colloid-affected mobility,(bio)degradation pathways etc.

Learning Outcomes
Implement different types of rate equation using the RATES data block
Use KINETICS for simulating kinetic processes
Set-up kinetic reaction paths using HPGeoChemistry
Include kinetic reactions in a reactive transport simulation using HPx

Numerical world - Implementing kinetics
Figure 65 shows the PHREEQC keywords for implementing kinetic processes. The keywords RATES, and
KINETICS are introduced in this section.
PHREEQC/HPGeochemistry/HPx do not have predefined and hard-coded rate equations as is the casein some other geochemical or reactive transport codes. The user needs to write a script that calculate
the rate using the keyword RATES. The outcome of the script is the mole transfer from the kinetic
reactant: a positive value indicates a decrease in the kinetic reactant; a negative value an increase.
To simulate kinetic reactants, the keyword KINETICS is used.

Figure 65 – Implementation of the model for kinetic processes.

Implementing the rate equation: RATESThe scripts defined in RATES define the rate equations. Code 48 shows the typical syntax of the RATES
data block:•Line 1: start of a RATES data block. Within a block, one or more rate equations can be definedby repeating lines 2 – 8.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

160
ISC: PublicKinetics•••

Line 2: the rate name. The rate name is used in KINETICS to include a rate in the geochemical
model. If the rate name corresponds to a phase name defined in PHASES, the stoichiometry of
the phase name is used to define which components are added or removed from the aqueous
phase. If the rate name does not correspond to a phase name defined in PHASES, thestoichiometry of the components must be defined in the kinetic block (using -formula).
Examples are given later in this section.
Line 3: indicates the scripting language. If omitted, the standard PHREEQC BASIC language isused. When –mybasic is present, the alternative scripting language is used.
Line 4 … Line 8: Between the –start and –end, the rate equation is defined. .

Typically, a number of statements leads to a rate (line 6) that is expressed in mol/T/L3 (instantaneousrate). T is by default seconds in PHREEQC 16. L3 is the volume of the considered geochemical system 17.
The instantaneous rate is integrated over a time interval (the time step) using the commandtime() 18 (line 6) to obtain the mole transfer of the kinetic reactant to the solution (mol/L3). As
mentioned above, a positive value will decrease the amount of the kinetic reactant whereas anegative value will increase the amount of the kinetic reactant. Depending on the stoichiometry of
the components, concentration of the components will decrease or increase. If the rate namecorresponds to a phase name, the mineral dissolves with a positive mole transfer will dissolve themineral and precipitates with a negative mole transfer.
Finally, the value of the mole transfer is returned with the save() 19 statement (line 8).

RATES<:txt:name>-mybasic-start<:txt:script>
mole = rate*time()save(mole)-end
Code 48 – Typical syntax of RATES.

More information on the keyword RATES is available on the online manual page.

Kinetic reactants: KINETICS RATESThe KINETICS block provides different type of information (Code 49):•••

To define which rate equations are included in the geochemical model (lines 2-5)
The reaction path (line 6)
The numerical aspects of solving the rate equations (not in Code 49).

Code 49 gives a basic syntax of KINETICS limiting the information to the aspects used in the examples.
To incorporate a given rate equation in the geochemical model, the rate name is defined in the
KINETICS block (line 2). The rate name is followed by a few identifiers (here lines 3-5) with specificdetails:•

Line 3: The identifier –m0 defines the initial amount of the kinetic reactant in moles/L3

16 In HPx, T is linked to the time unit defined in the HYDRUS graphical user interface.
17 In HPx, L3 is always 1 dm3.
18
19time in standard PHREEQC BASIC.
save in standard PHREEQC BASIC.

PHREEQC-HPx Course, Prague, September 3-4, 2020

SCK CEN/39560003

161
ISC: PublicSession 5••

Line 4: When a rate name is not linked to a phase name, the identifier -formula is required toindicate the chemical formula or elements and their stoichiometric coefficients that are involvedin the mass transfer . The amount of a given component is obtained by multiplying the moletransfer from the rate equation (i.e. the value returned by the save() statement) with thestoichiometric coefficient. The sign of the amount defines if the element is removed from theaqueous phase (negative sign) or added to it (positive sign).
Line 5: The identifier -parms defines a list of parameters that can be used in the rate equation.
by the function parm(<:int:index>) ()where index indicates the position in the list (startingfrom a value of 1).

Line 6 defines the reaction path (-step) by indicating the total time of the reaction path (in seconds)
and the number of equal time intervals.

KINETICS <:int:number><:txt:name>-m0 <:flt:amount>-formula {<:txt:formula> <:flt:stoichiometry>}-parms {<:flt:parm>}-steps <:total time:flt> in <:number of steps:int>
steps
Code 49 – Basic syntax of the KINETICS block
More information on the keyword KINETICS is available on the online manual page.

Examples of RATES and KINETICS
Reaction Path – Kinetic dissolution of a mineralIn this example, we treat the kinetic dissolution and precipitation of minerals. Rate equations are often
based on the so-called transition state theory (Aagaard and Helgeson 1982). An affinity term is included
to capture the decreasing dissolution rate when the solution is closer to equilibrium conditions (i.e. when
the saturation state is closer to 1). The rate equation is then complicated by including the effects solution
composition, presence of inhibitors or catalysts, temperature, change in reactive surface area andsurface- or transport-controlled (e.g. shrinking core model) processes (see e.g. Mayer et al. (2002)).
In a first exercise, the kinetic dissolution of quartz is simulated with a simple rate equation including theeffect of reactive surface area, pH of the solution and the affinity term.

Learning outcomes
Understand the different factors of a rate equation for mineral dissolution and precipitation basedon the transition state theory.
Implement a kinetic mineral dissolution rate equation based on the transition state theory.
Illustrate how to run different cases using the same input file.

Mathematical world: Mineral dissolution kinetics based on Transition State
Theory –based rate equation
A general rate equation for a surface-controlled mineral dissolution/precipitation reactions can be
written as:

R = S Fk ( C , T ) FT
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

(24)
162
ISC: PublicKinetics
where
R [mol/L³/T] is the dissolution/precipitation rate
SF [L²/L³] is a factor accounting for the reactive surface area in the sysemk [mol/L²/T] is the kinetic rate constant that may depend on the geochemical conditions (C) andtemperature (T)FT [-] is a chemical affinity factor which will be based on the transition state theory

Chemical affinity
The chemical affinity term expresses how the rate depends on the Gibbs free energy of the saturation
state. The chemical affinity A is related to the Gibbs free energy as:∆Gr =∆Gr° + RT ln Q
= RT ln K eq + RT ln Q

(25)where Q is the ion activity product, and Keq is the mass action equilibrium constant. Thus:

∆Gr =
RT ln

Q
=
−A
K eq

(26)

Figure 66 illustrates the dependency of the dissolution rate on the chemical affinity.

Figure 66 – (Left) Dissolution rate of kaolinite as function of DG (Nagy et al. 1991); (Right) Dissolution rate of KFeldspar as a function of DG (Gautier et al. 1994).

Affinity term based on transition state theory
Based on the Transition State Theory (Aagaard and Helgeson 1982), the affinity term in Eq. (24) can be
expressed in a general way as (Brantley 2003, Steefel et al. 2015): n∆G  =
FT  exp  − 1
 RT  = ( Ω n − 1)

m

(27)

m

where n and m are experimentally or theoretical derived exponential dependence and Ω is the saturation
ratio.
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

163
ISC: PublicSession 5
The linear form of Eq. (27) is (Brantley 2003):

FT = ( Ω n − 1)(28)

Figure 67 shows the relation where n = 1 or 2 for Ω smaller than 1, i.e., solution is undersaturated andmineral will dissolve. For far from equilibrium conditions (strongly undersaturated), the rate isindependent of the chemical affinity. Near equilibrium, the rate is linearly related to the chemical affinity
and the slope is proportional to n.
When the solution is oversaturated (Ω larger than 1), the rate is always linearly depending on thechemical affinity.

Figure 67 – FT as a function of saturation sate for 2 values of n (m=1 in Eq. (27)) for dissolution (Ω < 1).

Temperature dependency
Figure 68 shows that the dissolution rate of quartz on a semilog plot is linear with respect to thereciprocal of temperature (in Kelvin). This is typically described with Arrhenius equation:

k = Ae − Ea / RT

(29)

where
A is a pre-exponential factor (interpreted as a constant independent of temperature)
Ea is the activation energy (kJ/mol)
R is the universal gas constant (8.313 10-3 kJ/mol/Kelvin)
T is temperature (Kelvin)
If the rate constant at a given temperature is known (say k0 at T0), the rate constant at the othertemperature is given by:k = k0 e

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

 Ea  1 1   −  −  
 R  T T0  

(30)

164
ISC: PublicKinetics

Figure 68 – Dissolution rate of quartz as a function of temperature (Dove 1995)

Reactive surface area
As the rate constant is typically expressed in mol/m², the total reactive surface area of a mineral is thesystem needs to be known. During mineral dissolution and precipitation, the reactive surface areachange with reaction progress. Therefore, the factor SF in Eq. (24) accounts for properties of the mineralcrystals such as changes in crystal size, size distribution, ... and their evolution. For a monodisperse
population of uniformly dissolving/growing spheres, SF can be expressed as:2

m 3
S F = A0   m0 

(31)

where
A0 is the initial reactive surface area (m²)
m0 is the initial amount of the mineral in the system (moles)
m is the actual amount of the mineral in the system (moles)
Dependence of chemical composition
Steady-state dissolution rates of minerals are often dependent on the solution geochemistry (which is
linked to the speciation of the mineral surfaces). This dependency can be described with a relation
between solution pH and rate; examples are shown in Figure 69 for quartz and albite.. Table 5 lists somerate equations relating the steady-state dissolution rate to the solution pH.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

165
ISC: PublicSession 5
Figure 69 – Examples of pH-dependency of steady-state dissolution rates of quartz (left) and albite (right) (Walther
1996).
Table 5 – Equations expressing the pH dependency of the steady-state dissolution rate for some minerals.
MineralRate Equation(mole m-2 sec-1)

pH-range

Reference

AlbiteLog(k0) = -14.10 + 0.3 pH

7.5-12.5

Walther (1996)

K-feldspar

Log(k0) = -15.10 + 0.3 pH

Kaolinite

Log(k0) = -12.4 – 0.12 pH
Log(k0) = -15.68 + 0.3 pH

1-8
8-13

Walther (1996)

Na-Montmorillonite

k0 = 10-11.39 {H+}0.34
k0 = 10-12.31 {OH-}0.34

pH < 8
pH > 8.5

Huertas et al., 2001
Huertas et al., 2001 Huertas et al. (2001)

QuartzLog(k0) = -15.44 + 0.3 pH

7-12

Walther (1996)

Walther (1996)

Problem Definition – Kinetic quartz dissolution
In a batch type of experiment, 1 mole of quartz is dissolving in 1 litre of water for a period 100 days.
During the experiment, the solution is buffered to a pH of 8.5 by adding NaOH. A model predicts the Siconcentration and the saturation index of quartz in the. solution.
Quartz dissolution is described with following equation:=
R A0k (1 − Ω )(32)where k is given as (only accounting for pH dependency):(33)

k=−15.44 + 0.3pH

Conceptual model
The problem involves following components:
•
•
•
•
•

Equilibrium reactions
Two phases (aqueous, solids)
Kinetic dissolution of solid phaseWell-mixed
Buffer pH

We model a heterogeneous kinetic reaction path. .
A schematic representation of the model for kinetic dissolution of a mineral is show in Figure 70.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

166
ISC: PublicKinetics

Figure 70 – Schematic representation of the model to calculate kinetic dissolution of a mineral.

Hands-on
Step 1 – Create Case

Create a new case: Example K_min 01
Description: “Dissolution of Quartz”
Go to Database Panel
Button: “Components”
Select Si
Step 2 – Add Global Variables

Go to Global Definitions Panel
Item: “Global Variables” – Edit
Add global variables
We define a number of global variables to make the input more generic. We define two tabs:
One for “Experimental Conditions”
Global Variable: ReactionTime (10 days)
Global Variable: Sol_pH (8.5)
One for “Minerals” – Add group “Quartz”
Global Variable: mQuartz (10 moles)
Global Variable: AgQuartz (4 m2/g)
Global Variable: k0SlopeQuartz (0.3)
Global Variable: k0IntQuartz (-15.44)
Validate the table by using button Values
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

167
ISC: PublicSession 5
Step 3 – Add Rate Equation

Go to Additions to database Panel
Add the Fix_H phases
Help Tree -> HP Keywords -> Database -> Phases -> Templates
Add the rate equation

RATES
Quartz-mybasic-startA = parm(1) * gfw(phase_formula("Quartz")) * m()
k0log = parm(2)*(-la("H+")) + parm(3)k0 = 10^k0log
rate = A * k0 * (1-sr("quartz"))mole = rate*time()save(mole)-end
Code 50 – Script for modelling quartz dissolution using Eq. (32).

•Line 5,6: the script function parm(<:int:index>) indicates that parameter values defined in
the KINETIC block (­parms)-are used. The index is position of the parameter in the list (first
position has index 1). For example, for a list of parameters defined in kinetics as:-parms 2.5 1E-10 0.8
parm(1) will use the value 2.5, parm(2) the value 1E-10 and parm(3) the value of -5.3.

•••

Line 5 calculates the reactive surface area in m² in the system. The input variable (m²/g ) is
multiplied by the molar weight (g/mol) and the amount of quartz (mol). The function gfw(<:str:chemical formula>) gives the molar weight. phase_formula(<:str:phasename>) returns the chemical formula of the phase. m() gives the (current) amount of the
kinetic reactant (moles).
Line 8: the function sr(<:str:phase name>) returns the saturation index of the mineral.
Line 10: the sign of mole is determined by the affinity term in Line 8. In case the mineral isundersaturated, sr() is smaller than 1 and mole is a positive value indicating a mole transfer
to the solution and a decrease in the kinetic reactant quartz. For an oversaturated mineral, sr()is larger than 1; the chemical affinity term is negative and there is thus an increase in the kinetic
reactant (precipitation).

Step 4 – Add Geochemical model

Go to Geochemistry Panel
Add SOLUTION 1
Add EQUILIBRIUM_PHASES 1 with Fix_H to fix pH at @Sol_pH with NaOH as a buffer
Add KINETICS 1KINETICS 1
quartz-m0 @mQuartz-parms @AgQuartz @k0SlopeQuartz @k0IntQuartz-steps @<86400* @ReactionTime> in 100 steps
Code 51 – KINETIC block for Quartz dissolution

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

168
ISC: PublicKinetics••••

Line 2 refers to the corresponding Rate definition (Code 50). Because it is a phase from thedatabase, the reaction stoichiometry is taken from the database.
Line 3: The identifier –m0 defines the initial amount of the kinetic reactant (quartz) in moles
Line 4: The identifier -parms gives the list of parameter values that are used in the rate
definition (parm() in Code 50).
Line 5: The reaction path is defined with the identifier -steps. In this example, the final
calculation time is given (in seconds, calculated from the global variable in days multiplied by
86400 s/d). Output is given in 100 equal time intervals.
Add INCREMENTAL_REACTIONS true
This keyword indicates the time interval over which the reactions will be integrated. For thisparticular example, true means that for every step, a time step of 8640 s (0.1 d) will be used: 00.1, 0.1-0.2, 0.2-0.3, ... ; 9.9-10. If it is false (default option if INCREMENTAL_OPTIONS is notdefined), the time intervals are 0-0.1, 0-0.2, 0-0.3, ... , 0-10.
Add PRINT; -reset false to cancel printing in phreeqc.out.

Step 5 – Output and Plotting

Add SELECTED_OUTPUT
Ask for total Si concentration and the saturation index of quartz (-saturation_indices)
Add USER_PUNCH
Convert time to days (the current time is obtained with total_time())
Convert Si to mg/kw (using functions tot()()) and gfw())
Add SELECTED_GRAPH
Add three graphs:
Plot Si (mg/kgw) concentration as a function of time
Plot saturation index of quartz as a function of time (&si)
Plot the rate of quartz dissolution as a function of time (&dkin)
Insert Gnuplot commands for all graphs for a title of the plots:-start
set title "Dissolution of quartz at a pH of @Sol_pH"-end

OutputFigure 71 shows the concentration of Si, saturation index of quartz and dissolution rate of quartz as afunction of time at a pH of 8.5 (left figures).

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

169
ISC: PublicSession 5

Figure 71 – Quartz dissolution at pH 8.5 (left) and pH 10 (right). Si concentration (top), saturation index of quartz(middle) and dissolution rate of quartz (bottom) as a function of time.

Define different simulation cases
The same input file can be used to perform a simulation at a pH of 10. Instead of changing the values ofglobal variables within the input file as was illustrated previously, we use the Global Variables form
to define the cases (Figure 72).
Go to Global Definitions panel
Form “Global Variables” – Values
We define only one case for conditions with a pH of 8.5. The tool strip at the top of form (Figure 72)
allows for defining cases with different values of the global variables.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

170
ISC: PublicKinetics

Figure 72 – Defining different cases using the Global Variable form.

Add a case by using the first button
Define the Name and Description (optional)
Add another case by using again the first button
Define the Name and Description (optional)
Change the value of @Sol_pH to 10
Run the case
The plots are now for the case with the pH is buffered to 10 (right figures in Figure 71).

Reaction Path – Biodegradation using Monod kinetics
Learning outcomes
Understand the Monod rate laws and its parameters.
Implement the Monod rate equation.
Apply the Monod rate equation.

Mathematical world: Monod kinetics
Bekins et al. (1998) studied the degradation of phenol in a laboratory experiment. At higher
concentrations, there is a linearly decreasing concentration with time indicating a zero order degradation
rate. At small concentrations, the rate declines with time. This dependency of the rate on concentration
can be described by the Monod kinetic rate equation.

Figure 73 – Observed and modelled phenol concentrations during a methanogenic (bio)degradation experiment(Bekins et al. 1998)

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

171
ISC: PublicSession 5
The Monod kinetic rate equation is written asdSS= −kmaxdt
k0.5 + S(34)

where S is the concentration of the substrate (phenol) (M/L), kmax is the maximal rate (M/L/T); k0.5 is thehalf saturation constant (M/L).
The zero-order rate at very high concentrations (S >> k0.5) is equal to kmax as Eq. (34):kdS= − max S
dt
k0.5

(35)

On the other hand, at very low concentrations (S << k0.5), S can be neglected relative to k0.5 in the
denominator of Eq. (34) which then simplifies to a first-order rate equation:dS= −kmax
dt(36)

Figure 74 – The biodegradation rate as a function of substrate concentration defined by the Monod kinetic rateequation (from Appelo and Postma (2005). Parameters for the equation are: kmax = -1.39 mg phenol/L/day and k0.5
= 1.7 mg phenol / L.

Problem Definition – Biodegradation of phenol
Simulate the degradation of phenol using the rate parameters of the Monod equation (Eq. (34)) obtainedby Bekins et al. (1998): kmax = -1.39 mg phenol/L/day and k0.5 = 1.7 mg phenol / L. Consider an initialconcentration of 38.4 mg phenol/L.
Plot the degradation of phenol during 40 days.
This example is taken from example 10.7 in Appelo and Postma (2005).

Conceptual model
The problem involves following components:••••

Equilibrium reactions
One phases (aqueous)
Kinetic aqueous degradation
Well-mixed

We model a homogeneous (aqueous) kinetic reaction path. .
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

172
ISC: PublicKinetics
A schematic representation of the model for phenol degradation following Monod kinetics is show in
Figure 75.

Figure 75 – Schematic representation of the model to calculate biodegradation in the aqueous phase.

Hands-on: Biodegradation of phenol
Step 1 – Create Case

Create a new case: Example K_mod 01
Description: “Biodegradation of Phenol - Monod kinetics”
Step 2 – Define new component

A new aqueous component, Phenol, needs to be defined. This can be done with the
SOLUTION_MASTER_SPECIES and SOLUTION_SPECIES data blocks in the text input (e.g. in the Tab“Additions to the database”.
Alternatively, the dialogue window via the button “Add Information” in the Database tab can beused. The new component and species will be included automatically in the phreeqc.in input file, butthey will also appear in the help tree, autocompletion and syntax colouring when the component isselected in the component list.
Go to Database Panel
Button “Add Information”
The form for defining new aqueous components (solution master species) as part of the
Thermodynamic Data dialogue window is shown in Figure 76. The textbox Element name contains
a valid name of the element name (see SOLUTION_MASTER_SPECIES on page 34). For thisparticular example, Phenol has to be defined with arbitrary values for Alkalinity and the gramformula weight.
When the box “Include default identity reaction” is checked, the default definition of the primary
master species is included as well in the phreeqc.in input file. Otherwise, the user needs to definethat species with the solution_species data block (see page 35) or using the form in the tab
Solution Species in Figure 76.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

173
ISC: PublicSession 5

Figure 76 – Dialogue window to add Thermodynamic data – Aqueous Components

Element Name: “Phenol”
Alkalinity: 1
Gram formula weight for pure element: 1
Button: “Add”
Button: “OK”
Button: “Components”
Select Phenol
Step 2 – Add Global Variables

Go to Global Definitions Panel
Item: “Global Variables” – Edit
Add global variables (in one or more tabs or groups):
ExperimentalTime 10

day

Duration of the experiment

PhenolKMax

1.61E-8

g/L/s

Phenol – Monod: Maximum rate

PhenolKHalf

1.7E-3g/L

Phenol – Monod: Half-saturation constant

Validate the table by using button Values
Step 3 – Add Rate Equation

Go to Additions to database Panel
Add rate equation

RATES
Phenol-mybasic-start
S = tot("Phenol")if S < 1e-9 thenrate = 0
else
rate = - parm(1) * S /(parm(2) + S)
endifmole = rate * time()save(mole)-end
Code 52 – Script for modelling Phenol biodegradation using the Monod equation
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

174
ISC: PublicKinetics
The rate for this equation is negative (Line 9). This indicates that the kinetic reactant will increase.
However, because an aqueous degradation reactions is described here, the kinetic reactant is not aphysical entity (in contract to the example of quartz dissolution in which the kinetic reactant
represents the mineral quartz). This negative value has two consequences:••

The amount of the kinetic reactant will never be 0, thus the reaction will never stop because ofa zero amount of kinetic reactant (of course, when the Phenol concentration is zero, the rate ofdegradation is zero as well);
The negative value corresponds to removal from the aqueous phase. Note that thestoichiometry defined in –formula (Code 53) is important as well.

Step 4 – Add geochemical model

Go to Geochemistry Panel
Add SOLUTION with the concentration of Phenol
Add KINETICSKINETICS 1
Phenol-formula Phenol 1-m0 0-parms @PhenolKMax @PhenolKHalf-steps @<@ExperimentalTime * 86400> in 100 steps
Code 53 - KINETICS data block for Phenol degradation using Monod kinetics
Line 3 defines the stoichiometry of the kinetic reactions (-formula). The negative value from therate equation (Code 52) multiplied with the stoichiometry for Phenol (1, Code 53) gives a negativevalue for Phenol that will be removed from the solution.
Add INCREMENTAL_REACTIONS true
Add PRINT; -reset false
Step 5 - Output and Plotting

Add SELECTED_OUTPUT
Add USER_PUNCH
Convert time to days (the current time is obtained with total_time()()))
Convert Phenol to mg/kw (using functions tot())
Add SELECTED_GRAPH
Add Phenol concentrations as a function of time

OutputFigure 77 shows the evolution of Phenol as a function of time.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

175
ISC: PublicSession 5

Figure 77 – Phenol concentration as a function of timeFirst-order kinetic degradation network: Degradation of
PCE
In this exercise, we will introduce a first-order kinetic reaction equation and then extend the example toa degradation network with as sequential and parallel degradation path. As an example, the degradation
network of perchloroethylene (PCE) is used, but the rate equations can be used for other first-orderdegradation networks as well including radionuclide decay with branching. Note that PHREEQC and HPx
are flexible to include a wide range of reaction networks with sequential, parallel, reversible, andconvergent pathways .

Learning outcomes
Implement a first-order rate kinetic degradation network with sequential and parallel reactions(diverging and converging)Use input in external files to be used in different projects or cases (both in HPGeochemistry andHPx)
Implement boundary problems and initial value problems (contaminated soil) for single and dual
porosity (MIM) solute transport problems
Ask for output variables summed/integrated over a given soil domain

Background
Perchloroethylene (PCE, also called tetrachloroethylene) degrades slowly under reducing conditions,mainly by microbiological transformations. One of the most important pathways for anaerobic
biodegradation of PCE is by reductive dechlorination in a sequential way. Figure 78 (after Schaerlaekens
et al. (1999)) shows this pathway including six components: PCE, trichloroethylene (TCE), cis-1,2dichloroethylene (cis-DCE), trans-1,2-dichloroethylene (trans-DCE), 1,1-dichloroethylene (1,1-DCE), vinyl
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

176
ISC: PublicKinetics
chloride (VC). VC then eventually degrades to ethylene (ETH), which is environmentally acceptable anddoes not cause direct health effects
Although the reaction kinetics of the biodegradation depend on a variety of environmental conditions(such as redox potential, biomass, and compounds affecting solubility), the kinetics here are describedindependent of environmental conditions.

Figure 78 - Perchloroethylene (PCE) degradation pathway (after Schaerlaekens et al. (1999))

When the DCE (dichloroethylene) species are lumped, the transport problem can be formulated astransport of a single sequential degradation chain, which can then be solved using HYDRUS-1D (e.g.,
Schaerlaekens et al. (1999); Casey and Simunek (2001); see also one of the HYDRUS-1D tutorials). The
HP1 code can handle the full reaction network as defined in Figure 79 (i.e., with three separate DCE
species), considering different distribution factors α and yield coefficients γ.
α1PCE

k 1 , y1TCE

k 2 , y2α2α3cis-DCEtrans-DCE
1,1-DCE

k3
k4

y3

VC

k 6 , y4

ETH

k7

k5

Figure 79 - Degradation pathway of PCE using first-order rate constants.

In this example, the reaction network from Figure 79 is implemented in HPGeochemistry and HP1.
Parameter values are given in Table 6.
Table 6 - Definition of parameters and their values for the PCE biodegradation problem (from Case 1 and 2 in Sun
et al. (2004)). Rate parameters are for a reference temperature of 20°C.

Parameter

SymbolValues

Unit

Verification
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

177
ISC: PublicSession 5
Velocity

v

0.4

m d-1

Dispersion coefficientD

0.4

m² d-1First-order degradation rate 1

k1

0.075

d-1

First-order degradation rate 2

k2

0.070

d-1

First-order degradation rate 3

k3

0.020

d-1

First-order degradation rate 4

k4

0.035

d-1

First-order degradation rate 5

k5

0.055

d-1

First-order degradation rate 6

k6

0.030

d-1

First-order degradation rate 7k70.000001

d-1

Distribution factor, TCE to cis-DCE

α1

0.72

-

Distribution factor, TCE to trans-DCE

α2

0.15

-

Distribution factor, TCE to 1,1-DCE

α3

0.13

-

Yield coefficient, PCE to TCE

y1

0.79

-

Yield coefficient, TCE to DCE

y2

0.74

-

Yield coefficient, DCE to VC

y3

0.64

-

Yield coefficient, VC to ETH

y4

0.45

-

HP1 was verified against an analytical solution presented by Sun et al. (2004) that considers the threedaughter species of DCE. Their solution was derived for a homogeneous medium and no adsorption of
PCE or any of its daughter products. Table 1 lists parameters used in the benchmarking (Sun et al. 2004).
The upper boundary condition was a constant concentration of 1 mol/l for PCE and 0 mol/l for the othercomponents.
Numerical simulations were carried out for a spatial discretization Δx of 1 m (Pe = 1) in a 30-m longprofile, while the time steps were controlled by the performance index, ωs. Figure 80a compares theanalytical solution for PCE transport with numerical simulations using ωs = 1, 0.6 and 0.2. Results for
PCE’s daughter products with ωs = 0.1 are shown in Figure 80b-d. Notice relatively small differences in
Figure 80a between the analytical and numerical results for ωs values of 1 and 0.6. Accurate results wereobtained when smaller value of ωs were used for both TCE and its 6 daughter products (Figure 80b-d forωs = 0.1). This example shows that, as compared to the approximate rule of having ωs = 2 for accuratetransport calculations with HYDRUS-1D, ωs values smaller than 0.5 are required for the operator-splitting
approach in HP1.

11Concentration (mol / l)

0.8

(a)

HP1 (ωs = 1)
HP1 (ωs = 0.6)
HP1 (ωs = 0.2)

0.6

(b)
Concentration (mol / l)Analytical solution(Sun et al., 2004)

PCE

0.4

0.2

0

PCE

0.8

Analytical solution(Sun et al., 2004)
HP1

0.6

0.4

TCE
0.2

0
0

10

20

Distance (m)D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

30

0

10

20

30

Distance (m)178
ISC: PublicKinetics
0.25

0.1

Analytical solution(Sun et al., 2004)
HP1

0.2

cis-DCE

0.15
Analytical solution(Sun et al., 2004)
HP1

0.1

0.05

trans-DCE
1,1-DCE00

10Concentration (mol / l)Concentration (mol / l)

(c)
0.08

(d)

VC

0.06

0.04

ETH
0.02

0
20

30

0

10

Distance (m)2030Distance (m)Figure 80 - Profiles after 200d for the PCE biodegradation problem. (a) Comparison of HP1 results for different
values of ωs with the analytical solution for PCE, (b) – (d) comparison of HP1 simulations with ωs = 0.1 with theanalytical solution for 7 components

Possible extensions of this application that can be readily considered with HP1 are adsorption, otherreaction networks (e.g., TNT degradation to produce 2ADNT, 4ADNT, and TAT; Šimůnek et al., 2005a),transient variably saturated flow, spatial variability in the degradation coefficients (e.g., an anaerobic,
high carbon zone and an aerobic, low carbon zone, Aziz et al., 2000), and subsurface heterogeneity ingeneral (including layering).

Geochemical Modelling (HPg) – Modelling first-order degradation
Learning outcome
Implement first-order degradation reactions - PCE
Implement a sequential degradation network – PCE->DCE->

Background
A contaminant Cont is assumed to be subject to first-order degradation:

R = µ wC

(37)

where R is the decay rate [MaL-3T-1] and µw is the first-order rate constant for solutes in the liquid phase[T-1].
First-order decay (Eq. (37)) in PHREEQC is modelled using the keywords RATES and KINETICS. The
script statements in RATES result in a quantity moles which are transferred to or taken from the aqueous
phase. Thus, the expression to solve with the script statements in PHREEQC is:

=
RM

dM
= µwϑCdt(38)

where RM is the mole transfer (mole per time unit), M is the amount of component M (mol) and t is time
(time unit).

Problem Definition
In this exercise, 1 mmol PCE is introduced in 1 kg of water. PCE decay is simulated for 100 days using
the parameters in Table 6. We simulate the degradation of PCE and the sequential TCE degradation.
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

179
ISC: PublicSession 5

Conceptual model
The problem involves following components:••••

Equilibrium reactions
One phases (aqueous)
Kinetic aqueous degradation network
Well-mixed

We model a homogeneous (aqueous) kinetic reaction network path. . The schematic representation
was already shown in Figure 75.

Hands-on – First-order decay of PCE
Step 1 – Create Case

Create a new case: Example K_PCE 01
Description: “First-order degradation network - PCE”
Step 2 – Define global variables, new components and rates

PCE and its degradation products are not present in the database. Instead of defining them via the Add
Information button in the Database tab, we will add the information in a separate fil and include
the file in the project. The reason for this is that we can include the same file in the reactive transport
problem.
Create a new file via the “Create New Input File” button on the tool strip
Save the new input file:
Go to parent folder
Create new folder: Scripts
Go to Scripts
Save file as: PCE – First order degradation network – Rates.hsc
This file will contain global variables for the rates, alfa and yield coefficients, the definitions of the
new solution master species and the rates.
Add a global variable for the degradation rate (rate1 with a value of 0.075)
Add SOLUTION_MASTER_SPECIESPce Pce 1 Pce 1
Add SOLUTION_SPECIESPce= Pce ; log_k 0

Add RATES

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

180
ISC: PublicKinetics
PCEdegrad-mybasic-startrate = parm(1)*tot("water")*mol("Pce")mole = -rate*time()save(mole)-end

Code 54- Script for modelling PCE degradation using first-order degradation kinetics

Code 54 shows the implementation of Eq. (38), specifically Line 4. parm(1) the first­order rateequation that will be defined in the kinetic block. rate has the units (1/d)*(kgw / system) *(mass/kgwater) = mass/d/system.
Go to Additions to database Panel
Go to HP Keywords -> Scripting -> Include$ -> Select File – Only name
Browse to the script folder and select the file
Add ..\Scripts\ before the file name

INCLUDE$ ..\Scripts\PCE - First order degradation network Rates.hsc
A relative path starting from the project folder is created by “..\”. As mentioned on page 72,
INCLUDE$ starts from the project folder to look for the specified file. In this particular example, youhave
Project folder\..\Scripts\PCE - First order degradation network - Rates.hsc

You could also select HP Keywords -> Scripting -> Include$ -> Select File – Fullpath and replace the part before \Scripts with ..\.
Alternatively, INCLUDEPATH can be used to explicitly define folders that will be searched for files
defined by INCLUDE$. Via HP Keywords -> Scripting -> INCLUDEPATH -> Browse, you
can select the folder, and then include only the filename:
INCLUDEPATH D:\HP Geochemistry\Course\Scripts
INCLUDE$ PCE - First order degradation network - Rates.hsc

Step 3 – Geochemical model

Go to Geochemistry Panel
Add SOLUTION 1 with a concentration of 1 mg/I of PCE
Add END (to avoid printing of initial solution in output files)
Add USE solution 1
Add KINETICS 1
The definition of the kinetic reactants will be done in a separate file – as such, it can be used inother projects (e.g. the reactive transport problem) .
Create a new file via the “Create New Input File” button on the tool strip
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

181
ISC: PublicSession 5
Save the new input file:In the folder: Scripts
Save file as: PCE – First order degradation network – Kinetics.hsc
Add the PCEdegrad reactant
PCEDegrad-formula Pce 1-parms @rate1-m0 0
Go to Geochemistry Panel
Add under KINETICS the INCLUDE$ statement for PCE – First order degradation
network – Kinetics.hsc
Add the statement for a calculation of 100 days (-steps)
Add INCREMENTAL_ REACTIONS true
Add PRINT; -reset false
Step 4 – Output and plotting

Add SELECTED_OUTPUT for concentration of PCE
Add SELECTED_GRAPH
Figure 1: Plot PCE concentrations versus time
Figure 2: Plot Log10 PCE concentration versus time with GNUPLOT command set logscale:
set logscale y

OutputFigure 81 shows the decrease of PCE as a function of time ; the exponential decrease with time is the
typical behaviour of a first-order decay reaction. Plotting the concentration on a log10 scale gives alinear line (Figure 82). The slope is equal to the first-order rate constant/2.3 (plotting it on a In scale gives
directly the first-order rate constant). This property of first-order degradation can be used to verify the
implementation.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

182
ISC: PublicKinetics
Figure 81 – PCE concentration as a function of time.

Figure 82 – PCE concentration as a function of time on a semi-log plot.

Hands-on – Sequential degradation step PCE -> TCE
We expand the previous input to the degradation product TCE.
Go to the file PCE – First order degradation network – Rates.hsc
Add TCE in SOLUTION_MASTER_SPECIES and SOLUTION_SPECIESDefine the global variable yieldcoef1 for the yield of TCE from PCE degradation
Go to the file PCE – First order degradation network – Kinetis.hsc
In the kinetic reactant PCEDegrad, add the stoichiometric coefficient (yieldcoef1) for TCE in –formulaWhen PCE degrades (negative mass transfer, thus removed from the solution), TCE is produced. TCE
should have a positive mass transfer. Because the rate equation gives a negative value, thestoichiometry of TCE should be negative:-formula Pce 1 Tce -@yieldcoef1
Go to Geochemistry Panel
Add TCE to SELECTED_OUTPUT and SELECTED_GRAPH

OutputFigure 83 shows the time evolution of PCE and TCE. Note that TCE reaches a constant value because
degradation of TCE is not implemented. The plot gives an indication that the production of TCE iscorrectly implemented; its final concentration is 0.79 mg/I which corresponds to the yield coefficient.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

183
ISC: PublicSession 5

Figure 83 - PCE and TCE concentrations as a function of time (stable TCE).

Sequential degradation step PCE -> TCE ->Adding TCE degradation is straightforward:
Go to the file PCE – First order degradation network – Rates.hsc
Define the global variable rate2 for the TCE first-order rate constant
Add TCEdegrad to RATES
Go to the file PCE – First order degradation network – Kinetics.hsc
Add TCEdegrad to KINETICS

OutputFigure 84 shows the time evolution of PCE and TCE concentrations.

Figure 84 – PCE and TCE concentrations as a function of time.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

184
ISC: PublicKinetics

Geochemical Modelling (HPg) – Modelling first-order degradation
network
Learning outcome
Implement a diverging kinetic reaction network
TCE - > cis DCE->trans DCE-> 1,1 DCE
Implement a converging kinetic reaction network
cis DCE-> VC

trans DCE-> VC

1,1 DCE-> VC

Problem Definition
In this exercise, 1 mmol PCE is introduced in 1 kg of water. PCE decay is simulated for 100 days using
the parameters in Table 6. We implement the full PCE degradation network as illustrated in Figure 79.

Hands-on – Full first-order degradation network
Go to the file PCE – First order degradation network – Rates.hsc
Define all necessary global variables (rate, yieldcoef and alfa)
Add cisDCE, transDCE, 1.1 DCE, VC and ETH in SOLUTION_MASTER_SPECIES andSOLUTION_SPECIES (Dcec, Dcet, Dcee, Vc, Eth)
Add the rate equations RATES (DCEcdegrad, DCEtdegrad, DCEedegrad, VCdegrad andETHdegrad)
Go to the file PCE – First order degradation network – Kinetics.hsc
Add the kinetic reactants to KINETICS
For TCEDegrad, stoichiometric coefficients for Dcec, Dcet, and Dcee are obtained bymultiplying the yield coefficient yieldcoef2 with, respectively alfa1, alfa2 and alfa3
using the inline expressions operator @<> . For example, for Dcec (in –formula for TCEDegrad)Dcec -@<@yieldcoef2*@alfa1>
Go to Geochemistry Panel
Add Dcec, Dcet, Dcee, Vc, Eth in SELECTED_OUTPUT
Add figure in SELECTED_GRAPH with the concentrations of Dcec, Dcet, Dcee, Vc, EthOutput
The time evolutions of PCE and TCE are not changed (Figure 84). Figure 85 shows the concentration of
the remaining components.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

185
ISC: PublicSession 5

Figure 85 – DCE, VC and Eth concentrations as a function of time.

Reactive Transport Modelling (HP1) – Leaching and biodegradationof PCE
Learning outcome
Incorporate a previously defined geochemical model (HPGeochemistry) in HPx
Implement a boundary condition problem

Problem Definition
In this example we will simulate the transport and degradation of PCE and its daughter products in a soilcolumn. Degradation not only occurs as sequential reactions, but also partly as parallel degradation
reactions (see Figure 79). Degradation coefficients, yield factors and distribution factors are given in
Table 6.
Saturated flow conditions in a 2.0 m long soil column are maintained for 150 days. The inflowing solution
contains only 10 mmol of PCE. Will a steady-state profile be obtained within this timeframe? At which
depths do different components obtain a maximum value?
Physical properties of the soil are a porosity of 0.5, the saturated hydraulic conductivity of 1 cm/day anda dispersivity of 10 cm. Other soil hydraulic parameters are irrelevant for saturated conditions.

Conceptual model
The problem involves following components:
•
•
•
•
•
•
•
•

Equilibrium reactions
One phases (aqueous)
Aqueous kinetic degradation network
Well-mixed
Homogeneous domain
1D flow
1D advective-dispersive solute transport
Constant boundary condition
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

186
ISC: PublicKinetics
Figure 86 shows a schematic representation of the conceptual model.

Figure 86 – Schematic representation of a 1D reactive transport problem involving a kinetic reaction network(boundary condition problem).

Hands-out: Leaching and biodegradation of PCE
Step 1 – Define main settings and water flow
Project Manager
Button: "New"Name: "T-K-PCE-01"
Description: "PCE degradation network "Button: "OK"General information
Set soil depth of 200 cm
Set final time of 150 days
Only print output at regular time intervals of 1 day
Ask for 5 print times
Water Flow
Set porosity and saturated hydraulic conductivity
Upper and lower boundary conditions: constant pressure head
Step 2 – Solute transport and geochemistry
Geochemistry
Go to Reactive Transport Panel
Select the phreeqc.dat database
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

187
ISC: PublicSession 5
Go to Selection of Components item
Add the seven components in the table Components from Text Input
Go to Additions to database Panel
Include the file with the global definitions, SOLUTION_MASTER_SPECIES and the RATES
definitions (as an absolute or relative path)
Go to Solution definitions Panel
Define the initial solution (pure water) as SOLUTION 1001
Define the upper boundary solution (10 mmol/l PCE) as SOLUTION 3001
Go to Geochemistry Panel
Add KINETICS 1001 and link it to material 1
Include the file with the kinetic reactant definitions (as an absolute and relative path)
As we are only interested in concentrations, the results can be viewed using the post-processing
dialogue windows in the HYDRUS graphical user interface. When you want to use Gnuplot for
plotting, you need to define SELECTED_OUTPUT and SELECTED_GRAPH.
Go to Output Panel
Add SELECTED_OUTPUT
Add SELECTED_GRAPH – include a Time series plot with the concentrations of all components(HP Keywords -> Output -> SELECTED_GRAPH -> Templates -> Time Series).
If you want to plot only the last node, use –fig_multi last. Otherwise, you can use –fig_multi all(see also section Output – Creating time series plot on page 153).
Button “Next”
Transport
Change dispersivity
Give solution number for upper boundary condition (3001)
Step 3 – Profile Information
Set initial pressure head to 0
Include an observation point every 50 cm
Save and Run

OutputSome selected profiles and times series are shown in Figure 87. Figure 88 shows the outflow
concentrations of PCE and its degradation products.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

188
ISC: PublicKinetics

Figure 87 – Profiles (top) and time series (bottom) of some selected degradation products from PCE for steadystate flow conditions.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

189
ISC: PublicSession 5
Figure 88 – Outflow concentrations of PCE and its degradation products during steady-state flow conditions.

Reactive Transport Modelling (HP1) – Leaching of initial PCE source
Learning outcome
Implement an initial value problem – part of the soil profile is contaminatedAsk output variables that are summed over different zones in the soil profile

Problem definition
Instead of the boundary condition problem in the previous example, we consider now an initial value
problem with a PCE source in the top 50 cm. We assume that the amount of non-aqueous (immobile)
PCE is 0.01 mol/ dm³ soil. Plot the amount of PCE in the top 50 cm as a function of time.
Knauss et al. (2000) reported following partial molal thermodynamic quantities for the PCE liquidaqueous dissolution reaction:∆ r H º =− 1790 J/mol(39)∆ r S º =− 59 J/mol K(40)∆ rC °p =− 354.6 J/mol(41)

where R is the gas constant (8.31441 J/mol), and ∆rHº, ∆rSº, and ∆rCpº the change of entropy, enthalpy
and heat capacity of reaction at standard conditions (T0 = 298 K). The following equation is used in
PHREEQC to define the temperature dependence of solubility constants (Parkhurst and Appelo 1999):

log( K ) = A1 + A2T +A3
A+ A4 log(T ) + 52
T
T

(42)

Temperature dependence of the solubility coefficient can be described with a three-term extrapolation
expression with the coefficients A1, A3, and A4 using following relations (Puigdomènech et al. 1997):1
 ∆ r S ° − ∆ rC °p ( ln T 0 + 1) 
2.303R(43)

A3 = −

1( ∆ r H ° − ∆ rC °pT 0 )2.303R

(44)

=
A4

1
∆ rC °p ln (10 )2.303R(45)=
A1

The values for A1, A3 and A4 are -127.10, 5612 and 42.64.

Conceptual model
The problem involves following components:
•
•
•
•
•
•
•

Equilibrium reactions
Two phases (aqueous, solid)
Equilibrium dissolution of immobile phase
Aqueous kinetic degradation network
Well-mixed
Homogeneous domain
1D flow
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

190
ISC: PublicKinetics••

1D advective-dispersive solute transport
Initial value condition

Figure 89 shows a schematic representation of the conceptual model.

Figure 89 – Schematic representation of a 1D reactive transport problem involving a kinetic reaction network(Initial value problem).

Hands-on: Leaching of a PCE source
Step 1 – Define main settings and water flow
Project Manager
Select: "T-K-PCE-01"
Button: “Copy”
Name: "T-K-PCE-02”
Description: "PCE degradation network – Leaching of source"
Button: "OK"Open: “T-K-PCE-02”
General information
Increase the number of materials to 2. Material 1 is the top 50 cm with the PCE contamination.
Decrease maximum time step to 0.1 (controlling split-operator error, see Jacques et al. (2006))
Step 2 – Solute transport and geochemistry
Geochemistry
Open the file PCE – First order degradation network – Rates.hsc
Add a new phase that controls the solubility of PCE using the keyword PHASES
PCE_lqPce = Pce
-analytical_expression -127.10 0 5612 42.64 0
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

191
ISC: PublicSession 5
Go to Solution Definition Panel
Redefine SOLUTION 1001 to be initially in equilibrium with the PCE source (contaminated zone)
Define SOLUTION 1002 – pure water for uncontaminated zone
Redefine SOLUTION 3001 as a pure water
Go to Geochemistry Panel
Add KINETICS for material 2
Alternatively, you can define the same model for kinetic reactants for the whole soil as:
KINETICS 1-@_NODESINCLUDE$ ..\Scripts\PCE - First order degradation network- Kinetics.hsc
Add Pce_lq in EQUILIBRIUM_PHASES 1001 linked to material 1 (amount: 0.01 mol/dm³)

Output – Integrated variables
To have the total amount of PCE in the contaminated zone, we need to integrate the amounts in the
aqueous phase and the source with depth. The keyword REACTIVE_TRANSPORT_OUTPUT of HPx
defines over which materials, depths or nodes such an integration is required. When an integration isrequired over a given material type, the identifier –horizon is used. Code 55 shows the syntax of –horizon:•••

List of numbers: a single integer, or a list of different integers that will be grouped indicating thematerial numbers of the nodes that will be integrated.
Name: the name of this particular integrated output variables
Mult: integer indicating the time intervals at which the integration has to be executed:
o If not defined or zero, the integration is done at the print times of the profiles as definedin the HYDRUS graphical user interface;o If larger than zero, the integration is done at the time steps defined by Time interval
for regular print times defined in the Print Information dialogue window of
the HYDRUS graphical user interface, multiplied by the integer value.

These three inputs can be repeated in the same –horizon line for integration over different materials
or different –horizons can be defined. Integration is done for all (relevant) variables defined by
SELECTED_OUTPUT either by summing (e.g. for the amounts of a mineral or a kinetic reactant) or
weighted by water content (e.g. aqueous concentrations). Integration for the output variables definedby USER_PUNCH is controlled by a list of integers after the identifier –integralsum:ooo-99: No sum is made (-1E8 is printed in the output file)
1: a simple sum
2: sum weighted by water content
REACTIVE_TRANSPORT_OUTPUT-horizon <:ints:list of numbers> <:str:name> [t<:int:mult>]-horizon 1 "Source" t 1
Code 55 – Syntax of –horizon in REACTIVE_TRANSPORT_OUTPUT.

The integration can also be defined over depth intervals (-depth), nodes (-node) or (mass balance)
layers (-layer).
Line 3 in Code 55 indicates that an integration will be done over the nodes of material 1 at every
regular print time step. Output variables and files will contain the word “Source”.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

192
ISC: PublicKinetics
Go to Output Panel
Add Pce_lq to SELECTED_OUTPUT (-equilibrium_phases)
Add USER_PUNCH to calculate the total amount of PCE (i.e. the sum in the liquid and in the solid
phase)
USER_PUNCH-headings Pcetot-integralsum 1-mybasic-startpce = totmole("Pce") + equi("Pce_lq")
punch(pce)-end
The identifier –integralsum 1 indicates that the first variable (Pcetot) will be integrated as asimple sum.
totmole(“Pce”) returns the moles of a component in solution. Note the difference with tot();
tot() returns the total concentration of a component in solution, whereas totmole() returns
the number of moles in solution.
Add REACTIVE_TRANSPORT_OUTPUT
REACTIVE_TRANSPORT_OUTPUT-horizon 1 "Source" t 1-horizon 2 "Soil" t 1-horizon 1-2 "Total" t 1
The first –horizon will result in output and plots of the totals of PCE and its degradation productsin the source zone; the second one in the pristine soil under the source zone. The third one will
integrate the variables over the complete soil profile. Integration is done at every regular print
time.
Next
Transport
Adapt dispersivity of the second material
Step 3 – Profile information
Define the second material (50-200 cm) via Profile Information or Profile summary dialogue
windows
Define number 1002 to material 2 via Profile Information (Conditions->Initial Conditions>Concentrations; Concentration number 1) or Profile summary dialogue window.
Save and Run

OutputPHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

193
ISC: PublicSession 5

Figure 90 – Profiles of aqueous PCE (left) and immobile PCE (right) during leaching of an initially contaminated soil
profile.

Figure 91 – Outflow concentrations of PCE and its degradation products during leaching of an initiallycontaminated soil profile.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

194
ISC: PublicKinetics

Figure 92 – Total amount of PCE in the aqueous phase (top left), immobile phase (top right) and summed amount(bottom) in the source zone (0-50 cm), pristine soil (50-200 cm) and total soil profile during leaching of an initiallycontaminated soil profile.

Reactive Transport Modelling (HP1) – Leaching of initial PCE source –Source trapped in immobile water region
Learning outcome
Implement a reactive transport problem for a dual-porosity (Mobile-Immobile) solute transportmodel

Problem Definition
It is assumed that the solid PCE phase is trapped in immobile zones in the soil profile. Degradation
reactions occurs equally in the mobile and immobile phase. Assume following properties: immobile
water content: 0.1, first-order exchange coefficient: 0.01 day-1.

Conceptual model
The problem involves following components:
•
•
•
•
•

Equilibrium reactions
Two phases (aqueous, solid)
Equilibrium dissolution of immobile phase
Aqueous kinetic degradation network
Well-mixed
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

195
ISC: PublicSession 5
•
•
•
•
•

Homogeneous domain
1D flow
1D advective-dispersive solute transport
Physical non-equilibrium transport – Mobile-Immobile
Initial value condition

Figure 93 shows a schematic representation of the conceptual model.

Figure 93 – Schematic representation of a 1D reactive transport problem involving a kinetic reaction network(Initial value problem, mobile-immobile solute transport).

Hands-on: Leaching of a PCE source trapped in immobile region
Project Manager
Select: "T-K-PCE-02"
Button: “Copy”
Name: "T-K-PCE-03”
Description: "PCE degradation network – Leaching of source – MIM”
Button: "OK"Open: “T-K-PCE-02”

General information
Increase final time to 500 days.
Ask for 10 number of print times
Step 2 – Solute transport and geochemistryTransport – General Information
Select Dual Porosity Model
Button: Next
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

196
ISC: PublicKinetics
Geochemistry
Go to Solution Definitions Panel
Redefine SOLUTION 1001 (pure water): mobile region contaminated zone
Add SOLUTION 2001(equilibrium with Pce_lq): immobile region contaminated zone
Add SOLUTION 2002 (pure water): immobile region pristine soil
Go to Geochemistry Panel
Add KINETICS 2001 and 2002 for immobile regions. To indicate that the kinetic block is
assigned to the immobile region of a given material type, add im after -material:-material im 2
Change EQUILIBRIUM_PHASES 1001 to 2001, and link it to the immobile region
Go to Output Panel
Add –kinetic_reactants PceDegrad to SELECTED_OUTPUT. –kinetic_reactants will
print the amount of the kinetic reactant and the change in the kinetic reactant (in HPx in
mol/T/dm³, in HPGeochemistry/PHREEQC in mol/time step).
To USER_PUNCH
Add –totalsum 1 before –start
For each the output variables, you can define if a total (of mobile and immobile region) has to be
calculated and how it is calculated. Available options are:
0 No sum of mobile and immobile region is calculated
1 Simple sum of mobile and immobile region is calculated (e.g. for minerals)
2 Sum weighted by water contents (e.g. for concentrations)
Button: Next
Transport
Change mass transfer coefficient and immobile water content.
Step 3 – Profile Information
Go to Soil Profile – Summary
Adapt the numbers in the column Im. C. Comp. to 2001 (where material Mat is 1) and to 2002
(where material Mat is 2)
Save and Run

OutputSome selected results are shown in Figure 94 and Figure 95.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

197
ISC: PublicSession 5

Figure 94 – Profiles of PCE in the mobile region (top left), in the immobile region (top right), the averaged
concentration over both regions (bottom left), and the PCE in the immobile phase in the immobile region (bottom
right)

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

198
ISC: PublicKinetics

Figure 95 – Concentration at the bottom node in the mobile (top) and immobile (bottom) region.

Figure 96 plots the summed kinetic reactant PceDegrad. As its initial value is zero, it shows the total
amount of PCE degraded in a given soil compartment. The left bottom plot indicates that at 500 days,
about 0.5 mol/cm² is degraded in the complete soil profile, which is about equal to the initial amountpresent in the soil. About 75% of the PCE degrades in the contaminated zone. At the soil profile scale,
about 50% of the PCE degrades in the immobile zone (mostly in the contaminated zone). Checking PostProcess, Solute Transport – Actual and Cumulative Boundary fluxes – Solution 3 (PCE), you can see howmuch PCE has leached from the soil (about 2.2 mol/cm²) (see right bottom plot).

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

199
ISC: PublicSession 5

Figure 96 – Cumulative PCE degradation in different soil compartments during leaching of PCE trapped in theimmobile soil region (top plots and bottom left). The bottom right plot shows the cumulative leaching of PCE fromthe soil profile.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

200
ISC: PublicReferences
Aagaard, P. and Helgeson, H.C. (1982) Thermodynamic and kinetic constraints on reaction rates amongminerals and aqueous solutions, 1. Theoretical considerations. American Journal of Science 282, 237285.
Appelo, C.A.J. and Parkhurst, D.L. (2002) CALCULATING CATION EXCHANGE WITH PHREEQC (Version 2).
http://www.hydrochemistry.eu/pub/ap_pa02.pdf.
Appelo, C.A.J. and Postma, D. (2005) Geochemistry, groundwater and pollution, 2nd edition, A.A.
Balkema Publishers, Amsterdam, The Netherlands.
Bao, C., Li, L., Shi, Y.N. and Duffy, C. (2017) Understanding watershed hydrogeochemistry: 1. Developmentof RT-Flux-PIHM. Water Resources Research 53(3), 2328-2345.
Batlle-Aguilar, J., Brovelli, A., Porporato, A. and Barry, D.A. (2011) Modelling soil carbon and nitrogen
cycles during land use change : A Review. Agronomy Sustainable Development. 31(0), 251-274.
Bedekar, V., Morway, E.D., Langevin, C.D. and Tonkin, m. (2016) MT3D-USGS version 1: A U.S. Geological
Survey release of MT3DMS updated with new and expanded transport capabilities for use with
MODFLOW:U.S.
Geological
SurveyTechniquesandMethods
6-A53,
69p.,
http://dx.doi.org/10.3133/tm6A53.
Beegum, S., Simunek, J., Szymkiewicz, A., Sudheer, K.P. and Nambi, I.M. (2019) Implementation of Solute
Transport in the Vadose Zone into the "HYDRUS Package for MODFLOW". Groundwater 57(3), 392-408.
Bekins, B.A., Warren, E. and Godsy, E.M. (1998) A comparison of zero-order, first-order, and Monodbiotransformation models. Ground Water 36(2), 261-268.
Bildstein, O., Claret, F. and Frugier, P. (2019) RTM for Waste Repositories. Reviews in Mineralogy andGeochemistry 85(1), 419-457.
Bond, W.J. (1995) On the Rothmund-Kornfeld description of cation exchange. Soil Science Society ofAmerica Journal 59(2), 436-443.
Brantley, S., Goldhaber, M.B. and Ragnarsdottir, K.V. (2007) Crossing disciplines and scales to understandthe critical zone, pp. 307-314.
Brantley, S.L. (2003) Surface and Ground Water, Weathering, and Soils Drever, J.I., Holland, H.D. andTurekian, K.K. (eds), pp. 73-119, Volume 5, Elsevier.
Casey, F.X.M. and Simunek, J. (2001) Inverse Analyses of Transport of Chlorinated Hydrocarbons Subject
to Sequential Transformation Reactions, pp. 1354-1360.
Charlton, S.R. and Parkhurst, D.L. (2011) Modules based on the geochemical model PHREEQC for use inscripting and programming languages. Computers & Geosciences 37(10), 1653-1663.
Deutsch, W.J. and Siegel, R. (1997) Groundwater Geochemistry: Fundamentals and Applications to
Contamination, Lewis Publishers.
Dove, P.M. (1995) Chemical Weathering Rates of Silicate Minerals. White, A.F. and Brantley, S. (eds), pp.
235-290, De Gruyter.
Finke, P.A. and Hutson, J.L. (2008) Modelling soil genesis in calcareous loess. Geoderma 145(3-4), 462479.
Gautier, J.M., Oelkers, E.H. and Schott, J. (1994) Experimental-Study of K-Feldspar Dissolution Rates as a
Function of Chemical Affinity at 150-Degrees-C and Ph-9. Geochimica et Cosmochimica Acta 58(21),4549-4560.
Goldberg, S., Criscenti, L.J., Turner, D.R., Davis, J.A. and Cantrell, K.J. (2007) Adsorption - Desorption
processes in subsurface reactive transport modeling. Vadose Zone Journal 6(3), 407-435.
He, W., Beyer, C., Fleckenstein, J.H., Jang, E., Kolditz, O., Naumov, D. and Kalbacher, T. (2015) A
parallelization scheme to simulate reactive transport in the subsurface environment with OGS#IPhreeqc
5.5.7-3.1.2. Geosci. Model Dev. 8(10), 3333-3348.
Huertas, F.J., Caballero, E., de Cisneros, C.J., Huertas, F. and Linares, J. (2001) Kinetics of montmorillonite
dissolution in granitic solutions. Applied Geochemistry 16(4), 397-407.

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

201
ISC: Public

Jacques, D., Simunek, J., Mallants, D. and van Genuchten, M.T. (2018) The HPx software formulticomponent reactive transport during variably-saturated flow: Recent developments andapplications. JOURNAL OF HYDROLOGY AND HYDROMECHANICS 66(2), 211-226.
Jacques, D., Šimůnek, J., Mallants, D. and van Genuchten, M.T. (2006) Operator-splitting errors in coupled
reactive transport codes for transient variably saturated flow and contaminant transport in layered soil
profiles. Journal of Contaminant Hydrology 88, 197-218.
Jacques, D., Šimůnek, J., Mallants, D. and van Genuchten, M.T. (2008a) Modeling Coupled Hydrologic
and Chemical Processes: Long-Term Uranium Transport following Phosphorus Fertilization. Vadose Zone
Journal 7(2), 698-711.
Jacques, D., Šimůnek, J., Mallants, D. and van Genuchten, M.T. (2008b) Modelling coupled water flow,
solute transport and geochemical reactions affecting heavy metal migration in a podzol soil. Geoderma
145(3-4), 449-461.
Jang, E., Boog, J., He, W. and Kalbacher, T. (2018) Computational Hydrology III: OGS-PHREEQC coupled
reactive transport modeling. ISBN: 978-3-319-67152-9, Springer International Publishing, DOI:
10.1007/978-3-319-67153-6
Jansen, B. (2003) The mobility of aluminium, iron and organic matter in acidic sandy soils.
Jia, M., Jacques, D., Gérard, F., Su, D., Mayer, K.U. and Šimůnek, J. (2019) A benchmark for soil organic
matter degradation under variably saturated flow conditions. Computational Geosciences.
Johnson, J.W., Oelkers, E.H. and Helgeson, H.C. (1992) SUPCRT92: A software package for calculating thestandadard molal thermodynamic properties of minerals, gases, aqueous species and reactions from 1
to 5000 bar and 0 to 1000>C, pp. 899-947.
Keyvanshokouhi, S., Cornu, S., Samouëlian, A. and Finke, P. (2016) Evaluating SoilGen2 as a tool for
projecting soil evolution induced by global change. Science of the Total Environment 571, 110-123.
Kinniburgh, D.G. and Cooper, D.M. (2009) PhreePlot: Creating graphical output with PHREEQC.
Knauss, K.G., Dibley, M.J., Leif, R.N., Mew, D.A. and Aines, R.D. (2000) The aqueous solubility oftrichloroethene (TCE) and tetrachloroethene (PCE) as a function of temperature, pp. 501-512.
Leal, A.M.M., Kulik, D.A., Smith, W.R. and Saar, M.O. (2017) An overview of computational methods forchemical equilibrium and kinetic calculations for geochemical and reactive transport modeling. Pure andApplied Chemistry 89(5), 597-643.
Leterme, B., Blanc, P. and Jacques, D. (2014) A reactive transport model for mercury fate in soil—
application to different anthropogenic pollution sources. Environmental Science and Pollution Research,
12279-12293.
Leterme, B., Geukens, L., Jacques, D., Heyvaert, V.M.A., Huysmans, M., Neyens, C., Smolders, E., Springael,D., van Wesemael, B. and Walstra, J. (2019) RESPONSE project: Reactive transport of point source
contamination in soils and groundwater, pp. 54-61.
Levenspiel, O. (1999) Chemical reaction engineering.
Li, L., Maher, K., Navarre-Sitchler, A., Druhan, J., Meile, C., Lawrence, C., Moore, J., Perdrial, J., Sullivan, P.,
Thompson, A., Jin, L., Bolton, E.W., Brantley, S.L., Dietrich, W.E., Mayer, K.U., Steefel, C.I., Valocchi, A.,
Zachara, J., Kocar, B., McIntosh, J., Tutolo, B.M., Kumar, M., Sonnenthal, E., Bao, C. and Beisman, J. (2017)
Expanding the role of reactive transport models in critical zone processes. Earth-Science Reviews 165,
280-301.
Makselon, J., Zhou, D., Engelhardt, I., Jacques, D. and Klumpp, E. (2017) Experimental and NumericalInvestigations of Silver Nanoparticle Transport under Variable Flow and Ionic Strength in Soil.
Environmental Science & Technology 51(4), 2096-2104.
Mayer, K.U., Alt-Epping, P., Jacques, D., Arora, B. and Steefel, C.I. (2015) Benchmark problems for reactive
transport modeling of the generation and attenuation of acid rock drainage. Computational Geosciences
19(3), 599-611.
Mayer, K.U., Frind, E.O. and Blowes, D.W. (2002) Multicomponent reactive transport modeling in variably
saturated porous media using a generalized formulation for kinetically controlled reactions. Water
Resources Research 38, 1174, doi:1110.1029/2001WR000862.
Meeussen, J.C.L. (2003) ORCHESTRA: An object-oriented framework for implementing chemical
equilibrium models. Environmental Science & Technology 37(6), 1175-1182.
Meysman, F.J.R., Middelburg, J.J., Herman, P.M.J. and Heip, C.H.R. (2003a) Reactive transport in surface
sediments. I. Model complexity and software quality, pp. 291-300.
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

202
ISC: Public

Meysman, F.J.R., Middelburg, J.J., Herman, P.M.J. and Heip, C.H.R. (2003b) Reactive transport in surface
sediments. II. Media: an object-oriented problem-solving environment for early diagenesis. Computers& Geosciences 29(3), 301-318.
Nagy, K.L., Blum, A.E. and Lasaga, A.C. (1991) Dissolution and precipitation kinetics of kaolinite at 80
degrees C and pH 3; the dependence on solution saturation state, pp. 649-686.
Nardi, A., Idiart, A., Trinchero, P., De Vries, L.M. and Molinero, J. (2014) Interface COMSOL-PHREEQC (iCP),
an efﬁcient numerical framework for the solution of coupled multiphysics and geochemistry. Computers& Geosciences 69, 10–21.
Oostrom, M., Mehmani, Y., Romero-Gomez, P., Tang, Y., Liu, H., Yoon, H., Kang, Q., Joekar-Niasar, V.,
Balhoff, M.T., Dewers, T., Tartakovsky, G.D., Leist, E.A., Hess, N.J., Perkins, W.A., Rakowski, C.L., Richmond,
M.C., Serkowski, J.A., Werth, C.J., Valocchi, A.J., Wietsma, T.W. and Zhang, C. (2014) Pore-scale andcontinuum simulations of solute transport micromodel benchmark experiments. Computational
Geosciences, 1-23.
Pabalan, R.T. and Bertetti, F.P. (2001) Cation-exchange properties of natural zeolites. Reviews in
Mineralogy and Geochemistry 45, 453-518.
Parkhurst, D.L. and Appelo, C.A.J. (1999) User's guide to PHREEQC (version 2) - A computer program forspeciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations.
Parkhurst, D.L. and Appelo, C.A.J. (2013a) Description of input and examples for PHREEQC version 3—Acomputer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical
calculations.
Parkhurst, D.L. and Appelo, C.A.J. (2013b) Description of Input and Examples for PHREEQC Version 3 —A Computer Program for Speciation , Batch-Reaction , One-Dimensional Transport , and Inverse
Geochemical Calculations
Parkhurst, D.L. and Wissmeier, L. (2015) PhreeqcRM: A reaction module for transport simulators basedon the geochemical model PHREEQC. Advances in Water Resources 83(0), 176-189.
Patel, R., Perko, J., De Schutter, G., Van Breugel, K. and Ye, G. (2014) A versatile pore-scale
multicomponent reactive transport approach based on lattice Boltzmann method: Application to
portlandite dissolution. Physics and Chemistry of the Earth 70-71, 127-137.
Patel, R.A., Perko, J., Jacques, D., De Schutter, G., Ye, G. and Van Breugel, K. (2018) A three-dimensionallattice Boltzmann method based reactive transport model to simulate changes in cement paste
microstructure due to calcium leaching. Construction and Building Materials 166, 158-170.
Payne, T.E., Davis, J.A. and Waite, T.D. (1996) Uranium adsorption on ferrihydrite - Effects of phosphateand humic acid. Radiochimica Acta 74, 239-243.
Perko, J., Mayer, K.U., Kosakowski, G., De Windt, L., Govaerts, J., Jacques, D., Su, D. and Meeussen, J.C.L.
(2015) Decalcification of cracked cement structures. Computational Geosciences 19, 673-693.
Prommer, H., Barry, D.A. and Zheng, C. (2003) MODFLOW/MT3DMS-based reactive multicomponent
transport modeling. Ground Water 41(2), 247-257.
Puigdomènech, I., Rard, J.A., Plyasunov, A.V. and Grenthe, I. (1997). Grenthe, I. and PuigdomSnech, I.
(eds), pp. 427-493, OECD Nuclear Chemistry, Paris.
Schaerlaekens, J., Mallants, D., Simunek, J., van Genuchten, M.T. and Feyen, J. (1999) Numerical simulationof transport and sequential biodegradation of chlorinated aliphatic hydrocarbons using CHAIN_2D, pp.
2847-2859.
Seo, H.S., Simunek, J. and Poeter, E.P. (2007) Documentation of the HYDRUS Package for MODFLOW2000, the U.S. Geological Survey Modular Ground-Water Model. Colorado School of Mines, Golden, CO.
Seuntjens, P., Mallants, D., Simunek, J., Patyn, J. and Jacques, D. (2002) Sensitivity analysis of physical andchemical properties affecting field-scale cadmium transport in a heterogeneous soil profile. Journal of
Hydrology 264(1-4), 185-200.
Simunek, J., Jacques, D., Langergraber, G., Bradford, S.A., Sejna, M. and van Genuchten, M.T. (2013)
Numerical Modeling of Contaminant Transport Using HYDRUS and its Specialized Modules. Journal of
the Indian Institute of Science 93, 265-284.
Šimůnek, J., Jacques, D. and Šejna, M. (2013) HP2 / 3 : Extensions of the HP1 Reactive Transport Code to
Two and Three Dimensions, pp. 345-354.
PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

203
ISC: PublicŠimůnek, J. and Suarez, D.L. (1994) Two-dimensional transport model for variably saturated porous
media with major ion chemistry. Water Resources Research 30(4), 1115-1133.
Šimůnek, J., van Genuchten, M.T. and Šejna, M. (2016) Recent Developments and Applications of the
HYDRUS Computer Software Packages. Vadose Zone Journal 15(7).
Simunek, J., van Genuchten, M.T. and Senja, M. (2008) Development and Applications of the HYDRUS
and STANMOD Software Packages and Related Codes. Vadose Zone Journal 7, 587-600.
Steefel, C.I., Appelo, C.A.J., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., Lagneau, V., Lichtner, P.C.,
Mayer, K.U., Meeussen, J.C.L., Molins, S., Moulton, D., Shao, H., Šimůnek, J., Spycher, N., Yabusaki, S.B.
and Yeh, G.T. (2015) Reactive transport codes for subsurface environmental simulation. Computational
Geosciences 19, 445-478.
Steefel, C.I., DePaolo, D.J. and Lichtner, P.C. (2005) Reactive transport modeling: An essential tool and a
new research approach for the Earth sciences. Earth and Planetary Science Letters 240, 539-558.
Sun, Y.W., Lu, X.J., Petersen, J.N. and Buscheck, T.A. (2004) An analytical solution of tetrachloroethylenetransport and biodegradation. Transport in Porous Media 55(3), 301-308.
Vandenhove, H., Hurtgen, C. and Payne, T.E. (2011) Uranium: Radionuclides, Wiley.
Vandenhove, H., Van Hees, M., Wouters, K. and Wannijn, J. (2007) Can we predict uranium bioavailability
based on soil parameters? Part 1: Effect of soil parameters on soil solution uranium concentration.
Environ Pollut 145(2), 587-595.
Walther, J.V. (1996) Relation between rates of aluminosilicate mineral dissolution, pH, temperature, andsurface charge. American Journal of Science 296(7), 693-728.
Wissmeier, L. and Barry, D.a. (2011) Simulation tool for variably saturated flow with comprehensive
geochemical reactions in two- and three-dimensional domains. Environmental Modelling & Software
26(2), 210-218.
Xie, M., Mayer, K.U., Claret, F., Alt-Epping, P., Jacques, D., Steefel, C., Chiaberge, C. and Šimůnek, J. (2015)
Implementation and evaluation of permeability-porosity and tortuosity-porosity relationships linked tomineral dissolution-precipitation. Computational Geosciences 19, 655-671.
Zhu, C. and Anderson, G. (2002) Environmental application of geochemical modelling, Cambridge
University Press.

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

204
ISC: PublicIndexArrhenius equation, 164

B

-, 32
B-dot model, 24C- layer, 192- nodes, 192

"
", 32

#
#, 32
#{, 42#}, 42#endregion, 42#region, 42&

charge, 49
Chemical state, 17
Comments, 32
Component, 16
Concentration Units, 38
COUNTER, 71
COUNTERINCR, 71
Curve, 64Curve options, 76, 83
dt, 76, 84
i, 141k, 76
ls, 76, 83
style, 76vk, 76, 83&, 69, 75

;
;, 32

@
@, 63, 68@_NODES, 192@<>, 81{{}, 63<<>, 76A-active, 60, 141-activities, 59
Activity correction
aqueous species, 24
B-dot model, 24
Davies equation, 35
Extended Debye-Hückel model, 24
WATEQ equation, 36
Add components and species, 173

D
Database Manager, 28
Davies activity correction equation, 35
Default database for project group, 117-depth, 192
Dialogue window
Default database for project group, 117
Dialogue windows
Graph Settings, 63
Thermodynamic Data (adding), 173

E
Edit Graph Panel, 64-end, 74, 161
END, 33
Enhanced text format, 83, 147-equilibrium_phases, 59
EQUILIBRIUM_PHASES, 57
evaporation, 145
Evaporation reaction path with fixed activity, 145
EXCHANGE, 134
Exchange isotherm, 131, 140
EXCHANGE_MASTER_SPECIES, 132
EXCHANGE_SPECIES, 133
Exporting graphs, 64
Extended Debye-Hückel model, 24

PHREEQC-HPx Course, Prague, September 3-4, 2020
SCK CEN/39560003

205
ISC: Public

FFe(3), 66
Ferric iron, 66-fig_curve, 154-fig_curvt, 154-fig_end, 75-fig_multi, 154-fig_name, 74-fig_nr, 74-fig_scurve, 154-fig_start, 75-fig_sy, 147-fig_syt, 147-fig_type, 153-fig_x, 74-fig_xtr, 74-fig_y, 74-fig_yt, 74-file, 59
First-order degradation, 179
Fixing the pH, 57
Folding, 42
Forms
Aqueous Components (adding), 173-formula, 162, 175

G
Gaines-Thomas convention, 131, 133
Gapon convention, 131
Geochemical system, 16
gfw(), 168
Global Definition Panel
Global Variable, 70
Global variable Form, 70
Edit, 70
Multiple cases, 71, 170
Values, 71
GLOBAL_VARIABLE, 68
Gnuplot, 75
blanc lines, 141
Data file structure, 62
discontinuity in a curve, 141
Enhanced text format, 63, 83, 147
Multi-dataset file, 141
set key, 64
set logscale, 84, 182
set style line, 83, 84
set title, 83, 169
set xrange, 83
set xtics ... format, 121
set yrange, 83
set ytics, 83, 84
Graph Settings dialogue window, 63
Graph workspace, 64
GraphsClipboard, 64Create from data file, 62Curve, 64
D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

Edit graph, 64
Edit Graph Panel, 64
Gnuplot, 75
Graph Settings dialogue window, 63
Ribbon, 62
SELECTED_GRAPH, 74
Time Series Plot, 153
Workspace, 64

H
Heterogeneous kinetic reaction path, 166
Homogeneous kinetic reaction network path, 180
Homogeneous kinetic reaction path, 172-horizon, 192
How ToAdd a database, 30
Add a database group, 30
Add a project, 31
HPGeochemistry
Main Screen, 27

I
include, 73
INCLUDE$, 72, 181
INCLUDEPATH, 181
INCREMENTAL_REACTIONS, 58
inline expression, 185
inline expression operator, 81
Inline variables, 68@_NODES, 192
GLOBAL_VARIABLE, 79
inline expression, 81
Input Work Space, 41
Global Definition Panel, 70-integralsum, 192

K
Keyword
COUNTER, 71
COUNTERINCR, 71
END, 33
EQUILIBRIUM_PHASES, 57
EXCHANGE, 134
EXCHANGE_MASTER_SPECIES, 132
EXCHANGE_SPECIES, 133
GLOBAL_VARIABLE, 68
INCLUDE$, 72, 181
INCLUDEPATH, 181
INCREMENTAL_REACTIONS, 58
KINETICS, 161
MIX, 54
PHASES, 37
PUNCH, 141
RATES, 160
REACTION, 58
REACTIVE_TRANSPORT_OUTPUT, 192
REPEAT, 73
206
ISC: Public

SAVE, 135
SELECTED_GRAPH, 74
SELECTED_OUTPUT, 59
SOLUTION, 38
SOLUTION_MASTER_SPECIES, 34
SOLUTION_SPECIES, 35
USE, 33, 135
USER_PUNCH, 73-kinetic_reactants, 197
KINETICS, 161-formula, 162, 175-m0, 161-parms, 162-steps, 162

L

Saturation index, 26
PCE, 176
phase_formula(), 168
PHASES, 37
Phenol, 172
PhreePlot, 66
Phreeqc.out, 45-plot_titles, 123-plot_units, 123
Primary species, 16
PRINT-selected_output, 60
Project Manager, 30-punch, 60
PUNCH, 141
punch), 74

Qline continuation, 32
line separator, 32Mm(), 168-m0, 161
Manage databases, 28Manage Projects, 30
Manager
Database, 28
Projects, 30
Mass action law, 22
Master species, 16
Mercury, 18
Mineral kinetics
pH dependency, 165
Reactive surface area, 165
Temperature dependency, 164
Transition State Theory, 163
Mineral Kinetics, 162
MIX, 54
Mobile-immobile solute transport model, 195-molalities, 59
Monod kinetics, 171
Multi-Dataset File, 141
MuParserX, 69-mybasic, 74, 161
MyBasic, 74

O
Output file
File Structure, 62Graphs, 62
Output file workspace, 62Pparm(), 168-parms, 162
parms(), 162
Partial pressure of a gas, 26

Quartz, 166

Rrange of integers, 32
Rate EquationsFirst-order degradation, 179
Monod kinetics, 171
Transition State Theory, 163
RATE Examples
PCE biodegradation – First-order kinetics, 181
Phenol biodegradation - Monod, 174
Quartz dissolution, 168
RATES, 160-end, 161
-mybasic, 161-start, 161
REACTION, 58
Reaction path
Cation exchange isotherm, 140
Evaporation, 145
Heterogeneous kinetic reaction path, 166
Homogeneous kinetic reaction network path, 180
Homogeneous kinetic reaction path, 172
Sliding activity (pH), 66
Sliding concentration with fixed activity, 56
Reactive surface area, 165
REACTIVE_TRANSPORT_OUTPUT, 192- depth, 192-horizon, 192-layer, 192-nodes, 192
REPEAT, 73
include, 73
Ribbons
File Structure, 62Graph Export, 64
Graphs, 62
Rothmund-Kornfeld cation exchange model, 134

PHREEQC-HPx Course, Prague, September 3-4, 2020

SCK CEN/39560003

207
ISC: Public

S
Saturation index, 25
Saturation state, 25-saturation_indices, 169
SAVE, 135
save(), 161
Script functions
gfw(), 168m(), 168parm(), 168
parms(), 162
phase_formula(), 168
save(), 161
sr(), 168time(), 161
tot(), 169
total_time(), 169, 175
totmole(), 193
SELECTED_ GRAPH, 154-fig_curvet, 154-fig_scurve, 154-fig_syt, 147
SELECTED_GRAPH, 74&mol, 75
&none, 75&si, 75
Curve options, 76, 83, 141-fig_end, 75-fig_multi, 154-fig_name, 74-fig_nr, 74-fig_start, 75-fig_sy, 147-fig_type, 153-fig_x, 74-fig_xt, 74-fig_y, 74-fig_yt, 74
Multi-dataset file, 141
Time Series Plot, 153-selected_output, 60
SELECTED_OUTPUT, 59-active, 60, 141-activities, 59-equilibrium_phases, 59-file, 59-kinetic_reactants, 197-molalities, 59-saturation_indices, 169-total, 59, 60
Settings
Databases, 117

D. Jacques (djacques@sckcen.be)
SCK CEN/39560003

Defaults for HP1 Project Groups, 117
Sliding activity reaction path (pH), 66
Sliding concentration reaction path with fixed
activity, 56
Solubility
Temperature dependence, 190
SOLUTION, 38
SOLUTION_MASTER_SPECIES, 34
SOLUTION_SPECUES, 35
sr(), 168
-start, 74, 161-steps, 162

T
Temperature dependency
Analytical expression for species, 35
Mineral kinetics, 164
van 't Hoff equation, 35
Time Series Plot, 153time(), 161
tot(), 169-total, 59
total_time(), 169, 175-totalsum, 197
totmole(), 193
Transition State Theory, 163

U
Units, 38
USE, 33, 135
USER_PUNCH, 73-end, 74
-integralsum, 192-mybasic, 74-plot_titles, 123-plot_units, 123-start, 74-totalsum, 197Vvan ‘Hoff equation, 35
Vanselow convention, 131

W
WATEQ activity correction equation, 36
Work Spaces
Graph workspace, 64
Input, 41
Output file workspace, 62

208
ISC: Public