Comput Geosci (2015) 19:655–671
DOI 10.1007/s10596-014-9458-3

ORIGINAL PAPER

Implementation and evaluation of permeability-porosity
and tortuosity-porosity relationships linked to mineraldissolution-precipitation
Mingliang Xie · K. Ulrich Mayer · Francis Claret ·Peter Alt-Epping · Diederik Jacques · Carl Steefel ·Christophe Chiaberge · Jiri Simunek

Received: 20 May 2014 / Accepted: 29 October 2014 / Published online: 20 November 2014
© Springer International Publishing Switzerland 2014

Abstract Changes of porosity, permeability, and tortuosity
due to physical and geochemical processes are of vital
importance for a variety of hydrogeological systems, including passive treatment facilities for contaminated groundwater, engineered barrier systems (EBS), and host rocks
for high-level nuclear waste (HLW) repositories. Due tothe nonlinear nature and chemical complexity of the problem, in most cases, it is impossible to verify reactive
transport codes analytically, and code intercomparisons are
the most suitable method to assess code capabilities andmodel performance. This paper summarizes model intercomparisons for six hypothetical scenarios with generally
Electronic supplementary material The online version of thisarticle (doi:10.1007/s10596-014-9458-3) contains supplementarymaterial, which is available to authorized users.
M. Xie () · K. U. Mayer
Department of Earth, Ocean and Atmospheric Sciences,
The University of British Columbia, Vancouver,
BC, Canadae-mail: mxie@eos.ubc.ca
F. Claret · C. Chiaberge
BRGM, 45060 Orleans Cedex, France
P. Alt-Epping
Rock-Water Interaction Group, Institute of Geological Sciences,
University of Bern, 3012 Bern, Switzerland
D. Jacques
Belgian Nuclear Research Centre SCK. CEN,
B-2400 Mol, Belgium
C. Steefel
Lawrence Berkeley National Laboratory, Berkeley,
CA 94720, USA
J. Simunek
University of California, Riverside, CA 92521, USA

increasing geochemical or physical complexity using thereactive transport codes CrunchFlow, HP1, MIN3P, PFlotran, and TOUGHREACT. Benchmark problems include
the enhancement of porosity and permeability through mineral dissolution, as well as near complete clogging due
to localized mineral precipitation, leading to reduction of
permeability and tortuosity. Processes considered in thebenchmark simulations are advective-dispersive transport insaturated media, kinetically controlled mineral dissolutionprecipitation, and aqueous complexation. Porosity changes
are induced by mineral dissolution-precipitation reactions,and the Carman-Kozeny relationship is used to describechanges in permeability as a function of porosity. Archie’s
law is used to update the tortuosity and the pore diffusion
coefficient as a function of porosity. Results demonstrate
that, generally, good agreement is reached amongst thecomputer models despite significant differences in model
formulations. Some differences are observed, in particular
for the more complex scenarios involving clogging; however, these differences do not affect the interpretation ofsystem behavior and evolution.
Keywords Permeability-porosity relationship ·Tortuosity-porosity relationship · Mineraldissolution-precipitation · Benchmark · CrunchFlow ·HP1 · MIN3P · PFlotran · TOUGHREACT
1 Introduction
The simulation of permeability and tortuosity evolution due
to porosity changes is of vital importance for the description of various hydrogeological processes, such as fluid
circulation in geothermal systems [1], acid rock drainage(ARD) attenuation and treatment [2], well fouling due tooxide precipitation [3], aquifer storage and recovery (ASR)

656

Comput Geosci (2015) 19:655–671[4], and the long-term geochemical evolution of host rockconsidered for high-level nuclear waste (HLW) repositories[5–7]. Porosity changes can be the result of many processes including mechanical deformation, swelling of clay
minerals (e.g., montmorillonite), dissolution/precipitation,and thermal deformation. Mineral dissolution and precipitation reactions can also have a significant effect on porous
medium porosity and in turn permeability and tortuosity.
Porosity changes in natural porous media occur normallyvery slowly, but in the long term can be of pivotal importance for the migration of fluids and solutes. If porosity
increases are substantial, preferential fluid migration pathways may develop, accelerating solute transport. On theother hand, for significant porosity decreases, nearly complete clogging may occur, substantially inhibiting fluid andsolute migration. Clogging in particular is associated withnumerical challenges due to the development of very small
water-rock ratios, or the complete disappearance of theaqueous phase. It is therefore useful to evaluate the capabilities and performance of reactive transport codes to simulate
evolving porosity and permeability, both for conditions ofporosity enhancement and reduction. It is also of interest toevaluate the effect of changing tortuosity on pore diffusion,in particular for clogging-dominated problems. Ideally, theverification of numerical models is conducted by comparing results to analytical solutions. For simulations involving
porosity evolution, this is only possible for problems with alimited set of reactions and subject to simplifying assumptions [8–10]. However, considering the complex natureof realistic reactive transport problems, model verification
by means of code intercomparisons through well-definedbenchmarks is in many cases the only practical method [11].

This paper presents and analyzes the results of an intercomparison involving five reactive transport codes (i.e.,
CrunchFlow, HP1, MIN3P, PFlotran, and TOUGHREACT) used for the investigation of six hypotheticalscenarios to evaluate formulations and implementations
for permeability-porosity and tortuosity-porosity relationships. Transport and reaction processes considered include
multicomponent aqueous complexation, kinetic surfacecontrolled reversible mineral dissolution and precipitation,
advective and diffusive mass transport, and coupling of
dissolution-precipitation-induced porosity change and fluid
flow processes. The Carman-Kozeny relationship was used
to describe changes in permeability, and the tortuosity andpore diffusion coefficients were updated based on Archie’slaw, both as a function of porosity.

