J. Hydrol. Hydromech., 66, 2018, 2, 211–226
DOI: 10.1515/johh-2017-0049

The HPx software for multicomponent reactive transport during
variably-saturated flow: Recent developments and applications
Diederik Jacques1*, Jiří Šimůnek2, Dirk Mallants3, Martinus Th. van Genuchten4, 51

Engineered and Geosystems Analysis Unit, Belgian Nuclear Research Centre, 2400 Mol, Belgium.
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92521, USA. E-mail: jsimunek@ucr.edu
3
CSIRO, Urrbrae, SA 5064, Australia. E-mail: dirk.mallants@csiro.au
4
Center for Environmental Studies, CEA, São Paulo State University, Rio Claro, SP 13506-900, Brazil.
E-mail: rvangenuchten@hotmail.com
5
Department of Nuclear Engineering, Federal University of Rio de Janeiro, Rio de Janeiro, RJ, 21945-970, Brazil.
*
Corresponding author. E-mail: diederik.jacques@sckcen.be2

Abstract: HPx is a multicomponent reactive transport model which uses HYDRUS as the flow and transport solver and
PHREEQC-3 as the biogeochemical solver. Some recent adaptations have significantly increased the flexibility of thesoftware for different environmental and engineering applications. This paper gives an overview of the most significant
changes of HPx, such as coupling transport properties to geochemical state variables, gas diffusion, and transport in two
and three dimensions. OpenMP allows for parallel computing using shared memory. Enhancements for scripting may
eventually simplify input definitions and create possibilities for defining templates for generic (sub)problems. We included a discussion of root solute uptake and colloid-affected solute transport to show that most or all of the comprehensive features of HYDRUS can be extended with geochemical information. Finally, an example is used to demonstrate
how HPx, and similar reactive transport models, can be helpful in implementing different factors relevant for soil organic
matter dynamics in soils. HPx offers a unique framework to couple spatial-temporal variations in water contents, temperatures, and water fluxes, with dissolved organic matter and CO2 transport, as well as bioturbation processes.
Keywords: HYDRUS; HPx; Numerical modeling; Reactive transport; Organic matter; Bioturbation.
INTRODUCTION
Mathematical, process-based simulation models are indispensable tools for elucidating the role of coupled processes insoil systems, including their functioning and providing ecosystem services (Vereecken et al., 2016). Coupled multicomponent
reactive transport models for the so-called critical zone
(Brantley et al., 2007) are essential tools to analyze nonlinear
interactions between hydrological, geochemical and (micro)biological processes and to link the critical zone to otherspheres of the earth system such as the hydrosphere, atmosphere, biosphere, and geosphere (Li et al., 2017; Steefel et al.,2005). During the last two decades, a number of tools havebeen developed to simulate such coupled processes in natural
subsurface and engineered systems (Steefel et al., 2015). Numerous applications illustrate the enormous possibilities ofthese codes (e.g., see references in Li et al. (2017), Steefel et al.
(2015)). They typically couple a water flow equation (such as
the Richards equation for the vadose zone) with transport equations for heat and solutes based on the advection-dispersion
equation, a thermodynamically based geochemical solver, and akinetic reaction solver, while allowing for the coupling of multiple parameters with a series of state variables (Figure 1).
Because of the generic nature of such codes, they provide a
flexible (in terms of model formulations) and extendible (interms of model components) problem-solving environment totackle important environmental and engineering problems. This
paper focuses on the HPx simulator as developed within the
HYDRUS computer software packages (Šimůnek et al., 2016).
In this introduction we first introduce briefly two aspects which
illustrate the need for such modeling tools.
Reactive transport models allow for the direct coupling oftime-varying geochemical state variables with transport proper-

ties. In addition to the coupling of microstructural changes due
to mineral dissolution and precipitation to porosity and hydraulic properties (Freedman et al. (2004), Wissmeier and Barry(2009), Xie et al. (2015)), fluid properties such as viscosity anddensity, which are linked to the hydraulic properties of porous
media, depend upon the composition of the aqueous solution(Laliberté (2007), Laliberté and Cooper (2004), Sharqawy et al.
(2010)). Biomass accumulation or biofilm production during
microbiological reactions can also induce changes in hydraulicproperties, which may potentially influence the coupled hydrological-physical-geochemical-biological processes operative inporous media and, consequently, such processes as flow paths,
the performance of biofilters, and the transport of colloids (e.g.,
Bozorg et al. (2015a), Bozorg et al. (2015b), Carles Brangarí etal. (2017), Or et al. (2007), Rockhold et al. (2002), Yarwood etal. (2006)). The processes of precipitation and dissolution ofminerals may be affected by microbiological processes, such as
in the engineering practice of microbially induced calcium
carbonate precipitation (Hommel et al., 2015; Martinez et al.,2014). Another process that can produce changes in porosity
and the hydraulic properties is the migration of fines (solidmineral particles), which is influenced by both hydrodynamic
and chemical factors (Bennacer et al., 2017; Mays and Hunt,2007; Yu et al., 2013). Again, this process is relevant in manyengineering applications, but also within the field of soil development. Finally, under saturated conditions, biogenic processes
such as methanogenesis and denitrification may result in theformation of gas bubbles, which may result in a reduction of thehydraulic conductivity due to blocking of pore spaces (Amosand Mayer, 2006).
Organic matter dynamics in soils is a crucial process affecting many environmentally-related issues. Degradation produces
inorganic carbon, which can be released into the atmosphere.

211

Diederik Jacques, Jiří Šimůnek, Dirk Mallants, Martinus Th. van Genuchten

Flow
Transport
A
E
Bh1
Bh2
Bh/C

depth
(cm)
0
7
19
24

Moisture
Temperature
Concentration
Water flux
Root uptake

Biogeochemistry

28
50

C1

Precipitation75
Evaporation
C2
Transpiration
Advection
Diffusion
Dispersion
Heat advection&conduction
Heterogeneity
Macro-Micro domain
Root growth
Colloidal transport…

Moisture
Temperature
Concentration
Porosity
Hydraulic, solutetransport andheat transport
propertiesHYDRUS1/2/3

Aquatic speciation& complexation
Exchange
Sorption
Solid-Solutions
Minerals
Mobile surface
Kinetics
Microbiologicalmediated reactions
High-strength…

PHREEQC-3

Fig. 1. Reactive transport models coupling water flow and solute transport with biogeochemical processes.

Consequently, CO2 sequestration in soils may be a factor incontrolling atmospheric CO2, and thus climate change (ToddBrown et al., 2013). And vice versa, climate change, together
with land use, influences soil organic matter (SOM) dynamics(e.g. Batlle-Aguilar et al. (2011)). Dissolved organic matter(DOM) plays also an important role in many soil processes
including contaminant transport (e.g., Maes et al. (2006),
Paradelo et al. (2012)). For example, Leterme and Jacques(2015) showed that DOM concentrations are among the mostsensitive factors in simulating Hg fate in soil systems.
As reported by Manzoni and Porporato (2009), a large number of models exist for simulating SOM dynamics. Many SOM
models are so-called compartment models, which divide theheterogeneous nature of SOM into a number of homogeneous
pools, and then simulate the fluxes between these pools. Someof these models have a rigid structure in the sense that the number of pools, the abiotic factors that must be considered, theirfunctional forms, and the extent of biotic-abiotic interactions,are fixed. It is therefore not always possible to perform a structured comparison between different types of model formulations. Reactive transport models, on the other hand, provide a
flexible structure that enables the introduction of different SOM
degradation networks within the same simulation framework. In
addition, existing compartment models sometimes strongly
simplify flow and transport processes and/or the geochemistry
involved, such as simplifying or ignoring pH and buffering
effects on degradation or CO2 geochemistry. For predictions
that are more accurate, also within the framework of climate
predictions, SOM models should consider temporarily andspatially variable abiotic factors such as temperature and moisture, nutrient dynamics, microbial activity, solid-water interac-

212

tions, and transport processes of solid and dissolved organic
matter and/or inorganic carbon (e.g., Schmidt et al. (2011),
Tang and Riley (2015)). Reactive transport models should
account for these factors and effects, such as illustrated in recent studies by Tang et al. (2013) and Riley et al. (2014).
The main objective of this paper is to give an overview of
new capabilities of the HPx codes compared to the model descriptions given previously by Jacques et al. (2006) and Jacqueset al. (2008a). We also wanted to illustrate some of these capabilities as used in several recent studies, as well as present an
example involving different SOM degradation networks andbioturbation processes.
THE HPX REACTIVE TRANSPORT CODE
The HPx codes couple the HYDRUS-1D and HYDRUS(2D/3D) software packages (Šimůnek et al., 2016) with
PHREEQC (Parkhurst and Appelo, 2013) using a sequential
non-iterative coupling approach as explained by Jacques et al.
(2006). HYDRUS-1D or HYDRUS (2D/3D) act as the solver
for the hydrological and physical processes, including variablesaturated water flow, solute transport, diffusion in the gas
phase, and heat transport, whereas PHREEQC is the solver forthermodynamic and kinetic geochemical processes. While thefirst versions of HPx used PHREEQC-2.4 (Jacques andŠimůnek, 2005; Jacques et al., 2006), the most recent versions
2.4 of HPx use PHREEQC 3.1. Two- and three dimensional
flow and transport is included as well (Šimůnek et al., 2012).
Coupling with PHREEQC is performed using so called“hard coupling”, which is in principle similar to the procedures
available in the more recently developed wrappers iPhreeqc