2 Problem overview
Six benchmark levels with increasing complexity were
investigated in this exercise (Table 1)—the first level (B1)
considers mineral dissolution and associated permeability
enhancement, the second level (B2) focuses on a
precipitation-dominated (clogging) problem, and the thirdlevel (B3) extends the clogging problem of benchmark B2
by considering a significantly more complex geochemicalsystem (e.g., inclusion of additional chemical components
and minerals). Benchmarks B1–B3 have in common thatadvection is the only transport process considered. The
fourth level (B4) is based on benchmark B3, but considers
diffusive transport instead of advective transport to compare

Table 1 Overview of the benchmarks
Benchmark level

Description

B1
B2

Dissolution only
Clogging in a simple
geochemical system

B3

Dimension

1D homogeneous(2 m in length)

B4

B5

B6

# -Final simulation time

Clogging in a complex
geochemical system
2D heterogeneous(3 m in lengthand 2 m in height)

Processes

tFinal# [years]

Advection, kinetic dissolution
Advection, kinetic
dissolution/precipitation
Advection, kinetic
dissolution/precipitation,
redox reactions
Diffusion, no advection, kinetic
dissolution/precipitation,
redox reactions
Advection, diffusion, kinetic
dissolution/precipitation,
redox reactions
Advection, diffusion, kinetic
dissolution/ precipitation,
redox reactions

500
1000
300

3000

300

300

Comput Geosci (2015) 19:655–671

the codes under purely diffusive mass transfer conditions.
The fifth level (B5) combines benchmarks B3 and B4 andincludes both diffusive and advective transport. Benchmarks
B1–B5 were carried out for 1D flow and transport systems
with initially homogeneous medium properties. Benchmark
B6 further extends B5 into a 2D heterogeneous problem
with a randomly distributed hydraulic conductivity field.
Benchmark B1 considers the infiltration of a sulfuric acidsolution into a rock matrix containing calcite in an otherwise
nonreactive (inert) rock matrix. The only mineral phase considered is calcite; four components are included (H+ , Ca2+ ,2−
CO2−
3 , and SO4 ), leading to the formulation of 10 aqueous
complexes. Benchmark B2 also considers the infiltration ofa sulfuric acid solution into a calcite-containing rock matrix.
In addition to the primary mineral phase calcite, gypsum isconsidered as a secondary phase. The inflow solution contains a much higher SO2−
4 concentration. Benchmark B3
is based on B2 by adding Fe(II), Fe(III), K, Al, Na, andO2 (aq), thus increasing the number of aqueous components(10), secondary aqueous species (40) and mineral phases(6). Benchmarks B4 to B6 use the same reaction network
as benchmark B3. For benchmarks B4 to B6, the chemicalcomposition of the solution at the inflow boundary and thepore water in the domain, and the mineral composition of
the domain also remain the same as those in benchmark B3.
All simulations were conducted under standard pressure
and temperature (25 °C) and fully saturated conditions. The
geochemical database was derived from the database of
MINTEQA2 ([12], see Supporting Information).
This set of benchmarks was solved by five different reactive transport codes, namely CrunchFlow [13], HP1 [14, 15](restricted to B1–B5, because it is a 1D transport code),
MIN3P [16, 17], PFlotran [18], and TOUGHREACT [19–21]. An overview on the formulations and capabilities of theindividual codes is given in [13].

3 Theoretical backgroundThe generalized governing equations for multicomponent
reactive transport are provided in [13] and are therefore notrepeated here. This section focuses on the formulations used
to describe permeability-porosity and tortuosity-porosity
relationships, which are central to this benchmark. In
addition, the formulations for mineral dissolutionprecipitation reactions and activity corrections are provided,
considering that these relationships strongly affect theresults of the benchmarks.
3.1 Porosity, permeability, and tortuosity
The porosity change in each grid cell is calculated by summation over the volume fraction (ϕ) changes of the minerals

657

that occurred between the new time step (t + t) and that at
the previous time step (t). Assuming that porosity changes
are occurring slowly, porosity (φ) in each cell is explicitly
updated after completion of each time step:
φ t+t = φ t −Nm(ϕit+t − ϕit )

(1)

i=1