The HPx software for multicomponent reactive transport during variably-saturated flow: Recent developments and applications
Read phreeqc.in
file up to END

Data Block
TRANSPORTNo

Geochemical calculations
including
CALCULATE_VALUES(definition)

Yes
Read geometry, flow andtransport
Initial calculations including
CALCULATE_VALUES (initial)

Update transport
properties

Yes

Update w, h, s
T(0), C(0), G(0)No
Start new time step
H(n), R(n) T(n), C(n), G(n)Update boundary
conditions

Water Flow

Root water
uptakeRoot water
uptakeCalculate w using
T(n)Calculate H(n+1) using
R(n)

Yes

Calculate R(n+1) using
H(n+1)

No
Calculate h using
H(n+1)

No
Heat Transport

Heat transport

Yes

Update hproperties

Yes

Calculate h using H(n+1)
and G(n) by
CALCULATE_VALUES(Caph, Condh, Disph)No

Calculate s using
H(n+1)

No

n=n+1

Solute transport

Calculate T(n+1) using
H(n+1)

Update sproperties
Yes

Biogeochemical
Reactions

Calculate s using H(n+1),
T(n+1) and G(n) by
CALCULATE_VALUES(Tort
W, TortG, DispW)CALCULATE_VALUES (br)Update w
properties

Yes

Calculate C(n_T) using
H(n+1)

Calculate C(n+1) and G(n+1)
using C(n_T), H(n+1), T(n+1)

CALCULATED_VALUES(update)No

No

End ofsimulation
Yes
End

Fig. 2. Flow chart of HPx. The letters w, h, and s indicate parameters for water flow, heat transport, and solute transport, respectively. The
variables H, R, T, C, and G indicate state variables for water flow, root water uptake, heat transport, solute transport, and biogeochemical
processes, respectively. The boxes with blue background are steps run in parallel using OpenMP. CALCULATED_VALUE(type) indicatesthe execution of scripts labeled with type.

(Charlton and Parkhurst, 2011) and RMPhreeqc (Parkhurst andWissmeier, 2015). Due to the hard coupling, HPx does notalways include the latest version of PHREEQC. However, HPx
includes additional flexibilities compared to the original

PHREEQC, such as linking geochemical definitions (e.g., solution compositions, exchange sites, minerals, kinetic reactants)
to the material distribution, or exchanging information between
the flow model and the BASIC scripts (e.g., pressure heads,

213

Diederik Jacques, Jiří Šimůnek, Dirk Mallants, Martinus Th. van Genuchten

velocities, and scaling factors). The geochemical input is defined using the standard PHREEQC input conventions using
data blocks and identifiers, thereby allowing a direct implementation of the PHREEQC models into HPx. The scripting possibilities are now enhanced by introducing global variables ordata tables, which can be used as inline variables or as BASIC
scripting variables, using different types of scripting blocks toperform specific calculations, and introducing a specific version
of BASIC (such as MyBasic; https://github.com/paladint/my_basic) with structured grammar for some specific calculations). Templates for generic types of problems can be generated in advance and then imported into additional input files thatcan be created during project execution using the INCLUDE$
statement. Geochemical calculations are performed in parallelusing the OpenMP (www.openmp.org) shared memory approach and PHREEQC objects. The main object reads then the
input and performs all initial calculations. The information(thermodynamic data and geochemical models) is copied to auser-defined number of additional objects representing chunks in
the OpenMP formalism. The nodes are assigned evenly between
the chunks. A dynamic schedule controls the calculations of thenodes in a chunk on the available threads on a computer. Afterthe geochemical calculations, the information from different
chunks is assembled and sent to the transport calculations or tothe model output files. Because of the OpenMP implementationand the enhanced scripting possibilities, the current flow chart ofHPx (see Figure 2) is an extended version of the chart described in

Jacques et al. (2006).
Simulation output data now also contain variables, which areintegrated over a part of the domain, such as over a particularsoil material. HP2 makes it possible to obtain geochemical
information along mesh lines or cross-sections at different
times as defined in the standard HYDRUS (2D/3D) graphicaluser interface.
FLOW, TRANSPORT, AND BIOGEOCHEMICALPROCESSES
Conceptual models
Figure 3 shows schematically the different physical waterflow and solute transport models implemented in HPx. The
software considers flow, transport, and the prevailing biogeo-

Fig. 3. Conceptual physical and geochemical models implemented in HPx. Black arrows represent transport processes within a single
phase, and red arrows exchange processes between different phases. The top left, top right, and bottom plots indicate water flow, solutetransport, and biogeochemical reactions in a single-porosity continuum, in a single-porosity continuum with an immobile solute domain,and in a dual-porosity model, respectively (εs - the solid content [dm3 dm–3], θa - the air content [dm3 dm–3], θw - the water content [dm3 dm–3],
Γw - the water exchange between the immobile and mobile domains [dm3 dm–3 T–1]Γwa - the solute exchange between the aqueous phaseand the gas phase [mol dm–3T–1], Γws - the solute exchange between the aqueous phase and the solid phase [mol dm–3T–1], Γws,m - the soluteexchange between the mobile aqueous phase and the solid phase [mol dm–3T–1], Γws,im - the solute exchange between the immobile aqueous
phase and the solid phase [mol dm–3T–1], Γi - the solute exchange between the immobile and mobile domains [mol dm–3T–1], Sr,w - root
water uptake [dm3 dm–3T–1], Sr,s,i – root solute uptake [mol dm3 T–1], Sc,,i – kinetic transformations in the aqueous phase [mol dm–3 T–1]).

214

The HPx software for multicomponent reactive transport during variably-saturated flow: Recent developments and applications

chemical processes at the continuum scale with three main
phases being water (θw), air (θa) and the solid phase (εs) as
shown in the figure. The water and solid phases in the dual
porosity modeling options of the code are divided into twoparts, with one part comprising the mobile phase in whichwater, heat and solutes are allowed to move, and another partbeing the immobile domain where water is assumed to be stagnant.
Variably-saturated water flow, multicomponent solutetransport and heat transport are described using the Richards
equation, the advection-dispersion equations (for multiple
solutes), and (neglecting energy transported by vapor transport)
the convection-dispersion equation, respectively. Transport ofreactive gases, which is crucial for simulating biological processes or solid phase alterations such as carbonation of concrete(e.g., Phung et al., 2016) or weathering (Li et al., 2017), iseither only by diffusion or advection. The latter is based on theassumptions of negligible air compressibility due to the fastexchange with the atmosphere, zero air flux at the bottom of the
profile, and immediate matching of volume changes in the total
water content by volume changes in the gas content.
Possible sinks of water in the flow domain are uptake byplants roots (Sr,w) and exchange between the mobile and immobile water phases (Γw), the latter only if dual-porosity waterexchange is considered. Sinks or sources of solutes from theaqueous phase are via (i) exchange with roots (Sr,i with i referring to a solute component), (ii) kinetic transformations such asradioactive decay, degradation, or redox reactions (Sc,i), (iii)exchange between the mobile and immobile water phases (Γi),
and (iv) exchange with the gaseous (Γwa) or solid phases (Γws,mand Γws,im) via thermodynamic speciation reactions or kinetic
rate equations. The latter reactions represent processes such as
sorption, ion exchange, surface complexation, gas-aqueous
phase equilibria, and mineral dissolution and precipitation.
Mathematical models
Water flow in its most general form is described using a dual-porosity formulation of the Richards equation (Šimůnek andvan Genuchten, 2008):

∂θ m∂   A ∂h=
+ K izA   − Sr,w − Γ w
 K  K ij∂t∂xi  ∂x j


ration. Actual transpiration is calculated based on the spatialtemporal root distribution, while stress factors in the soil arelinked to local water contents and salinities. The codes may be
used to account for root water uptake compensation by allowing the roots to take up more water in less stressed soil zones(Šimůnek and Hopmans, 2009). Boundary conditions are either
system-independent (e.g., constant or variable pressure head
(Dirichlet), flux (Neumann), or gradient boundary conditions)or system-dependent by considering atmospheric conditions
involving time-dependent precipitation, potential evaporation,and potential transpiration rates.
Heat transport is described as:

Cs (θ )

∂T
∂ 
∂T 
∂T
=
 λij (θ ) − Cw qi∂t ∂xi ∂xi ∂xi

where Cs and Cw are volumetric heat capacities [ML–1T–2K–1] ofsoil and water, respectively, and λij is the apparent thermal
conductivity of the soil [MLT–3K–1]. The reactive module uses
temperature to correct equilibrium constants or change rate
coefficients as desired using user-defined relations. Equation(2) is valid for both single- and dual-porosity flow models whenassuming local temperature equilibrium, i.e., that temperature is
the same in the solid phase, as well as in the mobile and immobile liquid phases.
Solute concentrations in the aqueous phase are expressed interms of total concentrations for a given primary species Ck(mol/kgw), which includes the concentration of all aqueous
equilibrium species, cl (mol/kgw) or so-called secondary species:Na

Ck = ck + υlk cl

∂θim
= Γw
∂t
where θm is the volumetric water content in the mobile domain[L3L–3], θim is the volumetric water content in the immobile
domain [L3L–3], h is the pressure head [L], xi (i = 1, 2, 3) are thespatial coordinates [L], t is time [T], Sr,w is the root water uptake term [T–1], Γw is the water exchange term [T–1] driven by afirst-order rate equation proportional to the difference in watercontents or pressure heads between the mobile and immobile
domains, KAij are components of a dimensionless anisotropy
tensor KA, and K is the unsaturated hydraulic conductivity function [LT–1]. The single-porosity models assume that θm is equal
to the total water content θ, in which case θ im and Γm become
equal to zero. The second part of Eq. (1) for the immobile region then does not apply.
Several analytical forms are available to describe the unsaturated soil hydraulic properties θ(h) and K(h) including Brooks
and Corey (1964), van Genuchten (1980), Kosugi (1996), andDurner (1994). Root water uptake is driven by potential transpi-

(3)

l =1

where υlk is the stoichiometric coefficient of species k in secondary species l, k = 1,…, Nc (Nc is the number of primary
species) and l = 1,…, Na (Na is the number of secondary species). For each primary species, solute transport in the aqueous
phase is described using advection-dispersion-reaction equations of the form:∂θ m Cm, k
(1)

(2)

∂t

=

∂Cm, k  ∂qi Ck∂ w− Sr,w Cr,k θ m Dij−∂xi ∂xi  ∂xieqkin−Γ k − Γ wa, k − Rm,
k − Rm,k∂θim Cim, k

(4)

eq
kin
k − Rim,k − Rim,k=Γ∂tΓ k = ω ( Cm, k − Cim, k ) + Γ w ( H ( Γ w ) Cm, k + H (−Γ w )Cim, k )where Cm,k and Cim,k are the total concentrations of the k species
in the mobile and immobile domains [mol/dm3], respectively,θm and θim have the unit [dm3/dm3], Dwij is the dispersion coefficient tensor [L2T–1], qi is the ith component of the Darcian
fluid flux density [LT–1], Cr,k is the concentration of the sinkterm [mol/dm3], Γk represents solute exchange between theimmobile and mobile domains [mol/dm3T–1], Γwa,k accounts forthe exchange between the aqueous and the gas phase
eq
eq
kin
kin
, Rim,k, Rm,kand Rim,kare aqueous sink/source[mol/dm3T–1], Rm,k
3 –1
terms [mol/dm T ], with the superscripts eq and kin referring
to equilibrium and kinetic reactions, respectively, and subscripts m and im to the mobile and immobile domains, respectively, ω is the mass transfer coefficient [T–1] and H() is the
Heaviside step function.

215

Diederik Jacques, Jiří Šimůnek, Dirk Mallants, Martinus Th. van Genuchten

Reactive gases such as O2 and CO2 play important roles inmany soil processes (e.g., pyrite oxidation, Mayer et al. (2015),
and hence must form an integral part of reactive transport models. Diffusion in the gas phase is described as:∂θ a Ca, k∂t∂ a ∂Ca , k = θ a Dij
 + Γ wa ,k
∂xi 
∂xi 

(5)

where Ca,k is the concentration in the gas phase of component k[mol/dm3], Daij is the gaseous diffusion tensor [L2T–1], andother symbols as defined before.
The geochemical model calculates aqueous concentrationsof the primary and secondary species using an aqueous equilibrium chemical model written as:
Nc

Al  ν lk Ak

(6)

k =1

where Al and Ak are the chemical formulae of the species, νlk arestoichiometric coefficients and Nc is the number of primary
species. Mass action laws permit calculations of the aqueousconcentrations as follows:Nc∏ (γ c )

ν lk

k k

cl = k =1

Kl γ l(7)

where Kl is the equilibrium constant for the reaction, defined as
product of the activities of the reaction products divided by the
product of the activities of the reactants, and γi and γj are theactivity correction coefficients to obtain activity a [–] when
multiplied with c. The thermodynamic constants are adapted fortemperature using the enthalpy of the reaction in the van‘t Hoff
equation or by a six-term analytic equation as a function oftemperature (Parkhurst and Appelo (2013), see e.g. also Eq. (2)in the appendix). Some of the parameters of this equation canbe linked to thermodynamic variables, such as the change ofentropy, enthalpy and heat capacity of reaction, at least for
relatively small temperature variations (e.g., Puigdomènech etal. (1997), Jacques (2009)) that are expected for soil conditions.
Several activity correction models are implemented, such as
ion-association models (Davies equation, extended DebyeHûckel equation), Pitzer model, and a model based on specificion interaction theory (SIT). Ion exchange, surface complexation, mineral, gas-aqueous phase and solid solutions equilibrium models are available heterogeneous reactions to calculateeqeq
the Γwa,k, Rm,kand Rim,k, terms in Eq. (4).
Equilibrium between the gas and aqueous phases is given as:

Gl  Al(8)

where Gl is the chemical formula of the gaseous component.
The equilibrium constant of the reactions is the ratio of theactivity of Al over the fugacity of Gl. With the assumption of an
ideal gas phase, which is valid for low total gas phase pressures, the fugacity coefficient is 1 in which case the fugacity
equals the partial pressure of Gl. Fugacity correction factors canbe calculated with the Peng and Robinson equation of state(Peng and Robinson, 1976). However, for total gas phase pressures lower than 10 atm, the fugacity coefficient will be close to
1 (Appelo et al., 2014).
A general ion exchange reaction in which the site capacity isalways occupied by ions of the opposite charge (i.e., a charge

216

neutral reaction) is given as:

ν X X k + υl Al  ν X X l +ν k Ak
k

l

(9)

where X denotes the chemical formula of an exchange speciesand ν is the stoichiometric coefficient, which is linked to athermodynamic constant using different conventions for theactivity correction of the exchange species including GainesThomas, Gapon, and Vanselow convention. The RothmundKornfeld equations (Bloom and Mansell, 2001; Bond, 1995),
allows one to consider the effects of a changing selectivity withsite occupancy (e.g., Jacques et al. (2012)). Different exchangesites can be combined to define a multisite exchange model(e.g., Jacques et al. (2008b)).
A general equation for surface complexation is:

S zs ( Akzk )ν k +ν l Alzl  S zs ( Alzl )ν l +ν k Akzk

(10)

where S represents the chemical formula of a surface masterspecies, and z the charge of the surface master species, theaqueous species, or the complexed species. Chemical equilibrium is obtained when the ratio of the activity product is equal tothe equilibrium reaction constant, possibly corrected for electrostatic effects. These effects can be either neglected using a
non-electrostatic model or incorporated by using a diffuse
double layer model (Dzombak and Morel, 1990) or the CDMUSIC model (Hiemstra and VanRiemsdijk, 1996). The composition of the double layer to counterbalance the surface isestimated by explicitly calculating the charge distribution in thelayer (Borkovec and Westall, 1983) or by the Donnan model(Appelo and Wersin, 2007).
A general mineral dissolution/precipitation reaction is of the
form:
Nl

M l  ν kl Ak

(11)

k =1

where M represents the chemical formula of the mineral and Nl
the number of aqueous species in the reaction equation. Since
the activity of a pure phase (mineral) is assumed to be equal to1, the equilibrium constant is the activity product of the aqueous species. In ideal solid solutions, which are homogeneous
mixtures of different minerals, the activity of the phase is lessthan 1, leading to a more stable phase compared to the purephase. For non-ideal binary solid solutions, described using the
Guggenheim series expansion model (Guggenheim, 1937), a
mechanical mixture may sometimes be more stable than thehomogeneous mixture, resulting in an immiscibility gap orphase separation.
kin
kin
and Rim,kterms in Eq. (4) represents homogeneous
The Rm,kand heterogeneous kinetically controlled reactions, respectively. Homogeneous reaction networks may involve redox, degradation or decay reactions. Each species in such a degradation
network must be treated as primary species for which a separatetransport equation is required. Typical heterogeneous reactions
are dissolution and precipitation of minerals, for which different kinetic expressions exist. One group is based on transition
state theory (Aagaard and Helgeson, 1982), which allows forirreversible reactions; the rate parameter may depend on solution composition variables such as pH or activities of otherspecies, which could act as catalyzer or inhibitor. However,irreversible rate equations based on a shrinking core model(e.g., Liu et al. (2014)) or equations in which reaction rates arecontrolled by microbiological processes can be described as

The HPx software for multicomponent reactive transport during variably-saturated flow: Recent developments and applications

well. Monod terms for electron acceptors and donors can act as
catalytic and inhibitory terms in the microbially-mediated rateequations. Mayer et al. (2002) gave a general rate equation forboth surface and transport controlled kinetics for mineral dissolution and precipitation. The scripting possibilities in HPx and
PHREEQC permit users to incorporate these or alternative rateequations. In addition, rate equations can be defined for sorption processes when implemented as exchange or surface complexation reactions.
Flow-transport-geochemistry interactions
The concept of scaling hydraulic properties (Vogel et al.,1991), which is typically used in soil sciences to describe their
spatial variability, provides a flexible framework for describing
the temporal changes of these properties as well. Changes inporosity are then linked to a scaling factor of the water content,αθ [–], as:α θ (t , x) =θs (t , x) − θ r*
θs* − θ r*

(12)

where θs is the time- and space-variable saturated water content(or the porosity) [L3 L–3], and θs* [L3 L–3] and θr* [L3 L–3] are
the reference saturated and residual water contents, respectively, for a specific material. Variations in the hydraulic conductivity in time and space, K(x,t) [L T–1], are described by ascaling factor of the hydraulic conductivity, aK [–], as follows

K (t, x )α K (t , x) = s *
Ks

(13)

where K*s is the reference saturated hydraulic conductivity for a
specific material [L T–1]. This formulation enables one to include the effects of temperature and geochemistry on viscosityand density, as well as of precipitation/dissolution/clogging
processes on the microstructure. Also, changes in the apparent
air-entry value can be described using a scaling factor for thepressure head, αh [–], such that h/αh rather than h (for the reference hydraulic properties) is used to obtain the water content
and hydraulic conductivity. HPx does not provide any functional relation for calculating the scaling factors. Via the input file,users can implement their preferred model relating geochemicalvariables to the scaling factors. In the current implementation,change of porosity and scaling factor for the pressure head (inthe dual porosity model with water exchange driven by pressure
head gradient) are limited to the mobile domain. Assuming a
constant porosity of the immobile domain, the change in θs*
reflects a change in porosity of the mobile domain only. One
limitation is that the scaling factors only describe linear changes in the hydraulic properties. For example, they do not allow
for capturing a model for hydraulic property changes whichassumes that precipitation and dissolution reactions occur onlyin the fully water-saturated part of the pore space (e.g., themodel of Wissmeier and Barry (2009)). However, when thereactions take place in all of the pore space, the linear scaling
model is applicable (Wissmeier and Barry, 2010).
HPx can also update solute and heat transport properties
such as the tortuosity of the aqueous and gaseous phases, andthe dispersivity, heat capacity, heat conductivity, and heat dispersivity. Tortuosity is often linked to porosity using Archie’s
law or the model of Millington and Quirk (1961). For cement
systems, some models discriminate between the types of porosity (gel versus capillary porosity) and associated cement

phases (Patel et al., 2016). The heat capacity can be linked toindividual heat capacities of the minerals in the system. However, not much literature evidence exists for a model linking
dispersivity and conductivity to geochemical state variables.
Root solute uptake
Kinetic uptake and exudation processes can be implemented
in HPx using rate expressions defined in the input files. Examples of PHREEQC implementations of these types of processes
were given by Seuntjens et al. (2004) and Nowack et al. (2006).
Compared to the HYDRUS implementations, also available in
HPx, some additional features can be taken into account whenkin
in Eq. (4) are used, such as speciesthe rate equations Rm,kspecific uptake rather than uptake of an element irrespective ofthe distribution among different aqueous species, inclusion ofcounter-ions, root exudation, and uptake parameters depending
upon soil geochemistry. It is also possible to incorporate thekin. Following Silberbush et al.
term Sr,wCr,k in Eq. (4) in Rm,k(2005), a general formulation of root solute uptake that combines passive and active uptake is given byrl = pl + al
= pl + Aroot J l ,max (ωgc )pa,lcl − cl ,min

(14)

Kl ,m (ωgc ) + (cl − cl ,min )= min[ S w,r min cl , cmax,l  ,Aroot J p,max,l (ωgc ) ,pp,l ]= 0(15)

if cl > cr,lotherwise

where rl refers to the uptake of a species l, pl and al are passive
and active uptake rates [mol/dm3 T–1], cmax,l is the maximum
allowed concentration for passive uptake [mol/dm3], Aroot is theroot surface area [L2/dm3], Jp,max,l is the maximum allowed
passive nutrient uptake rate [mol L–2 T–1], ω is the current geochemical condition, pp,l is the potential passive nutrient uptake
rate [mol/dm3 T–1], cr,l is the critical concentration below which
passive root uptake is zero [mol/dm3], Jl,max is the maximum
active nutrient uptake which may depend on geochemical conditions or the maximum potential active uptake [T–1], Kl,m is the
Michaelis-Menten parameter [mol/dm3], and cl,min is the minimum nutrient concentration below which active uptake is zero.
The expression for passive uptake is generic as it limits uptakeby (i) a maximum allowable concentration for uptake, (ii) a
maximum nutrient uptake flux, and (iii) the maximum potential
root uptake rate. It is important to note that it is not straightforward to implement active and passive uptake processes as discussed by Šimůnek and Hopmans (2009), in which active uptake is only invoked when passive uptake does not meet a maximum nutrient plant demand.
OVERVIEW OF RECENT APPLICATIONSThe early studies of Jacques et al. (2008a, 2008b) illustrated
the coupled effects of water flow dynamics and geochemistryon the mobility and availability of various contaminants. Since
water content variations due to the temporal nature of precipitation, infiltration, and evapotranspiration continually change
concentrations in the aqueous phase, other chemical factors andprocesses such as ionic strength, pH, and aqueous complexation
are affected as well. Differences in component mobility resultin changes in aqueous speciation and, consequently, both217

Diederik Jacques, Jiří Šimůnek, Dirk Mallants, Martinus Th. van Genuchten

Fig. 4. Organic matter decomposition networks: (top) a schemewithout leachable soil organic matter, (bottom) a scheme withleachable soil organic matter. Orange lines indicate decomposition
to easily decomposable soil organic matter (SOM-P1), brown lines
indicate decomposition to slowly decomposable soil organic matter(SOM-P2), green lines indicate assimilation into biomass (Biomass), blue lines indicate respiration and production of inorganiccarbon species (CO2), and red lines indicate decomposition towards
leachable soil organic matter (SOM-L). The dashed black line
represents kinetic exchange between leachable soil organic matterand dissolved organic matter.

mineral dissolution/precipitation and sorption can be affectedby the temporal-spatial variability of water contents and water
fluxes. Comprehensive studies and/or assessments of fate andtransport of contaminants in soil systems need to take into
account these types of variations. HPx already has been appliedin many studies, including for evaluation of Hg fate in soil
systems (Bessinger and Marks, 2010; Leterme and Jacques,
2015; Leterme et al., 2014) or, more recently, scenario analyses
involving the transport of trace metals in a soil leached withcoal seam gas produced water (Mallants et al., 2017). The
benchmark study of Greskowiak et al. (2015), in which HPx
was used as one of several reactive transport models they evaluated, illustrated the applicability of HPx to describe Utransport using a multi-rate mass transfer model (MRMT, e.g.,
Haggerty and Gorelick (1995)).

218

Another benchmark study was used to verify the implementation of changes in the porosity, diffusion coefficient andhydraulic conductivity as a consequence of mineral precipitation and dissolution in multiple reactive transport models (Xieet al., 2015), including HP1. Valdes-Abellan et al. (2017) developed a model (using HP1) to simulate the spatial-temporal
evolution of soil hydraulic properties during non-conventional
water irrigation during non-isothermal transient flow. In
addition to linking the porosity, ah, and the saturated hydraulicconductivity to the amount of precipitated or dissolved minerals, soil hydraulic properties were also linked to changes in thesoil water composition, salinity, and temperature. ValdesAbellan et al. (2017) in their study indicated that a close interplay exist between gypsum-calcite dissolution/precipitation,
evaporation fronts, and the key role of the partial pressure of
CO2. While the partial pressure of CO2 was linked only to calcite equilibria, the role of organic matter degradation (andcorresponding CO2 production) could potentially be incorporated within the HPx framework as well (see below).
HYDRUS in its specialized C-Ride add-on module solves aset of equations describing the transport of colloid and colloidfacilitated solute transport (Šimůnek et al., 2006). A similar setof equations can be implemented in HPx, thus allowing one toconsider the effects of geochemistry on various parameters suchas attachment efficiency or competitive sorption on thetransport of colloids. Zhou et al. (2016) developed a model
simulating the transport of three antibiotics in the presence ofdifferent kinds of natural organics originating from manure.
Compared to the comprehensive capabilities in the standard
HYDRUS code, HPx can account for competitive sorption of
antibiotics to the soil using a Langmuir competitive kinetic
adsorption model in conjunction with competitive kinetic sorption on mobile and immobile organic colloids. The role ofgeochemistry on colloid mobility was described in a study byMakselon et al. (2017) who modeled the transport of silver
nanoparticles in an experiment with flow interruption, evaporation, and ionic strength changes. Ionic strength, and hence thegeneral aqueous geochemical conditions, determined the valueof the attachment efficiency, which can be calculated using
DVLO theory. A lower ionic strength results in a lower attachment efficiency because of a higher energy barrier and a smaller secondary minimum. BASIC scripts in the input file included
root finding algorithms to identify the primary and secondary
minima to obtain the attachment efficiency as a function of thechanging geochemical conditions during the transport experiment.
EXAMPLES
Organic matter dynamics
As an illustration of the flexible framework of HPx to define
reaction networks in general, and for soil organic matter (SOM)
degradation in particular, a reaction network loosely based on a
study by Porporato et al. (2003) was implemented for this example (HPx projects are available at https://www.pcprogress.com/en/Default.aspx?h1d-library). For convenience,the single-porosity model formulation is used in this example.
Figure 4 shows the network with three organic matter pools: (i)a fast decomposable SOM pool (SOM-P1) with concentration
CSOM-P1, (ii) a slow decomposable SOM pool (SOM-P2) withconcentration CSOM-P2, and (iii) a biomass pool (Biomass) withconcentration CB. In addition, inorganic carbon is released in
the aqueous and gaseous soil phases during decomposition, CI.
This decomposition is described in terms of four ordinary
differential equations:

The HPx software for multicomponent reactive transport during variably-saturated flow: Recent developments and applications

d CSOM-P1= iF + d B − dSOM-P1
dt

(16)

d CSOM-P2= rSOM-P1→P2 dSOM-P1 − dSOM-P2
dt

(17)

d CB
= rSOM-P1→B dSOM-P1 + rSOM-P2→ B dSOM-P2 − d B
dt

(18)

kin
Rm,C(4)= rSOM-P1→I dSOM-P1 + rSOM-P2→ I dSOM-P2(19)

where dX refers to the decomposition rates, rx→y denotes thefraction of the decomposition products of pool x going to pooly, and iF denotes the addition of fresh organic matter as litter
fall on top of the soil profile or as root decay. Wutzler andReichstein (2008) discussed different model formulations thatcan be used to include decomposer biomass into the decomposition rates: (i) no representation of the biomass using equations
that are based only on the substrate, (ii) a linear representationof biomass, and (iii) a nonlinear representation of biomass.
These models differ in the scaling behavior with the substrateand how they represent priming effects. Although all of theseformulations could be included in HPx, the example given
below is limited to a non-explicit representation of the biomass
and a linear model of the available substrate, i.e.,d X = k X CX

(20)

where kX is the first-order decomposition term [T–1], CX is theconcentration of the substrate [ML–3], and X denotes a soilorganic matter pool. The value of the rate parameter kX may
depend on temperature fT(T) [–] and soil moisture fM(Sθ) [–] asfollows:

k X = k X ,0 fT (T ) f M ( Sθ )(21)

where kX,0 is the reference rate [T–1]. Sierra et al. (2012) gave an
overview of different functions to account for the effects ofenvironmental factors on the degradation rate. Here we describefM(Sθ) using the equations of Porporato et al. (2003):

 Sθ
S
 fc
f M ( Sθ ) = 
 Sfc
 Sθ

Sθ ≤ Sfc(22)

Sθ > Sfc

where Sθ is the degree of saturation (=θ/θs) [–] and Sfc is thesaturation degree at field capacity [–]. The temperature dependency is described using the formulation from ROTC (Jenkinson
et al., 1990):

f T (T ) =47.9
 106 1 + exp  T + 18.3 (23)Still needed to complete the description of the organic matterdecomposition network is an aqueous speciation model for
inorganic C species. We refer to the Appendix for details.
Table 2 further gives the parameters used in a transient flow
simulation for a meteorological time series representative ofBelgium conditions.

Fig. 5. Profiles of SOM-P1 (top), SOM-P2 (middle), and Biomass(bottom) for network 1.

Figure 5 shows the computed distributions versus depth ofthe various pools of soil organic matter and biomass during fiveconsecutive years. Large changes during the first few years
indicate that selected initial and production terms are not in
long-term equilibrium with the degradation processes. The
current set of parameters and meteorological variables do not
produce much variability in the vertical direction. The relativeenrichment of the top 5 cm is due to SOM-P1 originating bothfrom litter fall and root decay whereas deeper in the profile,
only root decay contributes to the SOM-P1 pool. A morethroughout analysis of the model setup using a sensitivity analysis of the model parameters is beyond the scope of the current
paper but could be done straightforwardly with HPx as was
illustrated by Leterme and Jacques (2015).

219

Diederik Jacques, Jiří Šimůnek, Dirk Mallants, Martinus Th. van Genuchten

Fig. 6. Time series of SOM-P1 (top), SOM-P2 (middle), andBiomass (bottom) for network 1.

Fig. 7. Time series of fM (top), fT (middle), and fMfT (bottom) for
network 1.

Figure 6 shows time series of various pools of SOM and
Biomass at several depths. The strong seasonal fluctuations inthese pools are due to variations in soil water contents andtemperatures, leading to changes in the environmental factors
on the degradation rates as expressed by Eqs. (21)–(23). As
shown in Figure 7, variations in fM and fT are similar at different
soil depths, which explain the relatively uniform distribution ofvarious organic pools (Figure 5). HPx in general provides
geochemical information of the system such as pH (not further
shown here) and inorganic aqueous and gaseous C
concentrations, the latter shown in Figure 8 in terms of thepartial pressure of CO2 in the soil air phase. In this example C
was the only aqueous species and no pH controllingcomponents or processes were included. Figure 8 also shows
the cumulative flux of CO2 from the soil to the atmosphere,
which is a critical variable in soil-plant-atmosphere and regional or global climate studies (e.g., Thaysen et al. (2014)).

220

Dissolved organic matter
The previous organic matter degradation model is further extended in this example by including a leachable organic matterpool (SOM-L) and a dissolved organic matter pool (DOM),
which are both kinetically exchanged. The formulation, loosely
based on the model of Braakhekke et al. (2011), allows for

The HPx software for multicomponent reactive transport during variably-saturated flow: Recent developments and applications

Fig. 9. Time series of DOM (top) and SOM-L (bottom) for network 2.

d CSOM-L= rSOM-P1→SOM-L dSOM-P1 + rSOM-P2→SOM-L dSOM-P2
dt−dSOM-L − dS

(27)

kin
Rm,C(4)= rSOM-P1→I dSOM-P1 + rSOM-P2→I dSOM-P2 + rSOM-L→I dSOM-L (28)kinRm,DOM= ds(29)

where dS is the desorption term:
Fig. 8. Time series of C concentration (mol/kgw) (top), PCO2
(middle) and cumulative CO2 flux to the atmosphere (bottom)
for network 1.
advective-dispersive transport of DOM. The set of governing
ordinary-differential equations is as follows:d CSOM-P1= iF + d B − dSOM-P1
dt

(24)

d CSOM-P2= rSOM-P1→ P2 dSOM-P1 − dSOM-P2
dt

(25)

d CB
= rSOM-P1→B dSOM-P1 + rSOM-P2→ B dSOM-P2
dt+ rSOM-L→ B dSOM-L − d B

(26)

eq
ds = kS (CSOM-L − CSOM-L)= kS  CSOM-L −K L Smax CDOM 
1 + K LCDOM (30)

in which KL is the coefficient of binding affinity [dm3gC–1] andSmax the maximum adsorption capacity [gC dm–3]. Time seriesof the new two pools, SOM-L and DOM, are shown in Figure
9. DOM shows short-scale temporal variations due to changes
in the water content. These variations are less pronounced for
SOM-L due to the kinetics of the sorption process. The advection and dispersion of DOM result in an increase in both DOM
and SOM-L at depths below the initial source of soil organic
matter between 0 and 20 cm below the soil surface. Note thatthe other pools potentially will develop also below the 20-cm
depth when this degradation scheme is used, but their amounts
remain very small for the chosen set of parameters.

221

Diederik Jacques, Jiří Šimůnek, Dirk Mallants, Martinus Th. van Genuchten

As an example, a biodiffusion model was implemented within the input files of HPx:

J s, X = − DB, X

∂C X
∂x

(31)

where Js,X is the solid phase diffusive flux of the Xth pool [molT–1 L–2], DB,X is the effective solid phase diffusion coefficient[L2 T–1], and CX is the solid phase concentration [mol L–3]. The
biodiffusion model was applied to organic network 1 describedabove with SOM-P1, SOM-P2, and biomass. We assumed
rather arbitrarily that DB,X = 10–10 m2/s. As shown in Figure 10,
a smooth profile of the three organic pools developed over thecourse of 3 years to about 15 cm below the initial layer withorganic matter.
CONCLUSIONS

Fig. 10. Profiles of SOM-S1 (top), SOM-S2 (middle) and Biomass(bottom) for a simulation with network 1 including biodiffusion.

Bioturbation
The long term prediction of soil carbon cycling requires an
explicit representation of the vertical distribution of soil organic
matter, and hence a model able to represent processes leading
to the vertical distribution of SOM such as the root litter distribution and the transport of soil organic matter (Braakhekke etal., 2011). The biological transport of soil solid phases, often
referred to as bioturbation, depends strongly on the behavior of
the biological community species, which necessitated the development of different modeling approaches (Meysman et al.,2003a). Models based on Fick's law are local random models(Meysman et al., 2005), whereas non-local models need to be
employed when ingestion and egestion processes of SOM by
bio fauna are vertically separated (Jarvis et al., 2010).

222

The original HP1 model (Jacques and Šimůnek, 2005;
Jacques et al., 2006) has seen many extensions during the past
10 years or more. The most important changes were the development of two- and three-dimensional versions (HPx), and theincorporation of the software PHREEQC-3, OpenMP, enhanced
scripting possibilities, diffusion in the gas phase of reactive
gases in the software, and linking geochemical state variables totransport properties. Using the terminology of Meysman et al.
(2003b), HPx is now moving towards a more flexible and extendible problem-solving environment to allow users to incorporate a range of user-defined conceptual models to tackle their
specific projects. Flexibility and extendibility are discussed inthis manuscript with respect to coupling geochemistry withtransport properties, colloid and colloid-affected transport,
bioturbation, and root solute uptake. Recent applications showthat HPx, and similar reactive transport codes, are useful tools
to evaluate many environmental and engineering problems bothfrom a scientific as well as from an application point of view,thus underpinning the main statement from Li et al. (2017).
In this paper, we illustrated the capabilities and flexibility ofHPx by developing a reaction-transport model for soil organic
matter, in which we included different factors important in soilorganic matter degradation. Only a few SOM degradation models can similarly consider spatial-temporal variations in watercontents, water fluxes, and temperatures, as well as different
conceptual models for SOM degradation networks, advectiondispersive transport of DOM, CO2 transport, and bioturbation in
one modeling framework. Within the scope of this paper, therole of microbial processes together with the transport of oxygen were not addressed, but conceptually could be implemented
in a very similar way. We also acknowledge that further parameter sensitivity studies and validation of these models are crucial in further model applications.
Acknowledgement. The first author acknowledges the supportof the Belgian Nuclear Research Centre for developing, maintaining, and applying the HPx code.
REFERENCES
Aagaard, P., Helgeson, H.C., 1982. Thermodynamic and kinetic
constraints on reaction rates among minerals and aqueous
solutions, 1. Theoretical considerations. Am. J. Sci., 282,
237–285.
Amos, R.T., Mayer, K.U., 2006. Investigating the role of gas
bubble formation and entrapment in contaminated aquifers:
Reactive transport modelling. J. Cont. Hydrol., 87, 123–154.

The HPx software for multicomponent reactive transport during variably-saturated flow: Recent developments and applications

Appelo, C.A.J., Parkhurst, D.L., Post, V.E.A., 2014. Equations
for calculating hydrogeochemical reactions of minerals andgases such as CO2 at high pressures and temperatures.
Geochim. Cosmochim. Acta, 125, 49–67.
Appelo, C.A.J., Wersin, P., 2007. Multicomponent diffusion
modeling in clay systems with application to the diffusion oftritium, iodide, and sodium in Opalinus Clay. Env. Sci.
Tech., 41, 5002–5007.
Batlle-Aguilar, J., Brovelli, A., Porporato, A., Barry, D.A.,
2011. Modelling soil carbon and nitrogen cycles during land
use change: A review. Agron. Sust. Devel., 31, 251–274.
Bennacer, L., Ahfir, N.D., Alem, A., Wang, H.Q., 2017.
Coupled effects of ionic strength, particle size, and flow
velocity on transport and deposition of suspended particles
in saturated porous media. Transp. Por. Med., 118, 251–269.
Bessinger, B.A., Marks, C.D., 2010. Treatment of mercurycontaminated soils with activated carbon: A laboratory,
field, and modeling study. Remed., 21, 115–135.
Bloom, S.A., Mansell, R.S., 2001. An algorithm for generating
cation exchange isotherms from binary selectivity
coefficients. Soil Sci. Soc. Am. J., 65, 1426–1429.
Bond, W.J., 1995. On the Rothmund-Kornfeld description of
cation exchange. Soil Sc. Soc. Am. J., 59, 436–443.
Borkovec, M., Westall, J., 1983. Solution of the PoissonBoltzmann equation for surface excesses of ions in the
diffuse layer at the oxide electrolyte interface. J. Elect.
Chem., 150, 325–337.
Bozorg, A., Gates, I.D., Sen, A., 2015a. Impact of biofilm on
bacterial transport and deposition in porous media. J. Cont.
Hydrol.,183, 109–120.
Bozorg, A., Gates, I.D., Sen, A., 2015b. Using bacterial
bioluminescence to evaluate the impact of biofilm on porous
media hydraulic properties. J. Microb. Meth., 109, 84–92.
Braakhekke, M.C., Beer, C., Hoosbeek, M.R., Reichstein, M.,
Kruijt, B., Schrumpf, M., Kabat, P., 2011. SOMPROF: A
vertically explicit soil organic matter model. Ecol. Mod.,
222, 1712–1730.
Brantley, S., Goldhaber, M.B., Ragnarsdottir, K.V., 2007.
Crossing disciplines and scales to understand the criticalzone. Elements, 3, 307–314.
Brooks, R.H., Corey, A., 1964. Hydraulic properties of porous
media. Hydrol. Paper No. 3, Colorado State Univ., Fort
Collins, CO.
Carles Brangarí, A., Sanchez-Vila, X., Freixa, A., M. Romaní,A., Rubol, S., Fernàndez-Garcia, D., 2017. A mechanistic
model (BCC-PSSICO) to predict changes in the hydraulicproperties for bio-amended variably saturated soils. Water
Resour. Res., 53, 93–109.
Charlton, S.R., Parkhurst, D.L., 2011. Modules based on thegeochemical model PHREEQC for use in scripting andprogramming languages. Comp. & Geosc., 37, 1653–1663.
Durner, W., 1994. Hydraulic conductivity estimation for soils
with heterogeneous pore structure. Water Resour. Res., 30,
211–223.
Dzombak, D.A., Morel, F.M.M., 1990. Surface Complexation
Modeling – Hydrous Ferric Oxide. John Wiley, New York.
Freedman, V.L., Bacon, D.H., Saripalli, K.P., Meyer, P.D.,
2004. A film depositional model of permeability for mineralreactions in unsaturated media. Vadose Zone J., 3, 1414–
1424.
Greskowiak, J., Gwo, J., Jacques, D., Yin, J., Mayer, K.U.,
2015. A benchmark for multi-rate surface complexation and
1D dual-domain multi-component reactive transport ofU(VI). Comp. Geosc., 19, 585–597.

Guggenheim, E.A., 1937. Theoretical basis of Raoult's law.
Trans. Faraday Soc., 33, 151–159.
Haggerty, R., Gorelick, S.M., 1995. Multiple-rate mass transfer
for modeling diffusion and surface reactions in media withpore-scale heterogeneity. Water Resour. Res., 31, 2383–
2400.
Hiemstra, T., VanRiemsdijk, W.H., 1996. A surface structural
approach to ion adsorption: The charge distribution (CD)model. J. Col. Int. Sci., 179, 488–508.
Hommel, J., Lauchnor, E., Phillips, A., Gerlach, R.,
Cunningham, A.B., Helmig, R., Ebigbo, A., Class, H., 2015.
A revised model for microbially induced calcite
precipitation: Improvements and new insights based onrecent experiments. Water Resour. Res., 51, 3695–3715.
Jacques, D., 2009. Benchmarking of the cement model anddetrimental chemical reactions including temperaturedependent parameters. Project near surface disposal ofcategory A waste at Dessel, NIRAS-MP5-03 DATA-LT(NF)
Version 1.
Jacques, D., Šimůnek, J., 2005. User Manual of theMulticomponent Variably- Saturated Flow and Transport
Model HP1. SCK•CEN-BLG-998.
Jacques, D., Šimůnek, J., Mallants, D., van Genuchten, M.T.,
2006. Operator-splitting errors in coupled reactive transport
codes for transient variably saturated flow and contaminanttransport in layered soil profiles. J. Contam. Hydrol., 88,
197–218.
Jacques, D., Šimůnek, J., Mallants, D., van Genuchten, M.T.,
2008a. Modeling coupled hydrologic and chemical
processes: Long-term uranium transport following
phosphorus fertilization. Vadose Zone J., 7, 698–711.
Jacques, D., Šimůnek, J., Mallants, D., van Genuchten, M.T.,
2008b. Modelling coupled water flow, solute transport andgeochemical reactions affecting heavy metal migration in a
podzol soil. Geoderma, 145, 449–461.
Jacques, D., Smith, C., Simůnek, J., Smiles, D., 2012. Inverse
optimization of hydraulic, solute transport, and cationexchange parameters using HP1 and UCODE to simulate
cation exchange. J. Contam. Hydrol., 142–143, 109–125.
Jarvis, N.J., Taylor, A., Larsbo, M., Etana, A., Rosen, K., 2010.
Modelling the effects of bioturbation on the re-distributionof 137Cs in an undisturbed grassland soil. Eur. J. Soil Sci.,
61, 24–34.
Jenkinson, D.S., Andrew, S.P.S., Lynch, J.M., Goss, J.M.,
Tinker, P.B., 1990. The turnover of organic carbon andnitrogen in soil. Philosoph. Trans., 329, 361–368.
Kosugi, K., 1996. Lognormal distribution model forunsaturated soil hydraulic properties. Water Resour. Res.,32, 2697–2703.
Laliberté, M., 2007. Model for calculating the viscosity of
aqueous solutions. J. Chem. Eng. Data, 52, 321-335.
Laliberté, M., Cooper, W.E., 2004. Model for calculating thedensity of aqueous electrolyte solutions. J. Chem. Eng. Data,49, 1141–1151.
Leterme, B., Blanc, P., Jacques, D., 2014. A reactive transportmodel for mercury fate in soil—application to different
anthropogenic pollution sources. Environ. Sci. Pollut. Res.,
12279–12293.
Leterme, B., Jacques, D., 2015. A reactive transport model for
mercury fate in contaminated soil-sensitivity analysis.
Environ. Sci. Pollut. Res., 22, 16830–16842.
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.,

223

Diederik Jacques, Jiří Šimůnek, Dirk Mallants, Martinus Th. van Genuchten

Zachara, J., Kocar, B., McIntosh, J., Tutolo, B.M., Kumar,M., Sonnenthal, E., Bao, C., Beisman, J., 2017. Expanding
the role of reactive transport models in critical zone
processes. Earth Sci. Rev., 165, 280–301.
Liu, S., Jacques, D., Govaerts, J., Wang, L., 2014. Conceptualmodel analysis of interaction at a concrete-Boom Clayinterface. Phys. Chem. Earth, 70–71, 150–159.
Maes, N., Wang, L., Hicks, T., Bennett, D., Warwick, P., Hall,
T., Walker, G., Dierckx, A., 2006. The role of naturalorganic matter in the migration behaviour of americium inthe Boom Clay - Part I: Migration experiments. Phys. Chem.
Earth, 31, 541–547.
Makselon, J., Zhou, D., Engelhardt, I., Jacques, D., Klumpp, E.,
2017. Experimental and numerical investigations of silver
nanoparticle transport under variable flow and ionic strength
in soil. Envir. Sci. Techn., 51, 2096–2104.
Mallants, D., Šimůnek, J., van Genuchten, M.T., Jacques, D.,2017. Simulating the fate and transport of coal seam gas
chemicals in variably-saturated soils using HYDRUS.
Water, 9, 6, 385.
Manzoni, S., Porporato, A., 2009. Soil carbon and nitrogenmineralization: Theory and models across scales. Soil Biol.
Biochem., 41, 1355–1379.
Martinez, B.C., DeJong, J.T., Ginn, T.R., 2014. Biogeochemical reactive transport modeling of microbial
induced calcite precipitation to predict the treatment of sandin one-dimensional flow. Comp. Geotech., 58, 1–13.
Mayer, K.U., Alt-Epping, P., Jacques, D., Arora, B., Steefel,
C.I., 2015. Benchmark problems for reactive transportmodeling of the generation and attenuation of acid rockdrainage. Comp. Geosci., 19, 599–611.
Mayer, K.U., Frind, E.O., Blowes, D.W., 2002.
Multicomponent reactive transport modeling in variably
saturated porous media using a generalized formulation for
kinetically controlled reactions. Water Resour. Res., 38,
1174, DOI: 1110.1029/2001WR000862.
Mays, D.C., Hunt, J.R., 2007. Hydrodynamic and chemicalfactors in clogging by montmorillonite in porous media.
Envir. Sci. Techn., 41, 5666–5671.
Meysman, F.J.R., Boudreau, B.P., Middelburg, J.J., 2003a.
Relations between local, nonlocal, discrete and continuous
models of bioturbation. J. Mar. Res., 61, 391–410.
Meysman, F.J.R., Boudreau, B.P., Middelburg, J.J., 2005.
Modeling reactive transport in sediments subject tobioturbation and compaction. Geochim. Cosmochim. Acta,69, 3601–3617.
Meysman, F.J.R., Middelburg, J.J., Herman, P.M.J., Heip,
C.H.R., 2003b. Reactive transport in surface sediments. I.
Model complexity and software quality. Comp. Geosci., 29,
291–300.
Millington, R.J., Quirk, J.P., 1961. Permeability of porous
solids. Trans. Faraday Soc., 57, 1200–1206.
Nowack, B., Mayer, K.U., Oswald, S.E., van Beinum, W.,
Appelo, C.A.J., Jacques, D., Seuntjens, P., Gérard, F.,
Jaillard, B., Schnepf, A., Roose, T., 2006. Verification andintercomparison of reactive transport codes to describe rootuptake. Plant and Soil, 285, 305–321.
Or, D., Smets, B.F., Wraith, J.M., Dechesne, A., Friedman,
S.P., 2007. Physical constraints affecting bacterial habitats
and activity in unsaturated porous media – a review. Adv.
Wat. Res., 30, 1505–1527.
Paradelo, M., Perez-Rodriguez, P., Fernandez-Calvino, D.,
Arias-Estevez, M., Lopez-Periago, J.E., 2012. Coupledtransport of humic acids and copper through saturated
porous media. Eur. J. Soil Sci., 63, 708–716.

224

Parkhurst, D.L., Appelo, C.A.J., 2013. Description of Input andExamples for PHREEQC Version 3 – A Computer Programfor Speciation, Batch-Reaction, One-Dimensional Transport,
and Inverse Geochemical Calculations. Chapter 43 ofSection A. Groundwater Book 6, Modeling Techniques.
Parkhurst, D.L., Wissmeier, L., 2015. PhreeqcRM: A reaction
module for transport simulators based on the geochemical
model PHREEQC. Adv. Wat. Res., 83, 176–189.
Patel, R., Phung, Q.T., Seetharam, S.C., Perko, J., Jacques, D.,
Maes, N., De Schutter, G., Ye, G., van Breugel, K., 2016.
Diffusivity of saturated ordinary Portland cement- basedmaterials: A critical review of experimental and analyticalmodelling approaches. Cem. Con. Res., 90, 52–72.
Peng, D.-Y., Robinson, D.B., 1976. A new two-constantequation of state. Ind. Eng. Chem. Fund., 15, 59–64.
Phung, Q.T., Maes, N., Jacques, D., De Schutter, G., Ye, G.,
Perko, J., 2016. Modelling the carbonation of cement pastes
under a CO2 pressure gradient considering both diffusiveand convective transport. Const. Build. Mat., 114, 333–351.
Porporato, A., D'Odorico, P., Laio, F., Rodriguez-Iturbe, I.,
2003. Hydrologic controls on soil carbon and nitrogen
cycles. I. Modeling scheme. Adv. Wat. Res., 26, 45–58.
Puigdomènech, I., Rard, J.A., Plyasunov, A.V., Grenthe, I.,1997. Temperature corrections to thermodynamic data andenthalpy calculations. In: Grenthe, I., PuigdomSnech, I.
(Eds.): OECD Nuclear Chemistry, Paris, pp. 427–493.
Riley, W.J., Maggi, F., Kleber, M., Torn, M.S., Tang, J.Y.,
Dwivedi, D., Guerry, N., 2014. Long residence times of
rapidly decomposable soil organic matter: application of amulti-phase, multi-component, and vertically resolved model(BAMS1) to soil carbon dynamics. Geosci. Model Dev., 7,1335–1355.
Rockhold, M.L., Yarwood, R.R., Niemat, M.R., Bottomley,
P.J., Selker, J.S., 2002. Considerations for modelingbacterial-induced changes in hydraulic properties of variably
saturated porous media. Adv. Wat. Res., 25, 477–495.
Schmidt, M.W.I., Torn, M.S., Abiven, S., Dittmar, T.,
Guggenberger, G., Janssens, I.A., Kleber, M., KogelKnabner, I., Lehmann, J., Manning, D.A.C., Nannipieri, P.,
Rasse, D.P., Weiner, S., Trumbore, S.E., 2011. Persistenceof soil organic matter as an ecosystem property. Nature, 478,
49–56.
Seuntjens, P., Nowack, B., Schulin, R., 2004. Root-zone
modeling of heavy metal uptake and leaching in the presenceof organic ligands. Plant and Soil, 265, 61–73.
Sharqawy, M.H., Lienhard V, J.H., Zubair, S.M., 2010.
Thermophysical properties of seawater: A review of existingcorrelations and data. Desal. Wat. Treat., 16, 354–380.
Sierra, C.A., Müller, M., Trumbore, S.E., 2012. Models of soilorganic matter decomposition: the SoilR package, version
1.0. Geosci. Model Dev., 5, 1045–1060.
Silberbush, M., Ben-Asher, J., Ephrath, J.E., 2005. A model for
nutrient and water flow and their uptake by plants grown in asoilless culture. Plant and Soil, 271, 309–319.
Šimůnek, J., He, C., Pang, L., Bradford, S.A., 2006. Colloidfacilitated solute transport in variably saturated porous
media: Numerical model and experimental verification.
Vadose Zone J., 5, 1035–1047.
Šimůnek, J., Hopmans, J.W., 2009. Modeling compensated root
water and nutrient uptake. Ecol. Mod., 220, 505–521.
Šimůnek, J., Jacques, D., Šejna, M., van Genuchten, M.T.,
2012. The HP2 Program for HYDRUS (2D/3D). A Coupled
Code for Simulating Two-Dimensional Variably-Saturated
Water Flow, Head Transport, Solute Transport Flow, andThe HPx software for multicomponent reactive transport during variably-saturated flow: Recent developments and applications

Biogeochemistryin
PorousMedia.
(HYDRUS+PHREEQC+2D), Version 1.
Simunek, J., Sejna, M., Saito, H., Sakai, K., van Genuchten,
M.T., 2013. The HYDRUS-1D Software Package forSimulating the Movement of Water, Heat, and Multiple
Solutes in Variably Saturated Media. Version 4.17.
HYDRUS Software Series 3. Department of Environmental
Sciences, University of California Riverside, Riverside,
California, USA.
Šimůnek, J., van Genuchten, M.T., 2008. Modeling
nonequilibrium flow and transport processes using
HYDRUS. Vadose Zone J., 7, 782–797.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent
developments and applications of the HYDRUS Computer
Software Packages. Vadose Zone J., 15, DOI:
10.2136/vzj2016.04.0033.
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., Yeh, G.T., 2015.
Reactive transport codes for subsurface environmental
simulation. Comp. Geosci., 19, 445–478.
Steefel, C.I., DePaolo, D.J., Lichtner, P.C., 2005. Reactive
transport modeling: An essential tool and a new research
approach for the Earth sciences. Earth Plan. Sci. Let., 240,
539–558.
Tang, J., Riley, W.J., 2015. Weaker soil carbon-climate
feedbacks resulting from microbial and abiotic interactions.
Nature Clim. Change, 5, 56–60.
Tang, J.Y., Riley, W.J., Koven, C.D., Subin, Z.M., 2013.
CLM4-BeTR, a generic biogeochemical transport andreaction module for CLM4: model development, evaluation,
and application. Geosci. Model Dev., 6, 127–140.
Thaysen, E.M., Jacques, D., Jessen, S., Andersen, C.E., Laloy,
E., Ambus, P., Postma, D. and Jakobsen, I., 2014. Inorganiccarbon fluxes across the vadose zone of planted
andunplanted soil mesocosms. Biogeosci., 11, 7179–7192.
Todd-Brown, K.E.O., Randerson, J.T., Post, W.M., Hoffman,
F.M., Tarnocai, C., Schuur, E.A.G., Allison, S.D., 2013.
Causes of variation in soil carbon simulations from CMIP5
Earth system models and comparison with observations.
Biogeosci., 10, 1717–1736.
Valdes-Abellan, J., Jiménez-Martínez, J., Candela, L., Jacques,D., Kohfahl, C., Tamoh, K., 2017. Reactive transportmodelling to infer changes in soil hydraulic properties
induced by non-conventional water irrigation. J. Hydrol.,549, 114–124.
van Genuchten, M.T., 1980. Closed-form equation forpredicting the hydraulic conductivity of unsaturated soils.
Soil Sci. Soc. Am. J., 44, 892–898.
Vereecken, H., Schnepf, A., Hopmans, J.W., Javaux, M., Or,
D., Roose, T., Vanderborght, J., Young, M.H., Amelung,W., Aitkenhead, M., Allison, S.D., Assouline, S., Baveye,
P., Berli, M., Brüggemann, N., Finke, P., Flury, M., Gaiser,T., Govers, G., Ghezzehei, T., Hallett, P., Hendricks
Franssen, H.J., Heppell, J., Horn, R., Huisman, J.A.,
Jacques, D., Jonard, F., Kollet, S., Lafolie, F., Lamorski, K.,
Leitner, D., McBratney, A., Minasny, B., Montzka, C.,
Nowak, W., Pachepsky, Y., Padarian, J., Romano, N., Roth,K., Rothfuss, Y., Rowe, E.C., Schwen, A., Šimůnek, J.,
Tiktak, A., Van Dam, J., van der Zee, S.E.A.T.M., Vogel,
H.J., Vrugt, J.A., Wöhling, T., Young, I.M., 2016. Modelingsoil processes: Review, key challenges, and new
perspectives. Vadose Zone J., 15, 1–57.

Vogel, T., Cislerova, M., Hopmans, J.W., 1991. Porous media
with linearly variable hydraulic properties. Water Resour.
Res., 27, 2735–2741.
Wissmeier, L., Barry, D.A., 2009. Effect of mineral reactions
on the hydraulic properties of unsaturated soils: Modeldevelopment and application. Adv. Wat. Res., 32, 1241–
1254.
Wissmeier, L., Barry, D.A., 2010. Implementation of variably
saturated flow into PHREEQC for the simulation ofbiogeochemical reactions in the vadose zone. Env. Mod.
Soft., 25, 526–538.
Wutzler, T., Reichstein, M., 2008. Colimitation of
decomposition by substrate and decomposers - a comparisonof model formulations. Biogeosci., 5, 749–759.
Xie, M., Mayer, K.U., Claret, F., Alt-Epping, P., Jacques, D.,
Steefel, C., Chiaberge, C., Šimůnek, J., 2015.
Implementation and evaluation of permeability-porosity andtortuosity-porosity relationships linked to mineraldissolution-precipitation. Comp. Geosci., 19, 655–671.
Yarwood, R.R., Rockhold, M.L., Niemet, M.R., Selker, J.S.,
Bottomley, P.J., 2006. Impact of microbial growth on waterflow and solute transport in unsaturated porous media. Water
Resour. Res., 42, W10405, 1–11.
Yu, C., Muñoz-Carpena, R., Gao, B., Perez-Ovilla, O., 2013.
Effects of ionic strength, particle size, flow rate, andvegetation type on colloid transport through a dense
vegetation saturated soil system: Experiments and modeling.
J. Hydrol., 499, 316–323.
Zhou, D., Thiele-Bruhn, S., Arenz-Leufen, M.G., Jacques, D.,
Lichtner, P., Engelhardt, I., 2016. Impact of manure-related
DOM on sulfonamide transport in arable soils. J. Contam.
Hydrol., 192, 118–128.
Received 12 July 2017
Accepted 21 September 2017

APPENDIX – Details about the organic matter degradation
model.
The aquatic geochemical reactions within the adopted modeling approach are limited to the inorganic carbon reactions, as
given in Table 1. For all aqueous species, the activity correction
model of WATEQ is used:

log γ =− Az 2 µ+ bµ
1 + Ba µ(1)

where γ is the activity correction coefficient, z is the charge, µis the ionic strength, A and B are constants at a given temperature, and a and b are species-dependent parameters, except for
CO2 for which the activity coefficient is obtained as log γ = 0.1µ. The temperature dependence of the mass action constant isgiven as:log10 K = A1 + A2T +A3
A+ A4 log10 T + 52 + A6T 2
T
T

(2)

where A1, …, A6 are parameters and T is temperature in Kelvin.
The relevant temperature range is limited between ~0 to ~30°C.
Tables 2 through 5 provide a list of parameter definitionsand values used in the simulations of the organic matter degradation networks.

225

Diederik Jacques, Jiří Šimůnek, Dirk Mallants, Martinus Th. van Genuchten
Table 1. Overview of species, geochemical reactions, and thermodynamic parameters (values are from the phreeqc.dat database as incorporated in PHREEQC-3).
Species
H+
CO2-3
OHHCO-3
CO2
CO2(g)

Reaction

H2O = OH- + H+
CO2-3 + H+ = HCO-3
CO2-3 + 2H+ = CO2 + H2O
CO2 = CO2

Activity correction
a
b
9.0
0
5.4
0
3.5
0
5.4
0
–
–
–
–

A1

A2

293.29107.89464.20119.870.13610.0325
0.09342.19e-2

Temperature dependence
A3
A4A5–10576.9–5151.8–26986.2–7337.8–123.73–38.93–165.76–44.77A6

0563713.9
2248628.9669371–6.996e–5

Table 2. Parameter values for network 1 (Eqs. (16)–(19)). Parameter values are for illustrative purpose.
SymboliF
kB
rSOM-P1→P2
rSOM-P1→B
rSOM-P1→I
rSOM-P2→B
rSOM-P2→IρbkSOM-P1
kSOM-P2

Parameter description
Addition of fresh organic matter
Decay rate for biomass pool
Transformation factor from SOM-P1 to SOM-P2
Transformation factor from SOM-P1 to B
Transformation factor from SOM-P1 to I
Transformation factor from SOM-P2 to B
Transformation factor from SOM-P2 to I
Bulk density
Rate constant for SOM-P1 degradation
Rate constant for SOM-P2 degradation

Value
0.0015
8.5×10–3
0.1
0.3
0.6
0.4
0.6
1.43
3.00×10–2
5.00×10–3

Unit
gC/d/dm3 bulk
1/d
–
–
–
–
–
kg/dm3 bulk
1/d
1/dTable 3. Parameter values for network 2 (Eqs. (24)–(30)). Parameter values are for illustrative purpose.
SymboliF
kB
rSOM-P1→P2
rSOM-P1→B
rSOM-P1→I
rSOM-P1→L
rSOM-P2→B
rSOM-P2→I
rSOM-P2→L
rSOM-L→B
rSOM-L→I
KL
Smax

ρb

ksor
kSOM-P1
kSOM-P2
kSOM-L

Parameter description
Addition of fresh organic matter
Decay rate for biomass pool
Transformation factor from SOM-P1 to SOM-P2
Transformation factor from SOM-P1 to B
Transformation factor from SOM-P1 to I
Transformation factor from SOM-P1 to SOM-L
Transformation factor from SOM-P2 to B
Transformation factor from SOM-P2 to I
Transformation factor from SOM-P2 to SOM-L
Transformation factor from SOM-P2 to B
Transformation factor from SOM-P2 to I
Binding affinity of DOM
Binding capacity for DOM
Bulk density
First-order kinetic sorption term
Rate constant for SOM-P1 degradation
Rate constant for SOM-P2 degradation
Rate constant for SOM-L degradation

Value
0.0015
8.5×10–3
0.1
0.3
0.55
0.05
0.4
0.5
0.1
0.1
0.9
5
0.4
1.43
0.05
3.0×10–2
5.0×10–3
2.0×10–6

Unit
gC/d/dm3 bulk
gC/d/dm3 bulk
–
–
–
–
–

–
–
–
dm3/gCgC/dm3 bulk
kg/dm3
d–1
1/d
1/d
1/d

Table 4. Initial organic matter pools.
Organic Pool
SOM-P1
SOM-P2
SOM-L1*$
B
DOM$$
only for reaction network 2

Value
1.75×10–4
7.00×10–5
0.019
7.00×10–5
0.01*
Calculated in equilibrium with DOM

Unit
kgC/kg soil
kgC/kg soil
gC/dm3 bulk
kgC / kg soil
gC/dm3soil

Table 5. Soil hydraulic parameters for the van Genuchten (1980) soil hydraulic functions, solute transport parameters, and plant uptake
parameters (for the S-shape water uptake reduction model – see the HYDRUS-1D manual for more information Simunek et al. (2013)) for
single porosity model implementation.
Symbol

θr
θs
α
n
Ks
w
D
a
D

λ
P50
P3

226

Parameter description
Residual water content
Saturated water content
Shape parameter of the water retention function
Shape parameter of the water retention functions
Saturated hydraulic conductivity
Diffusion coefficient in the aqueous phase
Diffusion coefficient in the air phase
Dispersivity
Pressure head at which 50% reduction in water uptake occur
Exponent in the root water uptake response function

Value
0.08
0.45
0.05
1.6
25
2
20,000
10
–500
3

Unit
dm3/dm3
dm3/dm3
cm–1–
cm day–1
cm2 day–1
cm2 day–1cmcm–