where Nm defines the number of the minerals. The hydraulicconductivity is also updated at the end of each time step
according to the Carman-Kozeny relationship [22]:t+t )3(φ(1 − φ t )2
t+t
Kt
K
= 
(2)
2
(φ t )3
1 − (φ t+t )Archie’s law is used to calculate the tortuosity τ as afunction of porosity at the new time level:ατ t+t = φ t+t
(3)
with the exponent α = 1/3, consistent with the formulation
by Millington and Quirk [23]. The pore diffusion coefficientat the new time level (Dpt+t ) is calculated as
Dpt+t = τ t+t D0(4)

where D0 is the free phase diffusion coefficient in water(assumed to be identical for all species). The effective
diffusion coefficient De is calculated as
Det+t = φ t+t τ t+t D0 = φ t+t Dpt+t(5)

3.2 Mineral dissolution-precipitation reactions
Mineral dissolution-precipitation reactions are describedas kinetic surface-controlled reactions based on the rate
expression:


IAPm
i
m
m
Ri = −ki 1 −(6)
Kim
where Rim is the reaction rate, IAPm
i is the ion activity product, Kim is the equilibrium constant, and kim is the effective
rate constant. The rate constants kim of primary minerals
vary as a function of mineral abundance at each time step.
For primary mineral phases, a two-third power relationshipof the form
kim,t = −kim,0

ϕit
ϕi0

2/3

(7)

is used to update the effective rate constant [24]. In this relationship, kim,t and ϕit define the effective rate constant andmineral volume fraction at time t, respectively; while kim,0and ϕi0 define the initial rate constant and mineral volume
fraction, respectively. The rate coefficients for secondary
minerals are assumed to remain constant.

658

Comput Geosci (2015) 19:655–671

Fig. 1 2D geometry (left) and initial hydraulic conductivity distribution (in [m s−1 ]) (right)

3.3 Activity corrections
Activity coefficients γid for all charged dissolved species
Adi , where Adi can be either a component as species in solution or an aqueous complex, are calculated based on themodified Debye-Hückel equation if ion-specific parameters
are available [12, 25]:√−Ad Zi2 I
dlog γi =(8)√ + bi I
1 + Bd ai I
where I is the ionic strength, Zi is the charge of the i th ion,
Ad and Bd are constants, ai and bi are ion-specific parameters. If ai is available, but not bi , Eq. 8 is used with bi = 0.
When ai is not provided, the Davies equation is used as an
approximation [12]:
 √

I
d
2
(9)
log γi = −Ad Zi√ − 0.24I
1+ I
Activity correction for neutral species excluding water iscalculated as described by Allison et al. [12]:log γid = 0.1I

(10)

The activity correction for water is defined by [12]:γH2 0 = 1 − 0.017

Nd


Cid

(11)

i=1

Table 2 Reaction
stoichiometries and solubility
constants for minerals

where Cid are the concentrations of the dissolved charged
species and Nd defines the number of the dissolved charged
species [12].
3.4 Numerical considerations
Benchmark levels B1–B3 do neither consider diffusion,
nor dispersion, in order to allow for a better evaluationof numerical dispersion and operator-splitting errors as afunction of different model formulations. Because significant spatial differences develop for porosity, hydraulicconductivity, and tortuosity, the spatial weighting of theseparameters is important. Harmonic averaging of all flow andtransport parameters (i.e., hydraulic conductivity and tortuosity) were adopted as the recommended default option
for the simulations. For spatial weighting of the advection terms, upstream weighting was used; implicit timeweighting was used for time integration.
Marty et al. [26] reported that spatial discretization intervals can have a substantial effect on the simulated clogging
position, which could impact the model intercomparison. To
avoid this issue, all codes used the same spatial discretization for each of the benchmarks B1–B6.
Similar to spatial discretization, it would be ideal touse the same temporal discretization amongst all codes formodel intercomparison. However, phase disappearance ofminerals and near complete clogging events often lead toconvergence problems, which can be effectively addressedby automatic adjustment of the time steps. Under theseMineralReaction

logK25

CalciteGypsumFerrihydrite
Jarosite
Gibbsite
Siderite

Ca2+ + CO2−3 ↔ CaCO3
2−
2+
Ca + SO4 + 2H2 O ↔ CaSO4 :2H2 O
Fe3+ + 3H2 O ↔ Fe(OH)3 + 3H++
K+ + 3Fe3+ + 2SO2−
4 + 6H2 O ↔ KFe3 (SO4 )2 (OH)6 + 6H
3++
Al + 3H2 O ↔ Al(OH)3 + 3H
Fe2+ + CO2−
3 ↔ FeCO38.47504.5800-4.89109.2100-8.1100
10.4500

Comput Geosci (2015) 19:655–671

659

Table 3 Physical and kinetic mineral parameters
Mineralkim,0[m2 mineral L−1 bulk]

Density[g cm−3 ]Mol. weight[g mol−1 ]Molar volume[cm3 mol−1 ]Update type

CalciteGypsumFerrihydrite
Jarosite
Gibbsite
Siderite

5 × 10−8
5 × 10−8
5 × 10−9
5 × 10−9
5 × 10−10
5 × 10−9

2.71
2.32
4.37
3.20
2.35
3.96

100.09172.17104.87494.81
78.00
115.86

36.93
74.21
23.99
154.6333.19
29.26

Two third
Constant
Constant
Constant
Two third
Two third

circumstances, differences in model formulation make it
difficult to use the same temporal discretization for allcodes. To evaluate the effect of adaptive time stepping or
constant time stepping for a range of time increments, aseries of simulations were performed with the MIN3P code.
It was found that the time increment size had an insignificant effect on the clogging position and outflow response(results not shown). It was therefore deemed adequate tolet the various participants choose time increments that are
most appropriate for the corresponding code and model
formulation on numerical consideration.
A performance and efficiency comparison was not a goalof the exercise. Information on time stepping and Newton
iteration count is therefore not provided.

4.2 Physical properties

4 Benchmark definition

For the 1D problems B1 to B5, the porous medium
is initially homogeneous with an initial porosityof 0.35 and an initial hydraulic conductivity of
1.16 × 10−4 m s−1 (10 m day−1 ). An average free phasediffusion coefficient of 1.0 × 10−9 m−2 s−1 for all aqueous
species is considered for benchmarks B4–B6. For problems
including advection, dispersion is neglected to minimize thecomplexity of the problem and facilitate a focused comparison of the different codes. In addition, dispersion becomesinsignificant when clogging occurs, due to reduced flow
velocities and the localized nature of pore clogging due tomineral precipitation. For the 2D problem B6, advective
and diffusive transport is considered in a similar fashion
to B5; however, the porous medium is characterized by aheterogeneous initial hydraulic conductivity field (Fig. 1).

4.1 Dimensions and spatial discretization

4.3 Boundary conditions for flow

The 2-m-long 1D domain for benchmarks B1–B5 was discretized using a discretization interval of x = 0.025 m,corresponding to 80 cells for uniform discretization, or 81
cells with half cells on the boundaries. The 2D domain(3 m in length and 2 m in height) for benchmark B6 wasdiscretized into 31 × 41 = 1271 cells (Fig. 1).

Benchmarks B1–B4 considered constant hydraulic heads at
the inflow (0.007 m) and outflow (0.0 m) boundaries. For
benchmark B4, no-flow conditions are enforced by setting
equal hydraulic heads at the inflow and outflow boundaries. For benchmark B6, the hydraulic head at the inflow
boundary was set at 0.0105 m to maintain the same
hydraulic gradient as with all other benchmarks (except
B4).

Table 4 Initial and boundary geochemical compositions in aqueous
phase (B1)
Primary
components

Unit

pH
Ca2+
CO2−
3
SO2−
4

–
mol l−1 H2 O
mol l−1 H2 O
mol l−1 H2 O

Initialcondition

Boundary
condition

9.38
1.57 × 10−4
2.57 × 10−4
1.00 × 10−10

3.0
1.00 × 10−4
1.0 × 10−2
6.46 × 10−4

4.4 Mineralogical and chemical parameters4.4.1 Aqueous components and speciation reactions
The components for benchmarks B1 and B2 are H+ , CO2−3 ,+ . An enhanced set of components
Ca2+ , SO2−
,
and
Na
4
2−
2+
is used for benchmarks B3–B6 (H+ , CO2−
3 , Ca , SO4 ,
3+
2+
3+
+
+
Fe , Fe , Al , K , O2 (aq), and Na ). Relevant aqueous complexation reactions and activity model parameters
are tabulated in the Supporting Information.

660
Table 5 Initial mineralcomposition and saturationindex (SI) in the initial andboundary solutions (B2)

Table 6 Initial and boundary
geochemical compositions inaqueous phase (B2)

Table 7 Initial mineralcomposition and saturationindex (SI) in the initial andboundary solutions (B3–B6)Table 8 Initial and boundary
geochemical compositions inaqueous phase (B3–B6)

Comput Geosci (2015) 19:655–671

Mineral

Initial volume fraction [m3 m−3 ]SI (IC solution)

SI (BC solution)

CalciteGypsum

0.3
0.0

0.0
−3.1

−9.1
−1.9

Primary components

Unit

Initial condition (IC)Boundary condition (BC)

pH
Ca2+
CO2−
3
SO2−
4
Na+

pH
mol l−1 H2 O
mol l−1 H2 O
mol l−1 H2 O
mol l−1 H2 O

9.33
1.70 × 10−4
2.70 × 10−4
1.70 × 10−4
3.20 × 10−4

3.0
1.0 × 10−4
1.0 × 10−2
0.2
3.96 × 10−1

Mineral

Initial volume fraction [m3 m−3 ]SI (IC solution)

SI (BC solution)

CalciteGypsumFerrihydrite
Jarosite
Gibbsite
Siderite

0.22
0.00
0.00
0.00
0.05
0.05

0.00
−2.78
0.00
−12.97
0.00
0.00

−8.93
−2.02
−0.07
6.66
−2.94
−10.89

Primary components

Unit

Initial condition (IC)Boundary condition (BC)

pH
pe
Ca2+
CO2−
3
SO2−
4
Na+
Al3+
K+
Fe2+
Fe3+

–
–
mol l−1 H2 O
mol l−1 H2 O
mol l−1 H2 O
mol l−1 H2 O
mol l−1 H2 O
mol l−1 H2 O
mol l−1 H2 O
mol l−1 H2 O

8.01
−0.67
4.71 × 10−4
2.19 × 10−3
1.70 × 10−4
1.54 × 10−3
2.81 × 10−7
1.00 × 10−5
6.59 × 10−6
2.53 × 10−8

3.0
17.53
1.00 × 10−4
1.00 × 10−2
1.00 × 10−1
9.09 × 10−2
1.43 × 10−2
7.67 × 10−5
1.14 × 10−8
2.23 × 10−2Comput Geosci (2015) 19:655–671

661

Fig. 2 Profiles of porosity (a), volume fraction of calcite (b), and hydraulic head (c) at 10, 100, and 120 years, and flux time curve (d) from theoutflow boundary simulated by CrunchFlow, HP1, MIN3P, PFlotran, and TOUGHREACT

4.4.2 Minerals

4.4.3 Initial and boundary conditions for reactive transport

Benchmark B2 The mineral phases considered in the system are calcite and gypsum. The domain initially includes
only calcite with a volume fraction of 0.30 m3 m−3(Table 5). Resident pore water is in equilibrium with calcite
and undersaturated with respect to gypsum (Table 6).
With the infiltration of a sulfuric acid solution (pH = 3)2−2++
containing CO2−
3 , Ca , SO4 , and Na , calcite dissolves
and gypsum precipitates, followed by gypsum re-dissolution
near the inflow boundary. The initial and boundary geochemical compositions for the aqueous phases are listed inTable 6. Boundary conditions are set as in B1.

Benchmark B1 The only mineral phase considered is calcite
at a volume fraction of 0.30 m3 m−3 , the remaining fractionof the solid phase is assumed inert. The solution domain
is initially occupied by an alkaline solution in equilibriumwith calcite. The infiltrating water is acidic with a pH of
3 and is undersaturated with respect to calcite (Table 4). A
third type boundary condition (specified mass flux) is usedat the inflow boundary and a free exit boundary is used at
the outflow.

Benchmarks B3–B6 The domain initially contains calcitewith a volume fraction of 0.22 m3 m−3 , as well as gibbsite (0.05 m3 m−3 ) and siderite (0.05 m3 m−3 ) (Table 7).
The initial and boundary total component concentrationsin the aqueous phases are listed in Table 8. The incoming
solution is acidic and undersaturated with respect to all minerals, except jarosite. This boundary condition was specified
to induce jarosite precipitation near the inflow. The background solution is slightly alkaline and is in equilibrium

Mineral phases and their parameters are listed in Tables 2
and 3. Benchmarks B1 and B2 use a subset of these mineralphases. All minerals are treated as kinetic surface-controlled
dissolution/precipitation reactions. The only exception isjarosite, which was treated as an equilibrium reaction by
HP1 owing to numerical difficulties. The initial rate constants (kim,0 ) are shown in Table 3.

662

Comput Geosci (2015) 19:655–671

Fig. 3 Profiles (B2) of porosity (a), volume fraction of gypsum (c) and calcite (e), hydraulic head (b), hydraulic conductivity (d) at 10, 100, and120 years, and flux from the outflow boundary (f) simulated by CrunchFlow, HP1, MIN3P, PFlotran, and TOUGHREACT

with calcite, siderite, gibbsite, and ferrihydrite. The initial
volume fraction of calcite (0.22 m3 m−3 ) is lower than thatin B2. The lower calcite content limits secondary gypsum
formation, which is offset by additional precipitation of ferrihydrite, siderite, gibbsite, and jarosite. This approach istaken to provide a suitable benchmark that evaluates system evolution under near clogging conditions. Boundary
conditions are set as in B1 and B2, with the exception ofbenchmark B4 which uses first-type (Dirichlet) boundary
conditions at the inflow and outflow ends of the column.

Comput Geosci (2015) 19:655–671

663

Fig. 4 Comparison of profiles (B3) of porosity, hydraulic head, and hydraulic conductivity at 10, 100 and 300 years, and outflux time curve fromthe outflow boundary simulated by CrunchFlow, HP1, MIN3P, PFlotran, and TOUGHREACT

5 Results
Simulated results by all codes are discussed in the following subsections. Reference model results produced by
MIN3P are also provided in tabular format in the Supporting
Information.
5.1 B1—porosity and permeability change due
to dissolution
Simulation results for benchmark B1 show the gradualdissolution of calcite along an initially narrow dissolution front, which widens as permeability and flow velocities increase (Fig. 2). Although porosity and permeability
increase substantially near the inflow end already at earlytime, total flow rate increases are relatively moderate during
the first 100 years. However, a rapid increase in flow rate
is predicted, once calcite depletion approaches the outflow
end of the simulation domain. All profiles of porosity, permeability, and hydraulic head at 10 and 100 years calculatedby the five codes are nearly identical, demonstrating very
good agreement between the models, despite different transport schemes and coupling methods. Slight differences exist
for results at 120 years. This time correlates with the
breakthrough of the dissolution front and is naturally mostsensitive to deviations. The time curves of the flux at theoutflow boundary are identical for all codes up to 120 years.
5.2 B2—clogging due to dissolution/precipitation(simple reaction network)
Simulation results show that the intrusion of the sulfuric
acid solution causes the dissolution of calcite and formation of gypsum (Fig. 3). Although gypsum re-dissolves near
the inflow, a narrow region of gypsum accumulation causes
clogging of the flow path. The zone of gypsum accumulation is already visible after 10 years; however, significant
porosity remains at this time. After 100 years, clogging isnearly complete and flow rates are greatly reduced, resulting
in negligible clogging front migration over the remaining

664

Comput Geosci (2015) 19:655–671

Fig. 5 Comparison of volume fraction profiles (B3) of calcite (a), gypsum (b), ferrihydrite (c), gibbsite (d), siderite (e), and jarosite (f) at 10, 100,and 300 years simulated by CrunchFlow, HP1, MIN3P, PFlotran, and TOUGHREACT

900 years. Simulated results show good agreement between
all codes. The maximum volume fraction of gypsum reaches
62.6 % and results in the reduction of porosity with aminimum value of 0.0016 after 1000 years at x = 0.425 m(MIN3P, PFlotran, and TOUGHREACT), at x = 0.40 mwith a minimum porosity of 0.003 (HP1) and at x = 0.40 m

Comput Geosci (2015) 19:655–671

665

Fig. 6 Comparison of profiles (B4) of porosity (left) and effective diffusion coefficient (De ) (right) at 100, 1000 and 3000 years simulated by
CrunchFlow, MIN3P, PFlotran, TOUGHREACT, and HP1 (at 100 and 1000 years)

with a minimum porosity of 0.0018 (CrunchFlow). The simulated profiles of the hydraulic head show a sharp decline
at the clogging point. The calculated flux at the outflow
boundary decreased rapidly within the first 100 years andmore gradually after 500 years (Fig. 3), showing good
agreement amongst all codes with small deviations by HP1
owing to the higher minimum value of porosity by HP1.
5.3 B3—clogging due to dissolution/precipitation(complex reaction network—advection)
Simulated results (Figs. 4 and 5) show that with the infiltration of the acidic solution, calcite dissolves, while otherminerals like gypsum, siderite, gibbsite, and ferrihydrite dissolve or precipitate depending on the time and location, thus
forming various porosity modifications along the flow pathat different time levels (Fig. 5). Jarosite precipitates near
the inflow; consequently, the porosity decreases near the
infiltration boundary. However, the major porosity reduction occurs due to the dissolution of calcite and precipitationof gypsum, similar to benchmark level B2. At 10 years, theminimum porosity reduces from the initial value of 0.35 to
0.05 at 0.25 m, showing a tendency for clogging (Fig. 4).
The “clogging point” develops and moves further downstream to 0.45 m at 100 years. Beyond 100 years, however,the clogging point moves more slowly (Fig. 4, top left),as clogging is nearly complete and water and solute fluxes
are greatly reduced. At 300 years, the minimum porosity reached is 0.0044 (CrunchFlow), 0.0020 (HP1), 0.0028(MIN3P), 0.0022 (PFlotran), and 0.0048 (TOUGHREACT).
Overall, the model results compare well for all fivecodes; all models are in good agreement with respect tothe prediction of the mineral distribution and clogging location. However, there appears to be significant differencesin terms of outflux at the end of the simulation (ranging
over approximately 2 orders of magnitude). On the other

hand, this apparent mismatch has to be put into perspective. Comparing the residual outflux to the initial outflux at
T = 0 demonstrates that all models predict a reduction offlux by greater than 99.9 %, which again points to a very
good agreement in terms of predicting both hydrogeologic
and geochemical evolution. Reasons for residual differences
are likely due to the high sensitivity of fluxes to minor variations in porosity, which in turn is affected by the formationand dissolution of multiple mineral phases. These results
highlight the challenges that numerical models encounterdue to clogging and phase disappearance.
5.4 B4—clogging due to dissolution/precipitation(complex reactions network—diffusion)
Diffusive transport results in a similar geochemical evolution as simulated in the advective case (B3), but withmineral dissolution/precipitation fronts that are moving
more slowly. Consequently, the kinetic limitations have asmaller effect on the results, which can be seen by thecomparison of the profiles for jarosite. Jarosite was treatedas a kinetically controlled mineral dissolution/precipitation
reaction by all codes except HP1. The simulated results
at 100 and 1000 years by HP1 agree better to those by
CrunchFlow, MIN3P and PFlotran in comparison to benchmark B3 (compare Figs. 7 and 5). Figure 6 (left) shows aregion of reduced porosity, leading to reductions in tortuosity and effective diffusion coefficients (De ) (Fig. 6 right).
The point with minimal porosity is located at 0.325 m after
3000 years. Simulated geochemical and hydrogeologicalevolution (Fig. 6) as well as mineralogical composition(Fig. 7) by all participating codes showed good agreement.
Generally, the simulated profiles of mineral volume fractions by MIN3P CrunchFlow and PFlotran agree very wellwith slight difference for ferrihydrite and jarosite. Simulated mineral composition profiles by TOUGHREACT and666

Comput Geosci (2015) 19:655–671

Fig. 7 Comparison of mineralvolume fraction profiles (B4) of
calcite (a), gypsum (b),ferrihydrite (c), gibbsite (d),siderite (e) and jarosite (f) at
100, 1000 and 3000 years
simulated by CrunchFlow,
MIN3P, PFlotran and
TOUGHREACT, and HP1 (for
100 and 1000 years)

HP1 showed slightly faster moving fronts for all minerals,
but visual agreement amongst all codes can be considered
as good.
5.5 B5—clogging due to dissolution/precipitation (complex
chemical reactions—advection and diffusion)
Overall, the predicted geochemical and hydrogeologicalevolution (Fig. 8), as well as mineralogical composition(Fig. 9) are in agreement for all participating codes. The
clogging positions of all simulations are at 0.5 m, whichclosely coincides with the position determined in benchmark B3. This indicates that the clogging position is controlled by the advective mass transport for the benchmark
set considered here. The distributions of most minerals in
B5 resemble those of B3 with slight differences in thelocation of the moving fronts except jarosite simulatedby HP1. This can be explained because jarosite precipitation is treated as an equilibrium phase instead of kinetic
precipitation, as specified in the benchmark
In contrast to benchmark B3, larger differences canbe seen for outflux and minimum porosity (Fig. 8). The

Comput Geosci (2015) 19:655–671

667

Fig. 8 Comparison of profiles(B5) of porosity (a), hydraulichead (b) and hydraulicconductivity (c) at 10, 100 and
300 years, and outflux time
curve from the outflow
boundary (d) simulated by
CrunchFlow, HP1, MIN3P,
PFlotran, and TOUGHREACT

simulated fluxes are fairly close within the first 90 years. At
later time, simulated fluxes by MIN3P, HP1, PFlotran and
CrunchFlow decrease gradually (Fig. 8), while TOUGHREACT results suggest a substantial reduction of permeability
that remains constant after about 160 years. Although all
fluxes are small in relation to the initial flux through the system, the final predicted outfluxes vary by more than 4 ordersof magnitude. The predicted outfluxes by MIN3P, PFlotran
and CrunchFlow vary by less than 2 orders of magnitude.
These differences suggest that the numerical formulation ofdiffusion processes plays a key role for simulating porosity evolution, despite the fact that the problem is initiallyadvection-dominated. In comparison to B3, reduction offluxes occurs smoothly (compare Figs. 4 and 8). This indicates that diffusion tends to slow down clogging in this
scenario but leads to substantial reduction of the minimum
porosity and flux in the long term.
5.6 B6—clogging due to dissolution/precipitationin heterogeneous media
Simulated porosities for the 2D heterogeneous scenario are
depicted in Fig. 10 for two selected observation points P1(x = 0.1 m, z = 1.0 m) and P2 (x = 0.5 m, z = 1.0 m)
for MIN3P, CrunchFlow, TOUGHREACT and PFlotran.
At P2 all simulated results of porosity evolution show

good agreement—with a rapid decrease during the first
20 years and a slight increase thereafter (Fig. 10). At P1,
the simulated porosity show the same trend for the first
10 years, but the minimum porosity values differ substantially. At later time, the porosity increases again, gently
according to MIN3P and CrunchFlow results, but morestrongly according to the results obtained with PFlotran and
TOUGHREACT.
The simulated flux curves by MIN3P, CrunchFlow,TOUGHREACT and PFlotran show a gradually decreasing
trend with good agreement (Fig. 11). The flux declined very
rapidly over the first 10 years, during which all simulated
results are almost identical. At later times, the flux continues to decrease but at a reduced rate. These results seemto contradict the simulated porosity evolution at the observation points, which showed substantial differences at P1.
However, these results also indicate that the overall permeability evolution is reproduced in a similar fashion by allcodes, and that differences are local.
This observation is confirmed further by examining
simulated 2D porosity distributions at 300 years. Results
obtained with the different codes showed very similar pattern illustrating the effect of heterogeneity on clogging(Fig. 12). The low porosity zones (in blue), distributed
between x = 0.4 and 0.7 m, tend to form a continuous low
permeability belt, which prevents flow through the domain.

668

Comput Geosci (2015) 19:655–671

Fig. 9 Comparison of volume
fraction profiles (B5) of calcite(a), gypsum (b), ferrihydrite (c),gibbsite (d), siderite (e) andjarosite (f) at 10, 100 and
300 years simulated by
CrunchFlow, HP1, MIN3P,
PFlotran, and TOUGHREACT

The high porosity zones (in red) formed through the dissolution of calcite and are most prevalent in the regions of higher
initial hydraulic conductivity (compare Figs. 12 and 1).
Closer inspection of simulated results shown in Fig. 12
reveal different patterns between the codes for the porosity
close to the inlet: porosity fingering near the inlet is similar
for MIN3P and CrunchFlow, while a layer of higher porosity
formed at the inlet boundary for PFlotran and TOUGHREACT. These differences are due to the slightly different
treatment of the boundary conditions and domain discretization by participating codes. MIN3P and CrunchFlow use

half-sized cells at the boundary and the nodes are locatedat the boundary. PFlotran uses fullsized cells at the solution
boundary, which results in a slightly different redistributionof flow near the boundary. TOUGHREACT uses special
cells at the boundary—thin external layers but with verylarge volume to fix both hydraulic head and chemical conditions. Similar to the PFlotran results, this approach leads
to a slightly different redistribution of flow near the boundary cells. Nevertheless, such influences are limited to thecells near the boundaries and dissipate within the solution
domain.

Comput Geosci (2015) 19:655–671

Fig. 10 Comparison of porosity evolution (B6) at observation point
P1 (0.1 and 1.0 m) and P2 (0.5 and 1.0 m) simulated by CrunchFlow,
MIN3P, PFlotran and TOUGHREACT

669

Fig. 11 Comparison of flux evolution at outflow boundary (B6)
simulated by CrunchFlow, MIN3P, PFlotran and TOUGHREACT

6 Discussion
We acknowledge that the Carman-Kozeny relationship andArchie’s law are not the only relationships to describe thedependency of permeability and tortuosity on porosity. Multiple alternative relationships exist to describe permeability
and tortuosity evolution (e.g., [27]). The model comparison presented here showed good agreement between the fivecodes, but this does not imply that any of the codes—as
is—is applicable to simulate permeability or tortuosity evolution at a field site or in a laboratory experiment. It needs
to be evaluated whether the implemented relationships

Fig. 12 2D distribution ofporosity (B6) simulated by
MIN3P, CrunchFlow and
PFlotran and TOUGHREACT at
300 years, and the positions of
the observation points P1(x = 0.1 m, z = 1.0 m) and P2(x = 0.5 m, z = 1.0 m)

are suitable for the site or experimental conditions. However, most of the alternative relationships such as the
Brinkmann or Fair-Hatch model for permeability will also
be a function of porosity [28] and the general findings ofthe present investigation should therefore also be relevant
for other approaches. Lastly, it has to be recognized thatpermeability and tortuosity are not only functions of porosity, but may also depend on micro-structural transformations
that are taking place, i.e., different secondary minerals may
affect permeability in different ways due to their different
growth habit [29].

670Comput Geosci (2015) 19:655–671

7 Conclusions

References

A sixlevel benchmark problem with increasing complexity for flow, geochemical and mineralogical conditions wasanalyzed with five different reactive transport codes (i.e.
CrunchFlow, HP1, MIN3P, PFlotran and TOUGHREACT)
The purpose of the benchmarking exercise was to evaluate the robustness of various numerical implementationsof permeability-porosity and tortuosity-porosity relationships to assess the feedback between flow and reactive
transport with the potential for local phase disappearance,
and the tendency for complete clogging due to mineral
precipitation.
Generally speaking, the simulation results appear to berobust and not dependent on the employed code. Simulating porosity enhancement due to mineral dissolution proved
to be an easier task than simulating clogging. Increasing complexity due to the inclusion of additional transportmechanisms and a more complex mineralogical assemblage had an impact on the agreement of the model results
obtained by the various codes. Although all codes showed
good agreement of the predicted mineralogical assemblage
and the clogging locations, absolute values of mass fluxes
differed substantially. As discussed, these results have to
be seen in perspective, considering that the results of allcodes agree in the sense that flux reductions are near 100 %for the clogging scenarios. Although differences exist interms of the absolute flux values, the impact on interpretation of the results would be very limited or even negligible.
Nevertheless, these results highlight the difficulties to simulate problems with pore clogging and phase disappearance.
Interestingly, the 2D problem (B6) showed better agreement
amongst the codes than the equivalent 1D problem (B5),
despite the additional complexity in terms of dimensionality
and heterogeneity. Although local differences in porosity and mineral assemblage were significant, agreement
between overall outflow was good.
This benchmarking exercise also highlights that theimplementation of the permeability-porosity and tortuosityporosity relationships, the activity correction model, andthe numerical methods (e.g. spatial weighting schemes)
can significantly affect the benchmarking results. A consistent numerical treatment is crucial for the benchmarking ofreactive transport codes.

1. Wagner, R., Kühn, M., Meyn, V., Pape, H., Vath, U., Clauser, C.:
Numerical simulation of pore space clogging in geothermal reservoirs by precipitation of anhydrite. Int. J. Rock Mech. Min. Sci.
42, 1070–1081 (2005)
2. Annandale, J.G., Jovanovic, N.Z., Benadè, N., Tanner, P.D.: Modelling the long-term effect of irrigation with gypsiferous water onsoil and water resources. Agric. Ecosyst. Environ. 76, 109–119(1999)
3. Houben, G.J.: Corrigendum to iron oxide incrustations in wells—
part 1: Genesis, mineralogy and geochemistry. Appl. Geochem.
18(6), 927–939 (2003)
4. Rinck-Pfeiffer, S.M., Ragusa, S.R., Sztajnbok, P., Vandevelde, T.:Interrelationships between biological, chemical and physical processes as an analog to clogging in aquifer storage and recovery(ASR) wells. Water Res. 34(7), 2110–2118 (2000)
5. De Windt, L., Pellegrini, D., van der Lee, J.: Coupled modeling of cement/claystone interactions and radionuclide migration.
J. Contam. Hydrol. 68(3–4), 165–182 (2004)
6. Xie, M., Kolditz, O., Moog, H.C.: A geochemical transportmodel for thermo-hydro-chemical (THC) coupled processes
with saline water. Water Resour. Res. 47(2), W02545 (2011).
doi:10.1029/2010WR009270
7. Liu, S., Jacques, D., Govaerts, J., Wang, L.: Conceptual model
analysis of interaction at a concrete–Boom Clay interface. Phys.
Chem. Earth (2014). doi:10.1016/j.pce.2013.11.009
8. Lagneau, V., van der Lee, J.: Operator-splitting-based reactive
transport models in strong feedback of porosity change: the contribution of analytical solutions for accuracy validation and estimator
improvement. J. Contam. Hydrol. 112, 118–129 (2010)
9. Hayek, M., Kosakowski, G., Churakov, S.: Exact analyticalsolutions for a diffusion problem coupled with a precipitationdissolution reaction and feedback of porosity change. Water
Resour. Res. (2011). doi:10.1029/2010WR010321
10. Hayek, M., Kosakowski, G., Jakob, A., Churakov, S.: A classof analytical solutions for multidimensional multispecies diffusive transport coupled with precipitation-dissolution reactions and porosity changes. Water Resour. Res. (2012).
doi:10.1029/2011WR011663
11. Carrayrou, J., Hoffmann, J., Knabner, P., Kräutle, S., de
Dieuleveult, C., Erhel, J., van der Lee, J., Lagneau, V., Mayer,
K.U., Macquarrie, K.T.B.: Comparison of numerical methods for
simulating strongly nonlinear and heterogeneous reactive transport problems-the MoMaS benchmark case. Comput. Geosci.
14(3), 483–502 (2010)
12. Allison,
J.D.,
Brown,
D.S.,
Novo-Gradac,
K.J.:
MINTEQA2/PRODEFA2, a geochemical assessment modelfor environmental systems: version 3.0 user’s manual. UnitedStates Environmental Protection Agency, Office of Researchand Development, Washington, DC, EPA/600/3-91/021, pp. 106(1991)
13. 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., Šimu̇nek, J., Spycher,N., Yabusaki, S.B., Yeh, G.T.: Reactive transport codes for
subsurface environmental simulation. Comput. Geosci. (2014).
doi:10.1007/s10596-014-9443-x
14. Jacques, D., Šimu̇nek, J., Mallants, D., van Genuchten, M.T.h.:
Modeling Coupled Hydrologic and Chemical Processes: LongTerm Uranium Transport following Phosphorus Fertilization.
Vadose Zone J. 7(2), 698–711 (2008). doi:10.2136/vzj2007.0084

Acknowledgments Financial support for this work was provided by
the Nuclear Waste Management Organization (NWMO) and the Natural Sciences and Engineering Research Council of Canada (NSERC)
through research grants held by K.U. Mayer. This research was also
supported by an internal BRGM research project entitled PMME REP
CODSPE2013.

Comput Geosci (2015) 19:655–671
15. Jacques, D., Šimu̇nek, J., Mallants, D., van Genuchten, M.Th.:
Modelling coupled water flow, solute transport and geochemical
reactions affection heavy metal migration in a Podzol soil. Geoderma 145, 449–461 (2008)
16. Mayer, K.U., Frind, E.O., Blowes, D.W.: Multicomponent reactive transport modeling in variably saturated porous media using ageneralized formulation for kinetically controlled reactions. Water
Resour. Res. (2002). doi:10:1029/2001WR000862
17. Mayer, K.U., MacQuarrie, K.T.B.: Solution of the MoMaS
reactive transport benchmark with MIN3P—model formulation and simulation results. Comput. Geosci. (2010).
doi:10.1007/s10596-009-9158-6
18. Lichtner, P.C., Hammond, G.E., Lu, C., Karra, S., Bisht, G.,
Andre, B., Mills, R.T., Kumar, J.: PFLOTRAN User manual:
A Massively Parallel Reactive Flow and Transport Model for
Describing Surface and Subsurface Processes (2013)
19. Xu, T., Sonnenthal, E.L., Spycher, N., Pruess, K.: TOUGHREACT
user’s guide: a simulation program for non-isothermal multiphase reactive geochemical transport in variable saturated geologic media. In: Lawrence Berkeley National Laboratory Report
LBNL-55460 (Ed.) (2004)
20. Xu, T., Sonnenthal, E., Spycher, N., Pruess, K.:
TOUGHREACT—a simulation program for non-isothermal
multiphase reactive geochemical transport in variably saturated geologic media: applications to geothermal injectivity andCO2 geological sequestration. Comput. Geosci. 32, 145–165(2006)
21. Xu, T., Spycher, N., Sonnenthal, E., Zhang, G., Zheng, L., Pruess,K.: TOUGHREACT Version 2.0: a simulator for subsurface

671
reactive transport under non-isothermal multiphase flow conditions. Comput. Geosci. 37, 763–774 (2011)
22. Carrier, W.D.: Goodbye, Hazen; Hello Kozeny-Carman. J.
Geotech. Geoenviron. Eng. (2003). doi:10.1061/(ASCE)10900241(2003)129:11(1054)
23. Millington, R.J., Quirk, J.P.: Permeability of porous solids. Trans.
Faraday Soc. 57, 1200–1207 (1961)
24. Lichtner, P.C.: In: Lichtner, P.C., Steefel, C.I., Oelkers,
E.H. (eds.): Continuum formulation of multicomponentmultiphase reactive transport, vol. 34, pp. 1–81. MineralogicalSociety of America, Washington (1996). (Chapter 1)
25. Parkhurst, D.L., Appelo, C.A.J.: User’s guide to PHREEQC—a
computer program for speciation, batch-reaction, one-dimensionaltransport, and inverse geochemical calculations. In: Survey,
U.S.G. (ed.): Water-resources investigations, Denver (1999)
26. Marty, N.C.M., Tournassat, C., Burnol, A., Giffaut, E., Gaucher,
E.C.: Influence of reaction kinetics and mesh refinement on
the numerical modelling of concrete/clay interactions. J. Hydrol.
364(1–2), 58–72 (2009). doi:10.1016/j.jhydrol.2008.10.013
27. Saripalli, K.P., Meyer, P.D., Bacon, D.H., Freedman, V.L.:
Changes in hydrologic properties of aquifer media due to chemical reactions: a review. Crit. Rev. Environ. Sci. Technol. 31(2),311–349 (2001)
28. Le Gallo, Y., Bildstein, O., Brosse, E.: Coupled reaction-flow
modelling of diagenetic changes in reservoir permeability, porosity and mineral compositions. J. Hydrol. 209, 366–388 (1998)
29. Gouze, P., Coudrain-Ribstein, A.: Chemical reactions and porosity
changes during sedimentary diagenesis. Appl. Geochem. 17, 39–
47 (2002)