water
ArticleSimulating the Fate and Transport of Coal Seam Gas
Chemicals in Variably-Saturated Soils
Using HYDRUS
Dirk Mallants 1, *, Jirka Šimůnek 2 , Martinus Th. van Genuchten 3,4 and Diederik Jacques 5
1
2
3
4
5

*

CSIRO, Adelaide, SA 5000, Australia
Department of Environmental Sciences, University of California Riverside,
Riverside, CA 92521, USA; jsimunek@ucr.edu
Department of Earth Sciences, Utrecht University, 3512 JE Utrecht, The Netherlands;
rvangenuchten@hotmail.com
NIDES Interdisciplinary Center for Social Development, Federal University of Rio de Janeiro,
Rio de Janeiro 21941-901, Brazil
Belgian Nuclear Research Centre, 2400 Mol, Belgium; diederik.jacques@sckcen.be
Correspondence: dirk.mallants@csiro.au; Tel.: +61-8-8303-8595

Academic Editors: Boris Faybishenko and Keith Smettem
Received: 16 April 2017; Accepted: 22 May 2017; Published: 30 May 2017

Abstract: The HYDRUS-1D and HYDRUS (2D/3D) computer software packages are widelyused finite element models for simulating the one-, and two- or three-dimensional movementof water, heat, and multiple solutes in variably-saturated media, respectively. While the standard
HYDRUS models consider only the fate and transport of individual solutes or solutes subject tofirst-order degradation reactions, several specialized HYDRUS add-on modules can simulate far
more complex biogeochemical processes. The objective of this paper is to provide a brief overview
of the HYDRUS models and their add-on modules, and to demonstrate possible applications of
the software to the subsurface fate and transport of chemicals involved in coal seam gas extraction
and water management operations. One application uses the standard HYDRUS model to evaluate
the natural soil attenuation potential of hydraulic fracturing chemicals and their transformation
products in case of an accidental release. By coupling the processes of retardation, first-orderdegradation and convective-dispersive transport of the biocide bronopol and its degradation products,
we demonstrated how natural attenuation reduces initial concentrations by more than a factor ofhundred in the top 5 cm of the soil. A second application uses the UnsatChem module to explore thepossible use of coal seam gas produced water for sustainable irrigation. Simulations with different
irrigation waters (untreated, amended with surface water, and reverse osmosis treated) provided
detailed results regarding chemical indicators of soil and plant health, notably SAR, EC and sodium
concentrations. A third application uses the HP1 module to analyze trace metal transport involving
cation exchange and surface complexation sorption reactions in a soil leached with coal seam gas
produced water following some accidental water release scenario. Results show that the main process
responsible for trace metal migration in soil is complexation of naturally present trace metals withinorganic ligands such as (bi)carbonate that enter the soil upon infiltration with alkaline produced
water. The examples were selected to show how users can tailor the required model complexity to
specific needs, such as for rapid screening or risk assessments of various chemicals nder generic soil
conditions, or for more detailed site-specific analyses of actual subsurface pollution problems.
Keywords: HYDRUS; coal seam gas; contaminant transport; trace elements; HP1

Water 2017, 9, 385; doi:10.3390/w9060385

www.mdpi.com/journal/water

Water 2017, 9, 385

2 of 33

1. Introduction
Many agricultural, industrial, mining, cultural, environmental and other activities are leading to
the production and use of an unprecedented number of chemicals, many of which are intentionallyor accidently released into the environment. For example, agricultural operations now depend on
a broad range of chemicals in plant and animal production such as fertilizers, pesticides, hormones,
pharmaceuticals, pathogenic microorganisms, or fumigants, which have made agriculture one of themost important sources for non-point source pollution [1]. Examples of environmental problems
resulting from non-agricultural activities are the disposal of brine/saline waters produced during
extraction of coalbed or shale gas [2], releases from municipal or industrial landfills [3], acidmine drainage [4], contamination from nuclear activities [5,6], the formation of redox zones inorganic-contaminated aquifers [7], and pollutant release from reactive permeable barriers used for
aquifer remediation [8,9]). Advanced numerical models are required to evaluate the fate and transportof the different types of chemicals that are finding their way into the environment, including for risk
assessments and site-specific management or remediation analyses [10].
The HYDRUS software packages [11] are among many other mathematical models [12] of varying
degrees of complexity and dimensionality that have been developed over the last three to four decades
for evaluating the fate and transport of various chemicals involved. The HYDRUS-1D and HYDRUS(2D/3D) computer software packages are popularly used finite element models for simulating the one-,and two- or three-dimensional movement of water, heat, and multiple solutes in variably-saturated
media, respectively. The models, and their specialized add-on modules, have provided a great
flexibility in simulating the transport of various solutes.
The standard HYDRUS models can be used to simulate the transport of individual solutes whose
behavior is independent of that of other solutes that may be present in the soil solution, or of solutes
subject to sequential first-order degradation reactions. Typical examples of the latter include heavy
metals [13], radionuclides [14,15]), mineral nitrogen species [16], pesticides (e.g., [17]), chlorinatedaliphatic hydrocarbons (e.g., [18]), hormones [19], antibiotics [20], and explosives [21]. While thisapproach has proved to be suitable for many chemicals, important environmental problems often
require analyses of the coupled transport of multiple chemical species that may mutually interact,
create complexed species, precipitate, dissolve, and/or could compete with each other for sorption
sites. Several specialized HYDRUS add-on modules, such as UnsatChem, HP1, C-Ride, Wetland,or Fumigant, were developed to consider these more complex situations [11].
For example, the UnsatChem module [22,23] can simulate the transport and production of carbon
dioxide, and the fate and transport of major ions, while also considering precipitation/dissolutionof mineral phases and/or cation exchange reactions. On the other hand, the HP1 module, which
was developed by coupling the HYDRUS and PHREEQC [24] models, provides much flexibility inmodelling interactions of solutes with minerals, gases, exchangers and sorption surfaces based onthermodynamic equilibrium, kinetic, or mixed equilibrium-kinetic reactions [14,25]. Other specialized
HYDRUS add-on modules, such as Wetland, C-Ride, and Fumigant, allow additional flexibility
in addressing biogeochemical processes in wetlands [26–28], colloid-facilitated solute transport
problems [29–31]), or the fate of fumigants in agricultural soils [32,33], respectively.
An important environmental application for numerical models is the evaluation of the fate
and transport of various solutes involved in unconventional gas operations, such as the disposal ofbrine/saline waters produced during mining of coalbed (coal seam) or shale gas, the accidental releaseof produced water from water holding ponds, and the beneficial use of produced water for irrigation.
In the case of coal seam gas hydraulic fracturing operations, many different chemicals are involved,
the fate of which needs to be evaluated prior to their use to confirm negligible harm to humans andthe environment [2]. These chemicals vary tremendously in terms of their chemical properties andhow they may interact among themselves and/or with the environment. Their fate hence can beevaluated only using numerical models that have the capability to consider all pertinent reactions andinteractions affecting these chemicals.

Water 2017, 9, 385

3 of 33

Our main objective with this study was to demonstrate how the HYDRUS models can be used toevaluate complex processes and interactions affecting the subsurface transport of various chemicals
involved in coal seam gas extraction and water management operations. For this purpose we first
briefly review the capabilities of the HYDRUS models, including their specialized add-on modules,
to simulate the fate and transport of various solutes and their mutual interactions. We will limit
this review to only those add-on modules that will be used subsequently for several applications ofincreasing complexity related to coal seam gas extraction and water management. The first application
concerns an evaluation of the natural soil attenuation potential of hydraulic fracturing chemicals andtheir transformation products, for which we use the standard HYDRUS model. The second application
involves an evaluation of the safe use of coal seam gas produced water for sustainable irrigation;
for this application we used UnsatChem. The third application analyzes major ion geochemistry andtrace metal transport in a soil leached with coal seam gas produced water; this application requires
the HP1 module. We note that these three one-dimensional applications are intended to only briefly
demonstrate the capabilities of the HYDRUS models and their add-on modules.
2. Theory
2.1. Variably-Saturated Water FlowHYDRUS-1D numerically solves a modified form of the Richards equation, which describes
one-dimensional uniform water movement in a partially-saturated rigid porous medium. The approachassumes that water flow due to thermal gradients can be neglected and that the air phase plays an
insignificant role in the liquid flow process:
 

∂
∂h
∂θ
=
K
+ cos α − S
∂t
∂x
∂x

(1)

where h is the soil water pressure head [L], θ is the volumetric water content [L3 L−3 ], t is time [T],x is the spatial coordinate [L], positive upward, K is the unsaturated hydraulic conductivity function[LT−1 ], S is a sink term (usually representing root water uptake) [L3 L−3 T−1 ], and α is the angle between
the flow direction and the vertical axis (e.g., 1 for vertical flow, 0 for horizontal flow). The unsaturated
soil hydraulic properties, θ(h) and K(h), in Equation (1) can be described in HYDRUS using theanalytical models of Brooks and Corey [34], van Genuchten [35], Vogel and Císlerová [36], Kosugi [37],or Durner [38]. A macroscopic approach is used to model root water uptake, which can be either
non-compensated or compensated and which can account for various stresses, such as the osmotic andpressure head stresses [39].
To account for preferential nonequilibrium flow, HYDRUS-1D can consider not only uniform
variably-saturated flow, but also flow in dual-porosity and dual-permeability systems ([40]).
Additionally, HYDRUS-1D can consider flow in both liquid and gaseous (vapor) phases, as well
as flow driven by both capillary and gravitational forces as well as by thermal gradients [41]. Since
these processes are beyond the scope of this paper, they will not be further discussed here.
2.2. Solute Transport Applications of the Standard HYDRUS Model
The standard HYDRUS models can simulate the fate and transport of individual solutes, or ofsolutes subject to relatively simple sequential first-order degradation reactions. The models, in general,
assume that solutes can exist in all three phases (liquid, solid, and gaseous) and that the decay andproduction processes can be different in each phase. Solutes must be defined in the liquid phase(the liquid phase concentration is the primary variable in the advection-dispersion equation given
below), but can be present additionally also in the solid and/or gaseous phases. Interactions between
the solid and liquid phases may be described by linear or nonlinear (Freundlich or Langmuir) typeequations, and can be either equilibrium or nonequilibrium. Interactions between the liquid andgaseous phases are assumed to be always linear and instantaneous. The models further assume that

Water 2017, 9, 385

4 of 33

the solutes are transported by convective transport and dispersion in the liquid phase, as well as bydiffusion in the gas phase. Hence, both adsorbed and volatile solutes, such as pesticides, fumigants
and certain metals, can be considered.
The partial differential equations governing one-dimensional nonequilibrium chemical transport
involved in a sequential first-order decay chain during transient water flow in a variably saturated
rigid porous medium are taken as:




∂θci
∂ρsi
∂av gi
∂
∂c
∂
∂qci
g ∂g
+
+
=
θDiw i +a v Di i −
+ Pi
∂t
∂t
∂t
∂x
∂x
∂x
∂x
∂x

(2)

where c, s, and g are solute concentrations in the liquid [ML−3 ], solid [MM−1 ], and gaseous [ML−3 ],phases, respectively; ρ is the soil bulk density [ML−3 ], av is the air content [L3 L−3 ], q is the volumetricflux [LT−1 ], Dw is the dispersion coefficient [L2 T−1 ] for the liquid phase, Dg is the diffusion coefficient[L2 T−1 ] for the gas phase, and P is a reaction/production term accounting for various biochemical
reactions [ML−3 T−1 ]. The subscript i represents the ith chain number. The reaction term Pi in
Equation (2) in the standard HYDRUS model is defined as follows:
P1 = −(µw,1 + µ0w,1 )θc1 − (µs,1 + µ0s,1 )ρs1 − (µ g,1 + µ0g,1 ) av g1 + γw,1 θ + γs,1 ρ + γg,1 av − r a,1
Pi = −(µw,i + µ0w,i )θci − (µs,i + µ0s,i )ρsi − (µ g,i + µ0g,i ) av gi + µ0w,i−1 θci−1 + µ0s,i−1 ρsi−1+µ0g,i−1 av gi−1 + γw,i θ + γs,i ρ + γg,i av − r a,i , i e (2, ns )(3)

where µw , µs , and µg are first-order rate constants for solutes in the liquid, solid, and gas phases [T−1 ],respectively; µw ', µs ', and µg ' are similar first-order rate constants providing connections between
individual species in the decay chain, γw , γs , and γg are zero-order rate constants for the liquid[ML−3 T−1 ], solid [T−1 ], and gas [ML−3 T−1 ] phases, respectively, and ra is the root chemical uptake
term [ML−3 T−1 ], which for passive uptake is equal to the product of the sink term S in the water flow
Equation (1) and the concentration of the sink term cr [ML−3 ]. The subscripts w, s, and g correspondwith the liquid, solid and gas phases, respectively, and ns is the number of solutes involved in the
chain reaction. The nine zero- and first-order rate constants in Equation (3) may be used to represent
a variety of reactions or transformations including biodegradation, volatilization, and precipitation.
When the constants µw ', µs ', and µg ' that provide connections between individual chain species are
all set to zero, Equations (2) and (3) may be used to simulate the fate and transport of individual solutes,
the behaviour of which is then assumed to be independent of that of other solutes that may be present
in the soil solution. When these constants are nonzero, then Equation (2) for i = 1 is used to simulate
the transport of a parent compound and for i > 1 the transport of daughter products. The above
formulation has proved useful in modelling a broad range of chemicals, including radionuclides(e.g., [15,42]), mineral nitrogen species (e.g., [16,43,44]), pesticides (e.g., [17]), chlorinated aliphatichydrocarbons (e.g., [18,45]), hormones (e.g., [19,46–49]), antibiotics (e.g., [20,50]), and explosives [51,52].
HYDRUS-1D at present considers up to ten solutes (five for the dual-permeability model), which either
can be coupled in a unidirectional chain or are allowed to move independently of each other.
Physical nonequilibrium solute transport is accounted for in later versions of HYDRUS byassuming a two-region, dual-porosity (or dual-permeability) type formulation to partition theliquid phase into mobile (fast) and immobile (slow) regions. Chemical nonequilibrium solutetransport is accounted for by assuming kinetic interactions between solutes in the liquid and solid
phases. Attachment/detachment theories, including filtration theory, were included to simulate
the transport of viruses, colloids, bacteria, nanoparticles, and/or nanotubes. Details about thesenonequilibrium models can be found in Šimůnek and van Genuchten [53]. The above standard
solute transport module of HYDRUS will be used below to demonstrate its application to evaluationsof the natural soil attenuation potential of hydraulic fracturing chemicals (e.g., biocides) and theirtransformation products.

Water 2017, 9, 385

5 of 33

2.3. Solute Transport Applications of the UnsatChem Module
The UnsatChem module represents a typical example of a model with specific chemistry [54].
Such models usually are constrained to very specific applications since they are restricted to very
prescribed chemical systems. On the other hand, these models can be numerically more efficient thanmodels with generalized chemistry (e.g., the HP1 module as discussed below), since the numericalsolution can then be optimized for a particular chemical system.
The UnsatChem geochemical module, which has been implemented in all 1D, 2D and 3D HYDRUS
versions, considers the transport and production of CO2 and the transport of several main components(i.e., Ca2+ , Mg2+ , Na+ , K+ , SO4 2− , CO3 2− , and Cl− ). CO2 concentrations are modelled since they cansignificantly affect soil pH and multiple chemical reactions, such as dissolution/precipitation reactionsof calcite. The Unsatchem module considers most or all relevant equilibrium and kinetic geochemical
reactions in such chemical systems as complexation, cation exchange, and precipitation-dissolution(e.g., of calcite, gypsum, and/or dolomite). A complete list of species in the code is given in theUnsatChem manual and in Table 3 of Šimůnek et al. [11]).
The carbon dioxide module of UnsatChem solves the convection-dispersion equation with asource term accounting for the production of carbon dioxide in soils. The major ion module ofUnsatChem solves Equation (2) for the seven main components, while the gas-phase related terms
and the general reaction term of Equation (2) are equal to zero. The second term on the left side ofEquation (2) is also zero for components that do not undergo ion exchange or precipitation/dissolution
reactions. While the aqueous complexation reactions do not affect Equation (2), cation exchange andprecipitation/ dissolution reactions move a component from the liquid phase into the solid phase andvice versa. The solution of the chemical system is described in detail by [23].
Many possible applications of the Unsatchem module exist since saline waters are often used
for irrigating agricultural crops in regions which either have limited water resources or an excess of
saline waters (e.g., as produced in mining operations). Irrigation with saline water may potentiallycause salinization and sodification of irrigated agricultural lands. Efficient irrigation and leaching
management practices are thus critical in these regions to prevent soil salinization when rainfall is
not sufficient to leach accumulated salts during or following irrigation. The UnsatChem module
has been used in many studies to evaluate the sustainability of alternative irrigation schemes withrespect to salinization and sodification processes, to assess reclamation of saline or sodic soils, and toevaluate the movement of salts after accidental releases or possible beneficial applications of saline
waters resulting from mining operations (e.g., [55–64], among others). We will use the UnsatChem
module later to demonstrate its potential for evaluating the use of coal seam gas produced water forsustainable irrigation.
2.4. Solute Transport Applications of the HP1 Module
The HP1 module represents a typical example of a model with generalized chemistry [54,65].
These models provide users with much more flexibility in designing particular chemical systems thanmodels with specialized chemistry and allow them to carry out a much broader range of applications.
Users then either can select species and reactions from large geochemical databases, are able to define
their own species with particular chemical properties and reactions, and/or can build their own
equilibrium and kinetic reaction networks.
The one-dimensional HP1 model [14,25] couples the PHREEQC geochemical program [24] with
HYDRUS-1D. Its two-dimensional extension, HP2, was released in 2013 as an add-on module to
HYDRUS (2D/3D) [11]. HPx, which is an acronym for HYDRUS-PHREEQC-xD (1D or 2D), is avery comprehensive simulation module for simulating (1) transient water flow, (2) the transport ofmultiple components, (3) mixed equilibrium/kinetic biogeochemical reactions, and (4) heat transport
in one- and two-dimensional variably-saturated porous media. The HP1 and HP2 modules are suitable
for a broad range of low-temperature biogeochemical reactions in water, the vadose zone and/or

Water 2017, 9, 385

6 of 33

ground water systems, including interactions with minerals, gases, exchangers and sorption surfaces
based on thermodynamic equilibrium, kinetic, or mixed equilibrium-kinetic reactions.
HP1 and HP2 allow thermodynamic equilibrium calculations for multiple chemical reactions andother features such as (1) aqueous speciation with different activity correction models (i.e., Davies,
extended Truesdell-Jones, B-Dot, Pitzer, and SIT - Specific Ion Interaction Theory); (2) multi-site ion
exchange sites with exchange described using different models (Gaines-Thomas, Vanselow, or Gapon);(3) multi-site surface complexation sites with non-electrostatic (e.g., the Dzombak and Morel or
CD_MUSIC models), as well as different options to calculate compositions of the diffuse double layer;(4) mineralogical assemblages; (5) solid-solution reactions; and (6) gas exchange. Kinetic calculations
can be used furthermore to describe mineral dissolution/precipitation reactions, non-equilibrium
sorption, biogeochemical reactions (including first-order degradation networks), Monod kinetics
and/or Michaelis-Menten kinetics [11].
The capabilities of HP1 were extended recently by considering diffusion of components (e.g., O2 or
CO2 ) in the gas phase. This is important since concentrations of gas components may significantly affectvarious chemical reactions or their rates [66]. Additionally, an option to change the hydraulic and solutetransport properties as a function of evolving geochemical state variables has been implemented [67].
For example, bacterial growth and/or clogging can affect porosity and corresponding physicalproperties. Similarly, precipitation/dissolution may lead to changes in porosity and corresponding
changes in the soil water retention and hydraulic conductivity functions. It is now possible to accountin HP1 for changes in porosity, hydraulic conductivity, tortuosity, dispersivity, and thermal properties
such as thermal capacity, conductivity and dispersivity. The flexibility of the embedded BASIC
interpreter permits HP1 users to define user-specific relationships between geochemical state variables
and the transport properties [68].
Several recent applications have illustrated the versatility of HP1 to simulate a variety of complex
reaction networks. Early applications [14,25] simulated long-term transient flow and the transport ofmajor cations and heavy metals subject to a multi-site pH-dependent cation exchange complex and/or
surface complexation in a multi-layered soil profile. A complex reaction network of mercury, which
included kinetic dissolution from a mercury-containing solid or from an immobile non-aqueous liquid,
as well as redox reactions and exchange with the soil air phase, implemented in HP1, allowed theidentification of key processes and parameters for mercury transport in soils [69,70]). The network wasapplied to the treatment of mercury-contaminated soils with activated carbon [71]. CO2 production
and transport in bare and planted mesocosms [72], and the effects of lime and concrete waste on vadosezone carbon cycling [73] furthermore illustrate applications to the problem of CO2 sequestration insoils. Recent applications are also in the field of colloid and colloid-affected transport, including
competitive kinetic sorption on colloids [74] and accounting for the effects of water content and ionicstrength on attachment at the air-water and soil-water interfaces, respectively [75]. The HP1 module
will be used below to demonstrate its application to evaluating the fate and transport of major ionsand trace metals in soils leached with coal seam gas produced water.
3. Applications
3.1. Standard HYDRUS-1D: Fate of Hydraulic Fracturing Chemicals During Natural Attenuation
3.1.1. Problem Definition
Bronopol is one of many biocides used to prevent bacteria growth in the subsurface. Biocides are
essential components of hydraulic fracturing fluids used for unconventional gas extraction [76–78].
Bronopol also serves as a preservative in cosmetic products, liquid soaps, and cleaning agents [79].
Within well environments, bacteria often cause bioclogging and inhibit gas extraction, produce toxic
hydrogen sulfide, and induce corrosion, which may lead to downhole equipment failure. Within thecontext of unconventional gas extraction, [77] found that uncharged biocides will dominate in the
aqueous phase and be subject to degradation and transport, whereas charged biocides will sorb toWater 2017, 9, 385

7 of 33

solids and be less bioavailable. They also note that many biocides are short-lived or degradable throughabiotic and biotic processes, while some may transform into more toxic or persistent compounds.
The aquatic toxicity of bronopol, often expressed in terms of the half-maximal inhibitory
Water 2017, 9, 385 EC50 , or half-lethal concentration, LC50 , has been reported for many marine 7andof 34
concentration
freshwater biota, as well as for aquatic birds ([80]). Microalgae are the most sensitive to bronopol−1(EC50 == 0.020–0.41
by oysters(EC50(EC= 0.42=and
0.78
mg0.78L−1),mgwater(EC50fleas
= 1.6
(EC
0.020–0.41mg
mgL L),−1followed), followedby oysters
0.42
and
L−1fleas), water
5050−1), mysids (LC50 = 4.3 and 5.9 mg L−1), duckweed (EC50 = 38 mg L−1) and fish (LC50 = 7.5–59 mg
mg
L
−
1
−
1
−
1
(EC50 = 1.6 mg L ), mysids (LC50 = 4.3 and 5.9 mg L ), duckweed (EC50 = 38 mg L ) and fishL−1). Becauseof its low−1 octonal/water ratio and high solubility in water, bronopol is not expected to(LC
50 = 7.5–59 mg L ). Because of its low octonal/water ratio and high solubility in water, bronopol
bioaccumulate.
RapidhydrolysisRapid
occurshydrolysis
under warm
and/orhigh
pHand/or
conditions,in theis not expected to bioaccumulate.
occursunder
warmhigh resulting
pH conditions,formationofformaldehyde(CH2O) which should also be considered in environmental hazard
resulting in the formation of formaldehyde (CH2 O) which should also be considered in environmental
assessments[81]. [81].
hazardassessments
Studies
about human
to tobronopolarearenotnotwidelyavailable,with with
most
Studies abouthuman healthhealtheffects
effectsofofexposure
exposure
bronopolwidelyavailable,
or
all
data
on
bronopoltoxicity
being
limitedtostudiesonanimals.
Based
on
these
studies,estimatedmost or all data on bronopol toxicity being limited to studies on animals. Based on these studies,
occupational
short- andshortintermediate-term
exposure levels
andlevels
risks toworkers
beenhavecalculated
estimated
occupational
and intermediate-term
exposureandrisks tohave
workers
been
along
with
its
classification
as
a
Group
E
chemical(i.e.,
one
for
which
there
is
evidenceof noncalculated along with its classification as a Group E chemical (i.e., one for which there is evidenceofcarcinogenicity
tohumans;[81]).
non-carcinogenicity to humans [81]).
Substantial spills
spills ofof bronopol
bronopol into
into surface
surface waters
waters oror streams
streams bybyindustries
industries such
such asasthosethosethat
that
Substantial
involvehydrocarbon
extraction
maytherefore
have
noticeable
ecotoxicologicaleffectsonaquatic
involve hydrocarbon extraction may therefore have noticeable ecotoxicological effects on aquatic
species. These
These effects
effects will
willdepend,depend,among
amongother
otherfactors,
factors,ononbronopol’s
bronopol’spotential
potentialfor
fordegradation
degradationininspecies.
naturalenvironments.
natural environments.

3.1.2.Biocide
BiocideTransformation
TransformationPathways
Pathwaysand
andNatural
NaturalSoil
SoilAttenuation
AttenuationParameters
Parameters3.1.2.
Bronopol(BNP)(BNP)isisrelatively
relativelystablestableunder
undernormal/ambient
normal/ambient environmental
environmental conditions,conditions,but
but veryvery
Bronopolsusceptible toto degradation
degradation (hydrolysis)(hydrolysis) when
when exposed
exposed toto elevatedelevated temperaturestemperatures andand alkaline
alkalinesusceptible◦
conditions
[81].
For
example,bronopolhas
beenreportedtohydrolyze
within
3
h
at
60
°C
and
conditions [81]. For example, bronopol has been reported to hydrolyze within 3 h at 60 CpH
and8,
producingformaldehyde,
nitrosamines,
2-hydroxymethyl-2-nitropropane-1,3-diol(tris)
and2pH 8, producing formaldehyde, nitrosamines, 2-hydroxymethyl-2-nitropropane-1,3-diol (tris) andbromo-2-nitroethanol [81–85].
reported also,
also, including
including
2-bromo-2-nitroethanol[81–85]. Other
Otherdegradation
degradation products
products have
have been
been reported
nitromethane
and
the
more
persistentbromonitromethane[82,86,87].
Figure1shows
someof thethenitromethane and the more persistent bromonitromethane [82,86,87]. Figure 1 shows some ofdegradation
pathwaysbronopolaccording
toFor
[87].
the purposeof thisthestudy,
the degradation
degradation
pathwaysforforbronopolaccordingto [87].
theForpurposeof this study,
degradation
pathway
pathway
will betosimplified
to only threei.e.,
components,
i.e.,(BNP)
fromtobronopol(BNP) to 2-bromo-2will
be simplified
only three components,from bronopol
2-bromo-2-nitroethanol(BNE),
nitroethanol(BNE),
andthentobromonitromethane(BNM).
Importantly,
theecotoxicity
changesand then to bromonitromethane (BNM). Importantly, the ecotoxicity changes along the degradation
along theas
degradation
pathway,
as is2exemplified
Figure 2pyrenoidosa
for the algae
Chlorellapyrenoidosa[87].
pathway,
is exemplified
in Figurefor the algaeinChlorella[87].
EC50 values
decrease
EC
50valuesdecrease
from
BNP
to
BNE
to
BNM,
with
BNE
and
BNM
being
about
twice
and
five
times
from BNP to BNE to BNM, with BNE and BNM being about twice and five times more toxic as BNP at
more
as BNP
at 48 hrs exposure time.
48
hrs toxic
exposure
time.
OH

HCHO

HCHO

Br

Br
Br
HO

N+
O

O-

2-Bromo-2-nitro-1,3-propanediol (BNP)

HO

N+

N+
O

O-

2-Bromo-2-nitroethanol (BNE)OO-

Bromonitromethane (BNM)

FigureSimplified
degradation
pathway
for bronopol(2-bromo-2-nitro-1,3-propanediol)
as implemented
Figure1. 1.
Simplified
degradation
pathway
for bronopol(2-bromo-2-nitro-1,3-propanediol)asin
HYDRUS(modifiedfrom[87]).
implemented in HYDRUS (modified from [87]).

Water 2017, 9, 385

8 of 33

Figure 2. Algae (Chlorella pyrenoidosa) ecotoxicity for the BNP=>BNE=>BNM degradation pathway as
a function of exposure time (data from [87]).

The key parameters determining biogeochemical attenuation in the subsurface are the soil-waterpartitioning coefficient KD [L3M−1] and the first-order degradation constant µ = ln2/T1/2 (see Equation (3)),
where T1/2 is the chemical’s half-life [T]. For nonpolar organic chemicals, KD is strongly related to theorganic carbon in soil [88]. Organic compounds may also be sorbed onto mineral surfaces, especiallywhen dry. When soils are moist, however, water will generally displace organic compounds from themineral sorption sites; organic matter sorption then becomes more important than the small remaining
contribution of mineral sorption [89]. Polar and ionizable organic compounds may sorb on both themineral soil fractions and organic matter, with soil water content and pH often playing an importantrole [90].
1
For organic chemicals, KD is strongly related to the soil organic carbon content, f OC (expressed as
g carbon g−1 soil); it is common practice for organic chemicals to derive KD from the organic carbon
partition coefficient KOC , i.e., KD = f OC × KOC . The same dependency is true for the degradation
constant µ. Generally, the higher the organic carbon content, the higher is KD and the faster isbiodegradation. A higher organic carbon content favors microbiological activity, and as such also
biodegradation. Organic carbon is generally highest in the soil top layers and decreases with depth.
This implies that KD and µ generally decrease also with soil depth. For the polar chemicals BNP, BNE,
and BNM sorption is typically described through use of the KOC parameter.
For deriving soil specific liquid-solid partitioning coefficients KD , soil organic carbon data typicalof three soil types (Vertosols, Chromosols, and Sodosols) in New South Wales, Australia, were obtained
from [91]. Organic carbon data decreased exponentially with depth, with a geometric mean fraction(f OC ) of 0.0217 in the top 10 cm and a minimum background value of approximately 0.0014 at a depthof 150 cm. From there downwards, f OC remained more or less constant with depth. The exponentiallydecreasing organic carbon content with depth was approximated in HYDRUS-1D using a linearized
depth-dependent organic carbon function. A summary of organic carbon partitioning coefficients KOC
and aerobic degradation half-life values T1/2 obtained from a literature review is provided in Figure 3(for further details see Supplementary Materials Table S1).

Water 2017, 9, 385

9 of 33

Figure 3. Soil organic carbon partition coefficient (KOC ) and half-life (T1/2 ) for the BNP=>BNE=>BNM
degradation pathway.

The average value of f OC in the top 10 cm (0.0217 g carbon g−1 soil) was multiplied byliterature-based KOC values of the three compounds (BNP, BNE and BNM) to obtain their KD values(i.e., KD = f OC × KOC ). This procedure was repeated for the other depths, leading to depth-dependent
KD profiles for bronopol, 2-bromo-2-nitroethanol, and bromonitromethane. Finally, HYDRUS internallyconverted KD values into retardation factors Ri , given by:
Ri = 1 +ρK D,iθ(4)where ρ is dry bulk density [ML−3 ], θ is soil water content [L3 L−3 ] and KD,i the liquid-solid partitioning
coefficient for linear reversible sorption [L3 M−1 ]. Since in our case none of the compounds were volatile,and hence did not partition into the air phase, the solute transport equations given by Equation (2)simplify to


∂
∂qc1
∂θR1 c1
w ∂c1
=
θD1
−
− µ10 θc1
(5)
∂t
∂x
∂x
∂x


∂θRi ci
∂
∂qci
w ∂ci
=
θDi
−
− µi0 θci + µi0−1 θci−1
i = (2, 3)
(6)
∂t
∂x
∂x
∂x
which assumes that degradation occurred only in the liquid phase.
For the solute transport simulations we assumed a 100-cm-deep soil profile with three soil layers:
0–30 cm, 30–70 cm, and 70–100 cm. The spatial discretisation was 1 cm; the dispersion coefficient
was calculated based on a dispersivity of 10 cm or one tenth of the domain length, while molecular
diffusion was neglected. Bulk densities for the three soil layers were 1.32, 1.42 and 1.45 g cm−3 forlayers 1, 2 and 3, respectively. As an example, we used a θ value of 0.35 cm3 cm−3 and bulk density of
1.32 g cm−3 to calculate R values shown in Table 1 (representative for the top soil layer). Since the watercontent θ will vary with time and depth, R values are also time and depth-dependent. A summary of
the solute transport parameters R and aerobic degradation half-life values T1/2 is provided in Table 1(top soil layer values shown as example). The R and T1/2 values shown in Table 1 are based on aliterature review (see Supplementary Materials Table S1) and represent values for the top 0–10 cm
soil layer using the organic carbon value (f OC ) of 0.0217 g g−1 as a reference. Both R and µ decreased
linearly with depth according to the decrease in f OC . For this purpose, a new module was developed in
HYDRUS-1D to automatically update KD and the degradation constant µ as organic carbon decreases
with depth. The retardation and degradation constants are reduced by a depth-dependent function,5

Water 2017, 9, 385

10 of 33

which has a constant value of one (depth-reduction factor = 1) in the top of the soil profile (between
soil surface and depth = A), then decreases linearly down to the value C at depth B. Below the depth B
the reduction function has the value of C (Figure 4). The value for C is calculated as the ratio of organiccarbon at depth B (typically the depth where organic carbon becomes more or less constant) to theorganic carbon at depth A (i.e., the mean value between soil surface and depth A). In the hypotheticalexample of Figure 4, the mean organic carbon % from soil surface to depth A is 0.65, and the carbon
content at depth B is 0.12. As result, C = 0.12/0.65 = 0.185.
Table 1. Retardation factors R and aerobic degradation half-lives (T1/2 , days) of bronopol and two of
its daughter productions. N = number of data (see Supplementary Materials Table S1 for data sources).
ChemicalMean R

Min R

Max R

N

Mean T 1/2

Min T 1/2

Max T 1/2

N

Bronopol (BNP)
2-bromo-2-nitroethanol (BNE)
Bromonitromethane (BNM)

187
2
17

2
2
14

1170
2
18

11
1
3

87
28
18

0.1
11
8.7

548
51
30

7
5
3

Figure 4. Hypothetical example demonstrating organic carbon profile (a), the derived depth-reduction
factor (b), the depth-reduction factor applied to the distribution coefficient for three chemicals involved
in first-order transformation reactions (c) and the first-order degradation constant for the same threechemicals (d).

Both mean retardation and half-life decreased within the degradation pathway (Table 1): thedegradation products of BNP become more mobile (mean R values decreased from 187 to 2 and 17),and they degrade faster (mean T1/2 values decreased from 87 to 28 and 18) compared to BNP.
Simulations started with a 20 year infiltration period with solute-free water to achieve steady-stateflow conditions. Next, simulations of the leaching behavior of BNP and its degradation products
were carried out for a hypothetical leak of produced water from a water holding pond at a rateof 3.5 cm y−1 with a hypothetical concentration of 1 mg L−1 during a 30-year period (the average
lifespan of a gas project during which water holding ponds are present). After this period, solute-freewater was applied at a rate of 4 cm y−1 reflecting long-term mean recharge. For the unsaturated soil
hydraulic properties we used the functions of van Genuchten [35] based on measured soil data points
within the Cox’s Creek catchment, New South Wales (Supplementary Materials Table S2). Vertosols
are the predominant soil type observed in 84 out of 143 soil profiles (59%) ([92]). Breakthroughcurves for BNP, BNE, and BNM were produced at different depths to gain insight in the coupledprocesses of advective-dispersive transport, sorption, first-order degradation, and the generationof transformation products (notably BNE and BNM). The focus will be on the shallow soil depthsof 2 and 5 cm as deeper in the profile concentrations became extremely small. Three combinations

Water 2017, 9, 385

11 of 33of sorption and degradation parameters were considered using values from Table 1: {max R, min
T1/2 }, {min R, max T1/2 }, and {mean R, mean T1/2 }. The first parameter combination with a maximum
retardation and minimum half-life should result in maximum chemical attenuation of BNP as discussed
in the next section. The second parameter combination uses the minimum retardation and maximum
half-life, thus causing minimal attenuation for BNP. The third parameter set was expected to generate
an average attenuation for BNP. The degree of attenuation of the transformation product BNE should
not only depend on the value of its own parameter combination, but also on those of BNP. For BNM,
attenuation depends on parameter combinations from all three chemicals.
3.1.3. HYDRUS-1D Simulation Results for the Bronopol Degradation Chain
Calculated breakthrough curves of BNP and its degradation products BNE and BNM are shownin Figure 5. The plots reveal much uncertainty in BNP concentrations owing to the large range in
R and T1/2 values. The range in concentrations is much smaller for BNE and BNM owing to thesmaller range of their parameters R and T1/2 . For BNP, the maximum and minimum concentrations
were indeed determined by the parameter combinations {min R, max T1/2 }_BNP and {max R, min
T1/2 }_BNP, respectively (Figure 5, top). For the first degradation product of BNP (i.e., BNE), themaximum breakthrough concentrations is determined by the parameter combination {mean R, mean
T1/2 }_BNE (Figure 5, middle plots). This is a result of the parameter combinations of its parent BNP,
where {mean R, mean T1/2 }_BNP generates more BNE per unit of time than {min R, max T1/2 }_BNP.
Moreover, the latter parameters have a very high R value of 1170, thus causing BNP to move muchslower through the soil than for the former parameters with R = 187. For BNM the sequence of high,mean and low breakthrough concentrations and their respective parameter combinations is the same
as for BNE at 2 cm depth (Figure 5, bottom). For BNM at 5 cm depth, the {min R, max T1/2 }_BNM
parameter combination produced the lowest concentrations, while the {max R, min T1/2 }_BNM
parameter combination yielded the mean breakthrough curve concentrations (Figure 5, bottom). Thisis due to the considerably smaller T1/2 value (8.7 days) for the latter set compared to 30 days for theformer set, while the R values are very similar, leading to much faster BNM degradation.
Overall the calculations demonstrate that the degradation products have consistently smaller
concentrations than the parent compound, with maximum concentrations that are up to 100 times
smaller than the parent source concentrations at depths of 2 to 5 cm. This demonstrates thesignificant natural attenuation capacity of soils as a result of sorption and biodegradation. The riskof bioaccumulation in soil or leaching to groundwater is hence very small for the conditions of thisstudy, i.e., for a small infiltration flux approximately equal to the long-term recharge rate of several
tens of mm per year at the study area. We note that the methodology developed here can be easily
extended to a more systematic screening of the leaching risks of a broad range of organic compounds
and different soil types.

Water 2017, 9, 385
Water 2017, 9, 385

12 of 33
12 of 34

Figure 5. Breakthrough curves of BNP (top), BNE (middle), and BNM (bottom) at 2 and 5 cm soil

Figure 5. Breakthrough curves of BNP (top), BNE (middle), and BNM (bottom) at 2 and 5 cm soil
depth using three parameter combinations: min or low R and max or high half-life (L), mean R and
depth
using(M),threeparameter
combinations:min or(H).
low R and max or high half-life (L), mean R andhalf-lifemax
or high R and
min or low half-life
half-life (M), max or high R and min or low half-life (H).

Overall the calculations demonstrate that the degradation products have consistently smaller
concentrations thanthe parent
compound,maximumthatare up to 100 times
3.2. UNSATCHEM:
Optimizing
Coal
Seam Gaswith
Produced
Waterconcentrationsfor Sustainable
Irrigation
smaller than the parent source concentrations at depths of 2 to 5 cm. This demonstrates the significant
attenuation
capacity of soils as a result of sorption and biodegradation. The risk of
3.2.1.natural
ProblemDefinition

The second example considers coupled processes of variably-saturated water flow, plant water
uptake and the simultaneous transport of multiple major ions in soils irrigated with produced coalseam gas water featuring different water qualities. By coupling major ion soil chemistry to unsaturated

Water 2017, 9, 385

13 of 33

flow and plant water uptake, and by explicitly incorporating the effects of salt concentrations on
soil hydraulic properties and root water uptake (salinity stress), critical soil processes required fora salinity risk assessment associated with coal seam gas produced water can be included in theanalysis. Simulations with different irrigation water qualities provide detailed insight regarding
chemical indicators of soil and plant health, i.e. Sodium Adsorption Ratio (SAR), EC and sodium
concentrations. We compare such indicators in the soil profile with permissible soil quality guideline
values in Australia (ANZECC values—[93]) to assess the risk to soil and plant health. Insights frommodelling may be used to provide a scientific basis for defining safe water quality requirements forsustainable irrigation with coal seam gas produced water.
Produced water from unconventional gas extraction is generally unsuitable for direct surface
discharge or irrigation without any treatment or amendment [94–96]. While high sodium
concentrations cause soil particles to disperse, particularly if montmorillonite clays are present, mostions increase the aggregation of soil particles [96]. Irrigation water with a high SAR can therefore
lead to a decrease in infiltration and deterioration of the soil structure as the dispersed clay minerals,
once dry, cause soils to become dense, cloddy and structureless, thereby destroying natural particle
aggregation. For example, SAR values between 5 and 8 have been shown to cause irreversible plugging
of soil pores and swelling [97]. Furthermore, SAR values greater than 13 pose a risk to the soil
ecosystem [94]. Typical treatments to mitigate the effects of saline-sodic irrigation water include theaddition of gypsum and elemental sulphur [21,98].
The above problems can be addressed with HYDRUS using three approaches. One approach
would be to use the standard HYDRUS models (similar to the bronopol example) by assumingthat salinity behaves more or less like an inert tracer and hence is not subject to chemical reactions(e.g., [43]). A second approach is to use the UnsatChem module (the current example), which considers
the transport and reactions between major ions (e.g., [44,56]). While the former approach does not
permit such processes as cation exchange, dissolution of mineral amendments (e.g., gypsum or
calcite) or precipitation of these minerals when the soil solution becomes oversaturated, the latterapproach allows one to consider those geochemical processes and the effects of salts and soil water
quality on soil properties and plant water uptake. Finally, a third approach would use the coupled
HYDRUS-PHREEQC code [14,25], which provides users with an even greater suit of dissolved chemical
species, mineral phases, and biogeochemical reactions relevant to multi-component transport (see thethird example of this paper).
3.2.2. Soil Hydraulic and Chemical Model for Major Ion Transport
Three different chemical compositions of irrigation water were assumed (Table 2): (1) untreatedproduced water, (2) 3:1 produced water amended with surface water, and (3) produced water treatedwith reverse osmosis (RO). The untreated and treated compositions represent end-members that canbe mixed in different ratios to obtain fit-for-purpose irrigation water. Unlike the amended and treatedwaters, which have actually been field-tested and/or have been used in commercial irrigation schemes,
the untreated water would almost certainly not receive any permission for use as irrigation waterowing to its high potassium and SAR values (Table 2). The untreated water is included in this study todemonstrate the effects of several interacting processes and identification of tipping points, such as
what soil solution composition may trigger a reduction in the soil hydraulic conductivity, which thenwould modify water redistribution in the soil profile and impact on root water uptake. We note thatthe treated water is still much more saline than rainfall (e.g., Na concentrations are 60 times higher
than in rainfall, Table 2). Irrigation with RO treated produced water hence is likely to have some minor
effects on the soil’s chemical balance.

Water 2017, 9, 385

14 of 33

Table 2. Composition of untreated, amended, and reverse osmosis (RO) treated produced watersfrom Australian coal basins. Treated/rain water = a ratio of treated produced water to rain watercomposition * [99]. SAR = sodium adsorption ratio; ECT = total electrical conductivity.
Units

Alkalinity

Cl

SO4

Ca

K

Mg

Na

ECT

SAR

1406
61.1
3064
53.7

5705
100%

77
-

Untreated produced watermg L−1
meq L−1
ECi (µS cm−1 ) @%ECi @

1706
41.9
1245
21.8

593
16.7
1278
22.4

22.8
0.47
38
0.7

10.4
0.52
31
0.5

5.6
0.14
11
0.2

8.9
0.74
39
0.7

Amended produced water (3:1 ratio surface water:produced water)
mg L−1
meq L−1
ECi (µS cm−1 ) @%ECi @

183
4.5
134
8.5

164
4.6
353
22.5

337
7.0
561
35.7

10.3
0.51
31
1.9

34.8
0.89
65
4.2

6.5
0.53
28
1.8

184
8.0
401
25.5

1573
100%

11
-

291
100%

1.6
-

RO-Treated produced watermg L−1
meq L−1
ECi (µS cm−1 ) @%ECi @
Treated/rainwater*

25
0.61
18.2
6.3

67
1.89
144.3
49.6

1
0.02
1.67
0.6

10.33
0.52
30.7
10.5

0.43
0.011
0.81
0.3

7.6
0.62
33
11.3

28.7
1.25
62.5
21.5

417

103

3

43

4

95

60

Molar conductivity Λi for ions as trace concentration in water at 25 ◦ C
S cm2 mol−1

44.5

76.35

160

119

73.5

106

50.1

@ Calculated with Equation (8).

The same 100-cm-deep soil profile is used as in the previous example, i.e., having three soil
horizons (0–30 cm, 30–70 cm, and 70–100 cm), the same van Genuchten soil hydraulic properties, andthe same dispersion coefficient. Initial concentrations of the soil exchangeable cations were takenfrom a lysimeter located near Narrabri in northern New South Wales [100] of which the soil wascharacterised as a Haplic Grey Vertosol. The bulk density ranged from 1.32 to 1.45 g cm−3 for different
soil depths [100]. The cation exchange capacity and the concentration of exchangeable cations of thefirst layer were assigned to all three horizons. The initial soil solution composition was defined byusing the chemical composition of rainwater [99], which in the UnsatChem module is automaticallyequilibrated with the cation exchange complex.
The following hypothetical irrigation scheme was used for evaluating the potential build-up ofmajor ions for three different irrigation water qualities: over each 5-day period, 6 cm of irrigation wasapplied on the first day followed by 4 days without irrigation; this resulted in a total of 73 irrigation
days out of a simulation period of 365 days. Potential evapotranspiration (ET) was assumed to be
1 cm day−1 (based on the average from January 2014 in Narrabri, Australia [101]), and divided into
potential evaporation (0.2 cm day−1 ) and potential transpiration (0.8 cm day−1 ). For each 5-day period,
this yielded an irrigation excess of 20% above the potential ET for each day.
In our calculations we did not consider any rainfall. While this is a simplification of the soilhydrological balance and neglects the effects of rainfall salinity or salt leaching on the salt mass balance,
it does allow us to focus on the effects of irrigation only using coupled chemical-hydrological-biologicalprocesses. A 50-cm-deep root zone for pasture was assumed, with the root distribution decreasing
linearly with depth. Soil water and chemical mass balances were calculated for one year. Unsaturated
water flow was solved using Equation (1) and accounted for water uptake by grass. The volume
of water removed from a unit volume of soil per unit time due to plant water uptake, S(h) [day−1 ],is defined as [11]:S(h) = α(h)b(z) Tp(7)

Water 2017, 9, 385

15 of 33

where α(h) is the plant-water stress response function (between 0 and 1), b(z) is the normalized
water uptake distribution [cm−1 ] (here assumed linearly decreasing with depth), and Tp is potential
transpiration [cm day−1 ]. At the soil surface b = 0.04 cm−1 and S(h) = 0.032 day−1 .
Four modelling scenarios were considered using the Unsatchem major ion chemistry module of
HYDRUS to calculate concentrations of major ions (Ca, Mg, Na, K, alkalinity, Cl, SO4 ), SAR, and EC in
the soil profile. We used for these calculations the chemical compositions of the untreated, amended
and RO-treated produced waters from Australian coal basins shown in Table 2. The first scenario uses
untreated produced water, while the effects of solution composition on soil hydraulic conductivity
are ignored. This scenario allows us to demonstrate the errors made when using models that do notaccount for solution-driven impacts on hydraulic conductivity. In Table 2 the electrical conductivity
contribution of each ion (%ECi ) to the total EC of the solution was calculated as ([102]):%ECi =

100 × ECi
; ECi = Λi × mi
∑ ECi

(8)

where Λi is molar conductivity (S cm2 mol−1 ) and mi is concentration (mol L−1 ). Note this isapproximate as Equation (8) is valid for ions at trace concentrations.
For the second scenario we used untreated produced water at a lower 4 cm day−1 irrigation rate to
avoid water ponding due to reduced hydraulic conductivity. This in contrast to the initial calculations
for which the 6 cm day−1 irrigation rate did not lead to ponding (results not shown). This second
scenario and all further scenarios account for possible reductions in the hydraulic conductivity due tothe solution composition, based on the approach of McNeal [103,104]. The third scenario uses surface
water-amended produced water (at a 3:1 surface water: produced water ratio), again for a 6 cm day−1
irrigation rate. This last scenario considers RO treated produced water using the chemical composition
shown in Table 2. The second, third and fourth scenarios all involved the same coupled processes,
such as water balance calculations that account for root water uptake along with multiplicative waterand solute stresses.
3.2.3. UNSATCHEM Simulation Results
When untreated water was used for irrigation (Scenario 2), the simulated SAR and EC evolutions
showed high SAR values of over 100 (Figure 6a) and EC values exceeding 10 dS m−1 by the end of
the simulation period (Figure 7a). The ANZECC SAR guide value of 20 for sand-sandy loam soil was
exceeded across the entire soil profile, from the very start of the irrigation. The EC ANZECC trigger
value of 0.65 dS m−1 was also exceeded quickly after irrigation started, i.e., from the surface down
to 50 cm after 60 days. The rate at which SAR moves through soil is determined by the transport ofNa, Ca, and Mg. The electrical conductivity (Figure 7a) progressed much faster than SAR (Figure 6a)
because it is mainly determined by mono-charged sodium (54%) and the non-reactive and fastermoving chloride (22%) and alkalinity (22%) (Table 2). The EC reached nearly 6 dS m−1 at the bottom of
the profile after 365 days (Figure 7a), while after 365 days the SAR front (SAR up to 100) was between
30 and 40 cm depth (Figure 6a).
Using amended water (Scenario 3), the resulting SAR profiles shown in Figure 6b demonstratethat the SAR ANZECC trigger value (SAR = 20) for sand-sandy loam was exceeded after 365 days,
but only in a small section of the profile. The SAR ANZECC trigger value (SAR = 5) for light clay was
exceeded in the top 30 cm after 150 days, whereas after 365 days it was exceeded in the top 50 cm.
The EC ANZECC trigger value of 0.65 dS m−1 was exceeded in the entire soil profile after 90 days
(Figure 7b). The relatively high EC was due to the (still relatively high) EC of the irrigation water(1.57 dS m−1 or 1570 µS cm−1 , Table 2). In terms of contributions of individual ions in the amended
irrigation water to the EC, Table 2 shows that the most important contributors to EC were sulfate (a 36%
contribution to total EC), sodium (a 26% contribution to total EC) and chloride (a 23% contribution).
Optimization of irrigation water quality can be done by reducing concentrations of those ions that

Water 2017, 9, 385

16 of 33

contribute most to the total EC. In the above example the decreasing order of importance was sulfate,sodium,and chloride.
Water 2017, 9, 385
16 of 34

(A)

0

60 d

(B)

150 d

-20

365 d

-40

-60

ANZECC
guide value
for sand-sandy loam

-80
ANZECC
guide value
for light clay

-100
0

5

scenario 3

10

15

20

25

SAR

(C)

Figure 6. Simulated SAR values for Scenario 2 (A), Scenario 3 (B) and Scenario 4 (C). The ANZECC

Figure 6. Simulated SAR values for Scenario 2 (A), Scenario 3 (B) and Scenario 4 (C). The ANZECC [93][93] guide value for SAR is 5 for light clay and 20 for sand-sandy loam soil. Note the different scales
guide value for SAR is 5 for light clay and 20 for sand-sandy loam soil. Note the different scales inin the figures.
the figures.

Using amended water (Scenario 3), the resulting SAR profiles shown in Figure 6b demonstratethat theSAR ANZECC
triggervalue (Scenario(SAR = 20)4)
forthesand-sandy
loam
exceeded
after 365
days, theWhen
RO-treatedwater
was usedsimulated
SARwasprofilesremained
below
only
in a small trigger
section of
the profile.
SAR
ANZECCvalue(SARwas= 5) forlight clay
was thelowestbut
SAR
ANZECCvalue
of 5 forThe
light
clay
(Figuretrigger
6c). This
resultexpected
given
exceeded
in
the
top
30
cm
after
150
days,
whereas
after
365
days
it
was
exceeded
in
the
top
50
cm.soillow SAR value (i.e., 1.6) of the irrigation water and the absence of any initial salt build-up in the−1 was exceeded in the entire soil profile after 90 days
The
EC
ANZECC
trigger
value
of
0.65
dS
m
profile (the initial SAR value was 1.25 throughout the entire soil profile) or without other salt inputs(Figure 7b). The relatively high EC was due to the (still relatively high) EC of the irrigation water

Water 2017, 9, 385

17 of 33

into the profile (e.g. from saline shallow groundwater or lateral inflow from water discharging at the
break-of-slope [105]). The EC ANZECC trigger value of 0.65 dS m−1 was not exceeded across theentire soil profile, except at the very top few centimeters (Figure 7c). Although the EC of the irrigationwater was low (0.29 dS m−1 or 290 µS cm−1 , Table 2), the EC in the top of the soil profile increased
to a maximum value of 0.73 dS m−1 due to evapotranspiration, leading to dry soil conditions andtriggering salts to concentrate near the soil surface. In terms of contributions of individual ions in thetreated irrigation water to the EC, the most important contributors were chloride (a 49.6% contribution
to total
EC), sodium (21.5%) and magnesium (11.3%), as shown in Table 2.
Water 2017, 9, 385
18 of 34
(A)

(B)

(C)

Figure 7. Simulated electrical conductivity for Scenario 2 (A), 3 (B) and 4 (C). The ANZECC [93] EC

Figure 7. Simulated electricalconductivity for Scenario 2 (A), 3 (B) and 4 (C). The ANZECC [93] EC
trigger value is 0.65 dS m−1. Note the different scales in the figures.
trigger value is 0.65 dS m−1 . Note the different scales in the figures.

Cumulative actual transpiration is a useful metric to evaluate the combined effects of soil
chemical
physicalprocesseson theofsoil
water balance
andwaterplant qualities
water stress.
Thesoilcumulative
The nextand
sections
discuss
the effectdifferent
irrigationon thewater balance
actual transpirationof all scenarios
was less than(actual
its potential
value (292 cm);
4,componentsof evaporation,
plant transpirationand potential)
and fornetScenarios
drainage.3 andCumulativetranspiration reached 281 cm or 96% of the potential value (Figure 8b). The plants were found to
experience a very little water stress despite being irrigated with amended and treated produced water,which is most likely due to the sub-optimal irrigation regime generating soil water pressure heads
slightly outside the optimal range, especially near the more negative soil pressure heads where thelimiting pressure head becomes exceeded. Larger differences between actual and potential
transpiration occurred for Scenarios 1 (82% of the potential value) and 2 (64% of potential value). For

Water 2017, 9, 385

18 of 33

actual evaporation (Figure 8a) reached its maximum value of 73 cm (based on a 0.2 cm day−1 potential
evaporation rate) in all scenarios except Scenario 3. The latter scenario yielded 69 cm or 95% of thepotential value owing to the reduced infiltration, thereby generating drier soil moisture conditions andhence less evaporation.
Cumulative actual transpiration is a useful metric to evaluate the combined effects of soil
chemical and physical processes on the soil water balance and plant water stress. The cumulative
actual transpiration of all scenarios was less than its potential value (292 cm); for Scenarios 3 and4, transpiration reached 281 cm or 96% of the potential value (Figure 8b). The plants were found
to experience a very little water stress despite being irrigated with amended and treated produced
water, which is most likely due to the sub-optimal irrigation regime generating soil water pressure
heads slightly outside the optimal range, especially near the more negative soil pressure heads
where the limiting pressure head becomes exceeded. Larger differences between actual and potential
transpiration occurred for Scenarios 1 (82% of the potential value) and 2 (64% of potential value). For
Scenario 1 we used the combined effects of water stress and solute stress—in combination with the
major ion chemistry module—to control root water uptake. This modeling approach caused additional
stress on plant water uptake, mainly due to salinity stress (the irrigation water EC was 5700 µS cm−1 ).
The lower actual transpiration rate for Scenario 2 (188 cm) was due in part to (i) increased water
stress as a result of a drier soil (4 cm day−1 rather than 6 cm day−1 irrigation rate, every 5 days),(ii) sub-optimal pressure head conditions owing to a decrease in the hydraulic conductivity inducedby irrigation with high salinity, thus generating near-saturated conditions within the root zone, and(iii) salinity stress (irrigation water EC = 5700 µS cm−1 ).
Cumulative net drainage (the difference between cumulative inflow and outflow) was identical(−84 cm) for Scenarios 3 and 4 (Figure 8c). The higher drainage rate for Scenario 1 (123 cm) was due tothe reduced transpiration, thus leaving more water in the soil profile for possible gravity drainage.
The much lower drainage rate for Scenario 2 (29 cm) was a result of the smaller irrigation rate. Notethat the decrease in irrigation by ~33% resulted in a decrease in drainage of approximately 65%.
The cumulative actual surface and bottom fluxes shown in Figure 8 reveal that the use of lower
salinity irrigation water (Scenarios 3 and 4) produced significantly less salinity stress in that the actual
transpiration rate (280 cm) was nearly equal to its potential value (292 cm), 41 cm more than for
Scenario 1 (untreated water) and 93 cm more than for Scenario 2 (untreated water with possible
reductions in the hydraulic conductivity). As a result, drainage from the bottom of the soil profile
for Scenarios 3 and 4 was 39 cm lower than for Scenario 2 (for the same amount of irrigation water).
Scenarios 3 and 4 hence were more efficient in terms of applied irrigation water use.
Similar important conclusions about the practical implications of salinity management wereobtained in several earlier studies, such as by [43,61]. Ref. [61] used the UnsatChem module todemonstrate that leaching requirements would be lower when estimated with a transient modelling
approach than when using a more standard steady-state model. Adopting leaching requirements
based on the transient approach would hence lead to significant savings in terms of irrigation water
volumes. Ref. [43] showed that while the conventional or water balance approach for estimating
leaching fractions predicts little or no leaching when applied water levels are less than potential ET,
field data and HYDRUS modelling showed considerable leaching around the drip lines. Spatiallyvarying soil wetting patterns during drip irrigation causes localized leaching near the drip lines [43],thus allowing for more profitable production of various crops (e.g., processing tomato) as comparedwith other irrigation methods.

transpiration rate (280 cm) was nearly equal to its potential value (292 cm), 41 cm more than for
Scenario 1 (untreated water) and 93 cm more than for Scenario 2 (untreated water with possible
reductions in the hydraulic conductivity). As a result, drainage from the bottom of the soil profile for
Scenarios 3 and 4 was 39 cm lower than for Scenario 2 (for the same amount of irrigation water).
Water 2017, 9, 385
19 of 33
Scenarios 3 and 4 hence were more efficient in terms of applied irrigation water use.

Figure 8. Cumulative
evaporation(A),transpirationtranspiration (B)(B)
andanddrainage(C) rates
forrates
Scenario
Figure 8. Cumulative
actualactual
evaporation(A),drainage(C)
for2Scenario 2(Sc2),
3(Sc3),
and4(Sc4).
(Sc2), 3 (Sc3), and 4 (Sc4).

3.3. HP1: Trace Metal Transport in Soil Leached with Coal Seam Gas Produced Water
3.3.1. Problem Definition
Our third example uses the coupled one-dimensional HYDRUS-PHREEQC program, or HP1 [14,25],
to simulate the leaching through soil of an accidental spill of produced water containing a complex
mixture of inorganic chemicals and trace metals such as Cd, Cu, Pb, and Zn. The hypothetical example
uses the same soil hydraulic data and dispersion coefficient as the first example, supplemented
with appropriate soil chemical information to allow for the use of all geochemical capabilities of
PHREEQC. Sorption processes in the soil were modelled by a combination of ion exchange and surface
complexation. For details about coupling transport and reaction equations the reader is referredto [14,25].

Water 2017, 9, 385

20 of 33

3.3.2. Soil Geochemical Model for Trace Metal Transport
To estimate the soil’s sorption capacity associated with all relevant sorption processes, the following
assumptions were made. The total sorption capacity was split in equal amounts between adsorption
sites on the surface complex and adsorption sites on the cation exchanger (Table 3). This assumption
was made since no information was available about the contribution of surface complexation to thetotal sorption capacity.
Table 3. Cation exchange capacity and surface complexation capacity for 50% cation exchange sites
and 50% surface complexation sites (see Supplementary Materials for details). Data from [100].
Ca2+

Mg2+

Na+

K+

Total Sorption Capacity

Units

0.279
0.147
0.074

[eq kg−1 ] = [mol(+) kg−1 ]Cation exchange capacity

0.184(mole charge dm−3 soil)

Surface complexation capacity

0.00826(mole charge dm−3 soil)

0.17450.08750.044

0.0891
0.04450.022

0.00260.00260.00130.01290.01290.0065[mol kg−1 ][mol kg−1 ] 50%

The total sorption capacity was obtained from CEC measurements on a soil profile characterised asa Haplic, Self-mulching, Grey Vertosol, in northern New South Wales, Australia [100]. The adsorption
sites of the major cations were calculated based on analytical data of a lysimeter study by [100], and aresummarized in Table 3. The total sorption capacity of 0.147 mol kg−1 was divided in two equal
fractions (i.e., 0.074 mol kg−1 each), with fraction 1 assigned to the ion exchange sites and fraction
2 to the surface complexation sites associated with the weight% of Fe2 O3 in the soil. For use of thesorption capacity in HP1, units were converted from mol·kg−1 to mol dm−3 soil. The cation exchange
capacity of 0.184 mole charge dm−3 soil (for details, see Supplementary Materials) and the surface
complexation capacity of 0.098 mole charge dm−3 soil (obtained by multiplying 0.074 mol kg−1 soil
with bulk density) were used as input to the HP1 biogeochemical model.
The following steps were undertaken to build the soil geochemical model: (1) define initialconcentrations of the trace metals Cd, Cu, Pb, U, and Zn on the sorption complex of the soil (Cu andZn were included since they are potentially important competitors for sorption sites), and (2) define
initial pore-water concentrations of all relevant chemical species (major ions and trace metals).
Adsorbed trace metal concentrations prior to the start of the leaching process were based on
measurements of naturally occurring soil trace metals obtained from a geochemical survey in the study
area [106]. The mobile metal content was determined using the MMI™ technology, which strips the
MMI™ geochemistry mobile metal ions from the surface of soil particles using partial dissolution [106].
Figure 9 shows minimum, median, and maximum concentrations of the 17 most abundant elements(out of a total of 41 elements measured) based on 79 soil samples. Not all of the 17 most abundantelements were included in the current analysis; manganese and iron, for instance, were excluded tosimplify calculations, even though they exhibited the 3rd and 4th highest adsorbed concentration.
This implied a simplification in the soil geochemical model, which needs further corroboration in
future studies. Calculated concentrations of the main naturally occurring trace elements of interest(Cd2+ , Cu2+ , Pb2+ , UO2 2+ , Zn2+ ) on the adsorption complex are available from Table S3.

Water 2017, 9, 385

21 of 33

Figure 9. Mobile metal ion concentrations for 17 most abundant elements (out of 41). The top andbottom of each bar represents maximum and minimum concentrations, respectively; the horizontal
line is the median concentration. (Data from [106]).

Next, the initial concentrations of all relevant chemical species on the cation exchanger were used
to define an initial pore-water solution prior to infiltration of the leachate solution. The geochemical
model defines initial concentrations as follows:
Step 1: All initial concentrations of major cations (Table 3) and trace metals (Table S3) are assigned
to the cation exchanger;
Step 2: Mixing with an initially solute-free pore-water solution takes place to define an initial
pore-water solution with cations and trace metals;
Step 3: Rainwater is allowed to infiltrate, causing the pore-water solution and the rain water toequilibrate with the surface complex and the cation
2 exchanger.
The initial pore-water solution from Step 2 was mixed with a boundary solution typical ofrainwater [99] to achieve a steady-state chemical profile within the soil. In defining the chemicalcomposition of rainfall, with a total annual rainfall of 70 cm y−1 and an estimated recharge of 4 cm y−1 ,concentrations in rainwater were converted to concentrations in recharge water by multiplying theformer by the ratio of rainfall to recharge, i.e., 70/4 or 17.5. In this way, the effect of “concentration” ofsolutes in the top of the profile was accounted for.
The second boundary solution was defined by the composition of the produced water that wasconsidered to infiltrate in the soil as a result of a hypothetical accidental leak from a water holding
pond. Three different compositions of produced water were available from three different sites. For our
example we used the composition of Site 3 (Figure 10).

Water 2017, 9, 385

22 of 33

Figure 10. Simplified chemical composition of produced water used as input solution (only site 3 isreported here).

A total of 18 components (Cl, Ca, Mg, Na, K, Al, C(4), S(6), F, Cd, Zn, Pb, Cu(1), Cu(2) U(5), U(6),H, O) were transported in the model. Although this number simplifies the composition of produced
water, the simplification is justified since other components such as phosphate and nitrate had low tovery low concentrations in the produced water. Table 4 lists, as an example, the aqueous phase species
for lead since this element will be the focus of our discussion later; the remaining trace metal aqueous
phase species for Cd, Cu, U and Zn are available from Table S4. The cation exchange and surface
complexation species considered in the geochemical model are shown in Tables 5 and 6, respectively.
Table 4. Aqueous phase species of lead used in the geochemical model (the phreeqcU.dat database used).
Aqueous Phase Solution Species
Element

Pb

Reaction
Pb2+ + H2 O = PbOH+ + H+
Pb2+ + 2 H2 O = Pb(OH)2 + 2 H+
Pb2+ + 3 H2 O = Pb(OH)3− + 3 H+
Pb2+ + 4 H2 O = Pb(OH)4 2− + 4 H+
2 Pb2+ + H2 O = Pb2 OH3+ + H+
Pb2+ + Cl− = PbCl+
Pb2+ + 2 Cl− = PbCl2
Pb2+ + 3 Cl− 3= PbCl3−
Pb2+ + 4 Cl− = PbCl4 2−
Pb2+ + CO3 2− = PbCO3
Pb2+ + 2 CO3 2− = Pb(CO3 )2 2−
Pb2+ + HCO3 − = PbHCO3 +Pb2+ + SO4 2− = PbSO4
Pb2+ + 2 SO4 2− = Pb(SO4 )2 2−
Pb2+ + NO3 − = PbNO3
Pb2+ + F− = PbF+
Pb2+ + 2 F− = PbF2

logk
−7.71
−17.12−28.06
−39.7
−6.36
1.6
1.8
1.7
1.38
7.24
10.64
2.9
2.75
3.47
1.17
1.25
2.56

Water 2017, 9, 385

23 of 33

Table 5. Cation exchange species of the heavy metals included in the geochemical model(the phreeqcU.dat database used).
Cation Exchange
Element
Cd
Cu
Pb
UO2
Zn

Reaction
Cd2+ + 2X− = CdX2
Cu2+ + 2X− = CuX2
Pb2+ + 2X− = PbX2
UO2 2+ + 2X− = UO2 X2
Zn2+ + 2X− = ZnX2

logk
0.8
0.6
1.05
0.8
0.8

Table 6. Surface complexation species of the heavy metals included in the geochemical model(the phreeqcU.dat database used). Note that also H+ and other cation and anions may form
surface complexes.
Surface Complexation
Element
Cd
Cu
Pb
UO2
Zn

Reaction
Hfo_wOH + Cd2+ = Hfo_wOCd+ + H+
Hfo_wOH + Cu2+ = Hfo_wOCu+ + H+
Hfo_wOH + Pb2+ = Hfo_wOPb+ + H+
Hfo_wOH + UO2 2+ = Hfo_wOUO2 + + H+
Hfo_wOH + Zn2+ = Hfo_wOZn+ + H+

logk
−2.91
0.6
0.3
2.8
−1.99

The following three scenarios are considered illustrative examples of multi-component trace metaltransport in soil following a long-term leak from a water holding pond:
Scenario 1 assumes rainwater infiltration only with a flux of 4 cm y−1 over a 100 year period,
and with a rainwater chemical composition as given in Table S5. This allows a comparison ofsoil geochemical behaviour under natural conditions of rainwater infiltration with the scenariosof infiltration of produced water with a very specific hydrochemistry (Scenarios 2 and 3).
Scenario 2 considers infiltration of produced water from which trace metals were excluded, thus
keeping only major ions (site 3, Figure 10). This analysis was undertaken to evaluate trace metal
behaviour in soils where the naturally present adsorbed trace metals (the background concentration
shown in Figure 9) were exposed to an infiltrating alkaline solution (the produced water infiltrating
the soil has a pH of 9.14). This provides a basis for comparisons with trace metal behaviour whenmetal-containing produced water infiltrates.
Scenario 3 includes all measured trace metals in the produced water infiltrating the soil (site 3,
Figure 10).
To simplify the interpretation of trace metal chemical behavior, a 100-cm deep soil profile witha single soil material was considered, with the soil hydraulic properties of layer 1 from Table S2.
The profile was discretized into 100 equidistant finite elements. Soil sorption was hypothesized
to consist of 50% adsorption sites with surface complexation and 50% adsorption sites with cationexchange. The cation exchange and surface complexation capacities are summarised in Table 3.
Time-variable top boundary conditions were defined as follows: the first period consisted of 100 yearsof rainwater infiltration during which a water holding pond was absent, with the chemical compositionof the rainwater defined in Table S5. After this initial “warming up” period that allows chemical
equilibration between soil and rainwater to be established, the actual infiltration of produced water
took place over a 30-year period during which a water holding pond was present and leaked at a verysmall rate of 0.35 cm y−1 . The chemical composition of the Site 3 produced water was used as input
to the simulations (Figure 10). The last period was defined by infiltration of rainwater for another70 years, thus assuming that the water holding pond had been dismantled, using the same chemicalcomposition as the initial warming-up period. The rainwater infiltration flux was always 4 cm y−1 .

Water 2017, 9, 385

24 of 33

3.3.3. HP1 Simulation Results
The discussions below focus on the following aspects from the wealth of output provided by
the HP1 simulations: (i) changes in adsorbed concentrations of lead on the cation exchanger and thesurface complex, (ii) changes in liquid phase concentrations for the same trace metal, and (iii) changes
in trace metal fluxes at the bottom of the 1-m deep soil profile. The detailed analysis of the trace
metal behavior uses eight time steps at which adsorbed and aqueous phase concentrations are plotted.
The time corresponding to 100 years of rainwater infiltration, after which equilibrium in chemistry
between rainwater and soil water was established, was defined as t0 . Chemical concentrations within
the soil profile at t0 represented steady-state conditions for Scenario 1. Three time steps (10 [t10 ],20 [t20 ], and 30 years [t30 ]) were defined to capture changes during the 30-year infiltration period withproduced water. To evaluate the progression of chemical leaching once infiltration with produced
water terminated and rainwater infiltrated again (from t > 30 to t = 100 years), short (5-year) time steps
t35 , t40 and t45 were used to capture changes that would occur within 15 years of rainfall infiltration.
The final time level t100 , corresponded to the total simulation time of 200 years, i.e. 100 years after thestart of the infiltration of produced water.
Scenario 1 is represented by a single curve (i.e., the initial concentration at t0) in Figures 11 and 12).
The effect of infiltration with produced water (Scenarios 2 and 3) produced less adsorbed lead on the
cation exchanger and the surface complexation sites (Figure 12), which in turn decreased aqueous
lead owing to equilibrium between the solid and liquid phase concentrations (Figure 11). At the endof the simulation (t100 ), re-adsorption of lead occurred on both adsorption sites between depths of
0.4 and 1 m. The lower liquid phase concentration at t100 is a reflection of the mass loss of lead fromthe adsorption sites; and the equilibrium partitioning of lead between solid and liquid phase according
to the equations in Tables 5 and 6.

Figure 11. Liquid phase concentrations for lead in a 1-m deep soil profile for 30-y infiltration
with produced water (comparison between Scenario 2 [trace metal-free, solid lines] and Scenario 3[trace-metal containing water [TM, symbols]). Times are as follows: T = 0 y and T = 30 y are the startand end of infiltration with produced water; between T > 30 y and T = 100 y rainwater infiltrated.

4

Water 2017, 9, 385

25 of 33

Figure 12. Solid phase concentrations for lead in a 1-m deep soil profile for a 30-y infiltration withproduced water (comparison between trace metal-free and trace-metal containing [TM] water). Times
are as follows: T = 0 y and T = 30 y are the start and end of infiltration with produced water; between T> 30 y and T = 100 y rainwater infiltrated.

4
A comparison between Scenarios 2 (infiltration with metal-free water) and 3 (infiltration withmetal-containing water) shows no significant difference in lead concentrations in the soil profile,neither for the liquid (Figure 11) or adsorbed phase (Figure 12). Most of the total adsorbed lead resided
on the surface complex before and during infiltration of the produced water; some desorption fromthe surface complex is noticeable close to the surface (Figure 11). Strong desorption of Pb occurredon the cation exchanger, with only slow re-adsorption after t30 . The concentration in the liquid phase
initially increased to a maximum of about 10−7 mol L−1 up to 0.5 m depth at t30 (due to desorption ofPb from the cation exchanger), then decreased again due to leaching with rainwater such that at t100
Pb was removed from the pore-water with final concentrations smaller than at those at t0 . (Figure 11).
Our findings indicate that the main processes of trace metal sorption, desorption and transport
are not linked to metals present in the produced water. The observed concentration distributions are
the result of complexation reactions between inorganic compounds present in the alkaline produced
water and the trace metals present on the solid phases prior to infiltration. To evaluate which dissolved
species from the infiltrating solution are responsible for the desorption and transport of trace metals(for a full list of lead species considered, see Table 4), and how they affect the dynamics of trace metal
sorption/desorption during the leakage period, the most abundant liquid phase lead species at t0 andat t25 are displayed in Figure 13. To appreciate the results, changes in the pH also need consideration:
there is generally a significant increase in solution pH from initial values around 7.2–7.4 to nearly 9.
Prior to infiltration with produced water (at t0 ), about 60% of all dissolved lead species wererepresented by the cation Pb2+ , about 20% by the cation PbOH+ , and approximately 15–20% by theneutral PbCO3 (Figure 13). Smaller fractions were due to the species Pb(CO3 )2 2− , PbHCO3 + , PbSO4 ,and PbCl+ . After infiltration with produced water (species distribution shown at t25 ), the maximum
percentage for Pb2+ had decreased to 20%. The predominant species was then the neutral species
PbCO3 with a maximum percentage of approximately 80%. Furthermore, the abundance of two othercarbonate species (Pb(CO3 )2 2− , PbHCO3 + ) increased compared to the initial condition. These neutral
and negatively charged lead species have a higher mobility than the positively charged Pb2+ , causing
the latter to desorb and to be transported down the profile as predominantly PbCO3 , Pb(CO3 )2 2−
and PbHCO3 + .

dissolved species from the infiltrating solution are responsible for the desorption and transport oftrace metals (for a full list of lead species considered, see Table 4), and how they affect the dynamics
of trace metal sorption/desorption during the leakage period, the most abundant liquid phase lead
species at t0 and at t25 are displayed in Figure 13. To appreciate the results, changes in the pH also
needconsideration:
there is generally a significant increase in solution pH from initial values around
Water
2017,
9, 38526 of 33
7.2–7.4 to nearly 9.
0

0

T0

T25
-0.2

Soil depth (m)Soil depth (m)

-0.2

-0.4

-0.6

Pb2+
PbOH+
PbCO3
Pb(CO3)22PbHCO3+
PbSO4
PbCl+

-0.8

-1
0

20

-0.4

-0.6

-0.8

-1
40

60

80

Liquid phase lead species (%)

100

0

20

40

60

80

100

Liquid phase lead species (%)

Figure 13. Depth distribution of the most abundant lead species in the liquid phase at t0 (initial
Figure13. Depth distribution of the most abundant lead species in the liquid phase at t0 (initial
equilibrium condition, left) and t25 (infiltration of produced water, right).
equilibrium condition, left) and t25 (infiltration of produced water, right).

Prior to infiltration with produced water (at t0), about 60% of all dissolved lead species were
Further
for nearly
unretarded
transportof PbOHthe dominantliquid phase 15–20%
species is+, and approximately
representedevidenceby the cation
Pb2+, about
20% bythe cation
byfoundtheinneutral
a comparisonofthesolutefluxesoflead(totallead
concentration
based
on
all
dissolved
2−+
PbCO3 (Figure 13). Smaller fractions were due to the species Pb(CO3)2 , PbHCO3 , PbSOspecies)
4, and+. Afterandchloride
atinfiltration
the bottomwithof theproduced
soil profile.
Figure14 displays
the bottom
fluxesof total lead
PbCl
water
(species
distribution
shownsoluteat t25),the maximum2+ hadandchloride: for
thePb
peak
concentration
appeared
at predominant
the same timespecies
for bothchemicals,
35 years
after
percentage
decreased to20%. The
wasthen the i.e.,
neutralspeciesthebeginningoftheleakage.
This
confirms
thattheneutralor
negatively
chargedlead
species
PbCO3 with a maximum percentage of approximately 80%. Furthermore, the abundance of two otheraretransported
unretarded
throughthe soil3+)profile.
carbonatespecies(Pb(CO3)22−, PbHCO
increased compared to the initial condition. These neutral
Water
2017, 9, 385
27 of 34
and negatively charged lead species have a higher mobility than the positively charged Pb2+, causing
the latter to desorb and to be transported down the profile as predominantly PbCO3, Pb(CO3)22− andPbHCO3+.
Further evidence for nearly unretarded transport of the dominant liquid phase species is foundin a comparison of the solute fluxes of lead (total lead concentration based on all dissolved species)and chloride at the bottom of the soil profile. Figure 14 displays the bottom solute fluxes of total leadand chloride: the peak concentration appeared at the same time for both chemicals, i.e., 35 years after
the beginning of the leakage. This confirms that the neutral or negatively charged lead species aretransported unretarded through the soil profile.

Figure 14. Lead and chloride bottom solute fluxes as a function of time. Comparison between trace
Figure 14. Lead and chloride bottom solute fluxes as a function of time. Comparison between trace
metal-free and trace metal-containing produced water.
metal-free and trace metal-containing produced water.

We made a final comparison of lead fluxes obtained for Scenarios 2 (metal free) and 3 (metal
We madea final comparisonof leadfluxesobtained(metal free)
3 (metalcontaining)by considering
the bottom
fluxesof totallead atfor
1-mScenarios
depth. As2 expected,
no anddifferenceincontaining)byconsidering
thebottomfluxes
of
total
lead
at
1-m
depth.
As
expected,nodifference
inthe flux occurred for the two scenarios. The peak of the bottom solute flux at 1 m was reached 35
the flux
occurred
forofthescenarios.ofThepeak ofwater
the bottom
soluteflux at(Figure
1 m was14).
reached35 years
years
after
the start
thetwoinfiltration
produced
for bothscenarios
This confirmsthat lead leaching is due to desorption of naturally present lead species following the infiltration ofan alkaline solution; the contribution of trace metals, including those present in produced water, did
not materially contribute to the final lead flux.
4. Summary and Conclusions

Water 2017, 9, 385

27 of 33

after the start of the infiltration of produced water for both scenarios (Figure 14). This confirms thatlead leaching is due to desorption of naturally present lead species following the infiltration of analkaline solution; the contribution of trace metals, including those present in produced water, did notmaterially contribute to the final lead flux.
4. Summary and Conclusions
In this paper we demonstrated the use of the HYDRUS software and its specialized add-on
modules UnsatChem (for major ion chemistry) and HP1 (for more general multicomponent transport)
to evaluate the fate and subsurface transport of various chemicals involved in coal seam gas extraction
and water management operations after their accidental or intentional release at the soil surface.
The three examples we discussed represented increasing model complexity from the relatively standard
HYDRUS (example 1), to using the UnsatChem module (example 2), and finally HP1 (example 3).
The first example evaluated the fate and transport of organic chemicals (i.e. biocide) in soil as a resultof their accidental release from water holding ponds containing produced water. The standard HYDRUS
software was used to simulate the coupled processes of retardation, first-order degradation andconvective-dispersive transport of the biocide bronopol and its degradation products. We presented
a review of degradation pathways and relevant biogeochemical processes, including a collation ofphysico-chemical properties of bronopol and two of its degradation products (2-bromo-2-nitroethanoland bromonitromethane) as input to a set of generic simulations of transport and attenuation invariably-saturated soil profiles. We demonstrated the ability to model the coupled processes of fluid
flow and transport of multiple contaminants in soils, with a sensitivity analysis testing the robustnessof contaminant migration and attenuation with regards to chemical parameters.
In a second example we used the major ion chemistry module UnsatChem to simulate thesimultaneous movement of multiple major ions (Na, K, Ca, Mg, Cl, SO4 , alkalinity) present in irrigationwater sourced from coal seam gas produced water. This analysis illustrated the identification andpotential optimisation of irrigation water quality that will minimize long-term harmful effects ofdissolved ions on soil and plant health. Three different irrigation water qualities were evaluated:
untreated or treated produced water or produced water amended with surface water. Simulations
with the different irrigation waters provided detailed results regarding chemical indicators of soil andplant health, notably SAR, EC and sodium concentrations. By comparing these indicators in the soil
profile with water quality guideline values (i.e., the Australian and New Zealand ANZECC values),
an assessment could be made of the suitability of the applied produced water for long-term irrigation.
The major ion chemistry module UnsatChem implemented in HYDRUS-1D also allowed for testing
the effects of salinity on root-water uptake and thus on crop transpiration. The example additionallypermitted calculations of the effects of salinity on the soil hydraulic conductivity and how reductionsin this key soil property can lead to water logging and hence reduced root-water uptake. UnsatChem
offers a cost-effective way to optimise the quality of irrigation water derived from coal seam gas wateror other saline waters and thus to ensure that soils are managed in a sustainable manner.
The third example used HP1 to determine key processes contributing to the migration of naturally
occurring trace metals cadmium, lead, and uranium present in a particular soil profile. Based ona discussion that focused on the trace metal lead, chemical transport and attenuation processes insoil were simulated by including adsorption of trace metals through cation exchange and surface
complexation. The assessment was based on an analysis, in space and time, of adsorbed and aqueous
phase metal concentrations and the relative abundance of dissolved metal species. Results show thatthe main process responsible for trace metal migration in soil is complexation of naturally present
metals with inorganic ligands such as (bi)carbonate, chloride, and hydroxyl ions. These ligands enter
the soil upon infiltration with coal seam gas produced water, which is typically alkaline (pH = 8–9).
Trace metals naturally present on adsorption sites in the soil, even at low levels, may desorb due toformation of especially mobile carbonate metal-ligand compounds. These metal-ligand compounds
travel relatively unretarded through soil. The multi-component reactive transport simulator HP1 wasWater 2017, 9, 385

28 of 33

shown to be a powerful tool to assess contamination risk from release of industrial waters. The ability
to incorporate coupled physical, chemical, geological and biological processes in the simulations also
allows testing of mitigating measures as part of pollution prevention.
Supplementary Materials: The following are available online at www.mdpi.com/2073-4441/9/6/385/s1.
Table S1: Organic carbon partition coefficient KOC (L/kg) and aerobic degradation half-lives (T1/2 , days) of
bronopol and two of its degradation products; Table S2: van Genuchten soil hydraulic parameters and the bulk
density of Vertosols (averages across the Cox’s Creak Catchment area, NSW, Australia) (source: [7]); Table S3:
Conversion of the concentrations [ppm] measured by de Caritat and Lech [9] to the initial concentrations [mol/kg]
on the cation exchanger; Table S4: Aqueous phase species of the heavy metals included in the model (phreeqcU.dat
database used). Table S4: Rainfall data and chemical composition of rainwater at Wagga Wagga, Australia(Source: [9]). HP1 input data calculated for an infiltration flux of 40 mm/y.
Acknowledgments: We are grateful to Boris Faybishenko (Lawrence Berkeley National Laboratory), Guest Editorof the Special Issue “Water and Solute Transport in the Vadose Zone”, for the invitation to present our paper.
Thanks also to the three anonymous reviewers for their constructive comments. The first author acknowledges
CSIRO Land and Water for funding received through the strategic appropriation research project “Next generation
methods and capability for multi-scale cumulative impact and management”.
Author Contributions: Dirk Mallants collated the data and performed the simulations; Dirk Mallants, Jirka
Šimůnek, Martinus Th. van Genuchten and Diederik Jacques contributed to the writing of the paper.
Conflicts of Interest: The authors declare no conflict of interest.

References
1.
2.

3.
4.
5.

6.

7.
8.

9.
10.

11.
12.

Singh, G.; Kaur, G.; Williard, K.; Schoonover, J.; Kang, J. Monitoring of Water and Solute Transport in the
Vadose Zone: A Review. Vadose Zone J. 2017. [CrossRef]
US EPA. Hydraulic Fracturing for Oil and Gas: Impacts from the Hydraulic Fracturing Water Cycle on Drinking
Water Resources in the United States; Office of Research and Development: Washington, DC, USA, 2016;
EPA/600/R-16/236Fa.
Christensen, T.H.; Kjeldsen, P.; Bjerg, P.L.; Jensen, D.L.; Christensen, J.B.; Baun, A.; Albrechtsen, H.-J.;
Heron, G. Biogeochemistry of landfill leachate plumes. Appl. Geochem. 2001, 16, 659–718. [CrossRef]
Akcil, A.; Koldas, S. Acid mine drainage (AMD): Causes, treatment and case studies. J. Clean. Prod. 2006, 14,1139–1145. [CrossRef]
Fredrickson, J.K.; Zachara, J.M.; Balkwill, D.L.; Kennedy, D.; Li, S.-M.W.; Kostandarithes, H.M.; Daly, M.J.;
Romine, M.F.; Brockman, F.J. Geomicrobiology of High-Level Nuclear Waste-Contaminated Vadose
Sediments at the Hanford Site, Washington State. Appl. Environm. Microbiol. 2004, 70, 4230–4241. [CrossRef][PubMed]Vandenhove, H.; Sweeck, L.; Mallants, D.; Vanmarcke, H.; Aitkulov, A.; Sadyrov, O.; Savosin, M.;
Tolongutov, B.; Mirzachev, M.; Clerc, J.J.; et al. Assessment of radiation exposure in the uranium andmilling area of Mailuu Suu, Kyrgyzstan. J. Environ. Radioact. 2006, 88, 118–139. [CrossRef] [PubMed]
Lapworth, D.J.; Baran, N.; Stuart, M.E.; Ward, R.S. Emerging organic contaminants in groundwater: A review
of sources, fate and occurrence. Environ. Pollut. 2012, 163, 287–303. [CrossRef] [PubMed]
Mallants, D.; Diels, L.; Bastians, L.; Vos, J.; Moors, H.; Wang, L.; Maes, N.; Vandenhove, H. Removal of
Uranium and Arsenic from Groundwater Using Six Different Reactive Materials: Assessment of Removal
Efficiency. In Uranium in the Aquatic Environment; Merkel, B.J., Planer-Friedrich, B., Wolkersdolfer, C., Eds.;
Springer: Berlin, Germany, 2002; pp. 561–568.
Scherer, M.M.; Richter, S.; Valentine, R.L.; Alvarez, P.J.J. Chemistry and microbiology of permeable reactive
barriers for in situ groundwater clean up. Crit. Rev. Environ. Sci. Technol. 2000, 30, 363–411. [CrossRef]Šimůnek, J.; Jacques, D.; Langergraber, G.; Bradford, S.A.; Šejna, M.; van Genuchten, M.T. Numerical
Modelling of Contaminant Transport Using HYDRUS and Its Specialized Modules, Invited paper for theSpecial Issue “Water Management in Changing Environment”. J. Indian Inst. Sci. 2013, 93, 265–284.
Šimůnek, J.; van Genuchten, M.T.; Šejna, M. Recent developments and applications of the HYDRUS computer
software packages. Vadose Zone J. 2016, 15, 25. [CrossRef]Šimůnek, J.; Bradford, S. Vadose Zone Modelling: Introduction and importance. Vadose Zone J. 2008, 7,
581–586. [CrossRef]

Water 2017, 9, 385

13.

14.

15.

16.
17.

18.

19.
20.
21.
22.
23.
24.

25.

26.
27.
28.
29.
30.
31.

32.
33.

29 of 33

Seuntjens, P.; Mallants, D.; Šimůnek, J.; Patyn, J.; Jacques, D. Sensitivity analysis of physical and chemical
properties affecting field-scale cadmium transport in a heterogeneous soil profile. J. Hydrol. 2002, 264,
185–200. [CrossRef]
Jacques, D.; Šimůnek, J.; Mallants, D.; van Genuchten, M.T. Modelling coupled hydrologic and chemical
processes: Long-term Uranium transport following Phosphorus fertilization. Vadose Zone J. 2008, 7, 698–711.
[CrossRef]
Mallants, D.; Jacques, D.; Zeevaert, T. Modelling 226 Ra, 222 Rn, and 210 Pb Migration in a Proposed Surface
Repository of Very Low-Level Long-Lived Radioactive Waste. In Proceedings of the ICEM 2003, the 9thInternational Conference on Radioactive Waste Management and Environmental Remediation, Oxford, UK,
21–25 September 2003; icem03-4632, p. 8.
Li, Y.; Šimůnek, J.; Zhang, Z.; Jing, L.; Ni, L. Evaluation of nitrogen balance in a direct-seeded-rice field
experiment using Hydrus-1D. Agric. Water Manag. 2015, 148, 213–222. [CrossRef]
Papiernik, S.K.; Yates, S.R.; Koskinen, W.C.; Barber, B. Processes affecting the dissipation of the herbicide
isoxaflutole and its diketonitrile metabolite in agricultural soils under field conditions. J. Agric. Food Chem.
2007, 55, 8630–8639. [CrossRef] [PubMed]
Schaerlaekens, J.; Mallants, D.; Šimůnek, J.; van Genuchten, M.T.; Feyen, J. Numerical simulation of transportand sequential biodegradation of chlorinated aliphatic hydrocarbons using CHAIN_2D. Hydrol. Proc. 1999,
13, 2847–2859. [CrossRef]
Fan, Z.; Casey, F.X.M.; Hakk, H.; Larsen, G.L. Discerning and modelling the fate and transport of testosterone
in undisturbed soil. J. Environ. Qual. 2007, 36, 864–873. [CrossRef] [PubMed]
Unold, M.; Šimůnek, J.; Kasteel, R.; Groeneweg, J.; Vereecken, H. Transport of manure-based applied
sulfadiazine and its main transformation products in soil columns. Vadose Zone J. 2009, 8, 677–689. [CrossRef]Šimůnek, J.; Jacques, D.; van Genuchten, M.T.; Mallants, D. Multicomponent geochemical transport modellingusing Hydrus-1D and HP1. J. Am. Water Resour. Assoc. 2006, 42, 1537–1547. [CrossRef]Šimůnek, J.; Suarez, D.L. Two-dimensional transport model for variably saturated porous media with major
ion chemistry. Water Resour. Res. 1994, 30, 1115–1133. [CrossRef]Šimůnek, J.; Suarez, D.L. Sodic soil reclamation using multicomponent transport modelling. J. Irrig.
Drain. Eng. 1997, 123, 367–376. [CrossRef]
Parkhurst, D.L.; Appelo, C.A.J. User’s Guide to PHREEQ C (Version 2)—A Computer Program for Speciation,
Batch-Reaction, One-Dimensional Transport and Inverse Geochemical Calculations. In Water-Resources
Investigations; Report 99-4259; U.S. Geological Survey: Earth Science Information Center: Denver, CO, USA,
1999; p. 312.
Jacques, D.; Šimůnek, J.; Mallants, D.; van Genuchten, M.T. Modelling coupled water flow, solute transportand geochemical reactions affecting heavy metal migration in a podzol soil. Geoderma 2008, 145, 449–461.
[CrossRef]
Langergraber, G.; Šimůnek, J. Modelling variably-saturated water flow and multi-component reactive
transport in constructed wetlands. Vadose Zone J. 2005, 4, 924–938. [CrossRef]
Langergraber, G.; Šimůnek, J. Reactive transport modelling of subsurface flow constructed wetlands. Vadose
Zone J. 2012, 11. [CrossRef]
Langergraber, G.; Rousseau, D.; García, J.; Mena, J. CWM1—A general model to describe biokinetic processes
in subsurface flow constructed wetlands. Water Sci. Technol. 2009, 59, 1687–1697. [CrossRef] [PubMed]Šimůnek, J.; Changming, H.; Pang, L.; Bradford, S.A. Colloid-facilitated transport in variably-saturated
porous media: Numerical model and experimental verification. Vadose Zone J. 2006, 5, 1035–1047. [CrossRef]
Pang, L.; Šimůnek, J. Evaluation of bacteria-facilitated cadmium transport in gravel columns using the
HYDRUS colloid-facilitated solute transport model. Water Resour. Res. 2006, 42, W12S10. [CrossRef]
Zhang, M.; Engelhardt, I.; Šimůnek, J.; Bradford, S.A.; Kasel, D.; Berns, A.E.; Vereecken, H.; Klumpp, E.
Co-transport of chlordecone and sulfadiazine in the presence of functionalized multi-walled carbon
nanotubes in soils. Environ. Pollut. 2017, 221, 470–479. [CrossRef] [PubMed]
Spurlock, F.; Šimůnek, J.; Johnson, B.; Tuli, A. Sensitivity analysis of vadose zone fumigant transport andvolatilization. Vadose Zone J. 2013, 12, 12. [CrossRef]
Spurlock, F.; Johnson, B.; Tuli, A.; Gao, S.; Tao, J.; Sartori, F.; Qin, R.; Sullivan, D.; Stanghellini, M.; Ajwa, H.
Simulation of fumigant transport and volatilization from tarped broadcast applications. Vadose Zone J. 2013,12, 10. [CrossRef]

Water 2017, 9, 385

34.
35.
36.
37.
38.
39.
40.

41.
42.

43.
44.

45.
46.
47.
48.
49.

50.

51.
52.
53.
54.
55.
56.

30 of 33

Brooks, R.H.; Corey, A.T. Hydraulic Properties of Porous Media. In Hydrology Paper No. 3; Colorado State
University: Fort Collins, CO, USA, 1964.
van Genuchten, M.T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.
Soil Sci. Soc. Am. J. 1980, 44, 892–898. [CrossRef]
Vogel, T.; Císlerová, M. On the reliability of unsaturated hydraulic conductivity calculated from the moistureretention curve. Transp. Porous Med. 1988, 3, 1–15. [CrossRef]
Kosugi, K. Lognormal distribution model for unsaturated soil hydraulic properties. Water Resour. Res. 1996,
32, 2697–2703. [CrossRef]
Durner, W. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour. Res.
1994, 30, 211–233. [CrossRef]Šimůnek, J.; Hopmans, J.W. Modelling compensated root water and nutrient uptake. Ecol. Model. 2009, 220,505–521. [CrossRef]Šimůnek, J.; Jarvis, N.J.; van Genuchten, M.T.; Gärdenäs, A. Review and comparison of models for describing
non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. 2003, 272, 14–35.
[CrossRef]
Saito, H.; Šimůnek, J.; Mohanty, B. Numerical analyses of coupled water, vapor, and heat transport in thevadose zone. Vadose Zone J. 2006, 5, 784–800. [CrossRef]
Pontedeiro, E.M.; van Genuchten, M.T.; Cotta, R.M.; Šimůnek, J. The effects of preferential flow and soil
texture on risk assessments of a NORM waste disposal site. J. Hazard. Mater. 2010, 174, 648–655. [CrossRef][PubMed]Hanson, B.R.; Šimůnek, J.; Hopmans, J.W. Numerical modelling of urea-ammonium-nitrate fertigation under
micro-irrigation. Agric. Water Manag. 2006, 86, 102–113. [CrossRef]
Ramos, T.B.; Šimůnek, J.; Gonçalves, M.C.; Martins, J.C.; Prazeres, A.; Pereira, L.S. Two-dimensionalmodelling of water and nitrogen fate from sweet sorghum irrigated with fresh and blended saline waters.
Agric. Water Manag. 2012, 111, 87–104. [CrossRef]
Casey, F.X.M.; Šimůnek, J. Inverse analyses of the transport of chlorinated hydrocarbons subject to sequential
transformation reactions. J. Environ. Qual. 2001, 30, 1354–1360. [CrossRef] [PubMed]
Casey, F.X.M.; Larsen, G.L.; Hakk, H.; Šimůnek, J. Fate and transport of 17β-Estradiol in soil-water systems.
Environ. Sci. Technol. 2003, 37, 2400–2409. [CrossRef] [PubMed]
Casey, F.X.M.; Larsen, G.L.; Hakk, H.; Šimůnek, J. Fate and transport of testosterone in agriculturally
significant soils. Environ. Sci. Technol. 2004, 38, 790–798. [CrossRef] [PubMed]
Casey, F.X.M.; Šimůnek, J.; Lee, J.; Larsen, G.L.; Hakk, H. Sorption, mobility, and transformation of estrogenic
hormones in natural soil. J. Environ. Qual. 2005, 34, 1372–1379. [CrossRef] [PubMed]
Das, B.S.; Lee, L.S.; Rao, P.S.C.; Hultgren, R.P. Sorption and degradation of steroid hormones in soils duringtransport: Column studies and model evaluation. Environ. Sci. Technol. 2004, 38, 1460–1470. [CrossRef][PubMed]
Wehrhan, A.; Kasteel, R.; Šimůnek, J.; Groeneweg, J.; Vereecken, H. Transport of sulfadiazine in soilcolumns—experiments and modelling approaches. J. Contam. Hydrol. 2007, 89, 107–135. [CrossRef][PubMed]Dontsova, K.M.; Yost, S.L.; Šimůnek, J.; Pennington, J.C.; Williford, C. Dissolution and transport of TNT, RDX,
and Composition B in saturated soil columns. J. Environ. Qual. 2006, 35, 2043–2054. [CrossRef] [PubMed]
Dontsova, K.M.; Pennington, J.C.; Hayes, C.; Šimůnek, J.; Williford, C.W. Dissolution and transport of
2,4-DNT and 2,6-DNT from M1 propellant in soil. Chemosphere 2009, 77, 597–603. [CrossRef] [PubMed]Šimůnek, J.; van Genuchten, M.T. Modeling Nonequilibrium Flow and Transport Processes Using HYDRUS.
Vadose Zone J. 2008, 7, 782–797. [CrossRef]Šimůnek, J.; Valocchi, A.J. Geochemical Transport. In Methods of Soil Analysis, Part 1, Physical Methods,
Chapter 6.9, 3rd ed.; Dane, J.H., Topp, G.C., Eds.; SSSA: Madison, WI, USA, 2002; pp. 1511–1536.
Kaledhonkar, M.J.; Keshari, A.K. Modelling the effects of saline water use in agriculture. Irrig. Drain. 2006,
55, 177–190. [CrossRef]
Gonçalves, M.C.; Šimůnek, J.; Ramos, T.B.; Martins, J.C.; Neves, M.J.; Pires, F.P. Multicomponent solutetransport in soil lysimeters irrigated with waters of different quality. Water Resour. Res. 2006, 42, W08401.
[CrossRef]

Water 2017, 9, 385

57.

58.
59.
60.

61.
62.

63.
64.

65.

66.

67.

68.

69.
70.
71.
72.

73.
74.

75.

76.

31 of 33

Kaledhonkar, M.J.; Keshari, A.K.S.; van der Zee, E.A.T.M. Relative sensitivity of ESP profile to spatial andtemporal variability in cation exchange capacity and pore water velocity under simulated field conditions.
Agric. Water Manag. 2006, 83, 58–68. [CrossRef]
Kaledhonkar, M.J.; Sharma, D.R.; Tyagi, N.K.; Kumar, A.; van der Zee, S.E.A.T.M. Modelling for conjunctive
use irrigation planning in sodic groundwater areas. Agric. Water Manag. 2012, 107, 14–22. [CrossRef]
Skaggs, T.H.; Shouse, P.J.; Poss, J.A. Irrigation of forage crops with saline drainage waters: 2. Modellingdrainage and root water uptake. Vadose Zone J. 2006, 5, 824–837. [CrossRef]
Schoups, G.; Hopmans, J.W.; Tanji, K.K. Evaluation of model complexity and space–time resolution on theprediction of long-term soil salinity dynamics, western San Joaquin Valley, California. Hydrol. Process. 2006,
20, 2647–2668. [CrossRef]
Corwin, D.L.; Rhoades, J.D.; Šimůnek, J. Leaching requirement for soil salinity control: Steady-state vs.
transient-state models. Agric. Water Manag. 2007, 90, 165–180. [CrossRef]
Ramos, T.B.; Šimůnek, J.; Gonçalves, M.C.; Martins, J.C.; Prazeres, A.; Castanheira, N.L.; Pereira, L.S. Field
evaluation of a multicomponent solute transport model in soils irrigated with saline waters. J. Hydrol. 2011,
407, 129–144. [CrossRef]
Rasouli, F.; Pouya, A.K.; Šimůnek, J. Modelling the effects of saline water use in wheat-cultivated lands using
the UNSATCHEM model. Irrig. Sci. 2013, 31, 1009–1024. [CrossRef]
Jakubowski, R.; Haws, N.; Ellerbrock, D.; Murtagh, J.; MacFarlane, D. Development of a management
tool to support the beneficial use of treated coal seam gas water for irrigation in Eastern Australia.
Mine Water Environ. 2014, 33, 133–145. [CrossRef]
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.; et al. Reactive transport codes for subsurface environmental simulation.
Comput. Geosci. 2015, 19, 445–447. [CrossRef]
Mayer, K.U.; Alt-Epping, P.; Jacques, D.; Arora, B.; Steefel, C.I. Benchmark problems for reactive transportmodeling of the generation and attenuation of acid rock drainage. Comput. Geosci. 2015, 19, 599–611.
[CrossRef]
Xie, M.; Mayer, K.U.; Claret, F.; Alt-Epping, P.; Jacques, D.; Steefel, C.; Chiaberge, C.; Simunek, J.
Implementation and evaluation of permeability-porosity and tortuosity-porosity relationships linked tomineral dissolution-precipitation. Comput. Geosci. 2015, 19, 655–671. [CrossRef]
Jacques, D.; Šimůnek, J.; Mallants, D.; van Genuchten, M.T. The HPx reactive transport models: Summaryof recent developments and applications. In Proceedings of the 4th International Conference “HYDRUS
Software Applications to Subsurface Flow and Contaminant Transport Problems”, Prague, Czech Republic,
21–22 March 2013; Šimůnek, J., van Genuchten, M.T., Kodešová, R., Eds.; Czech University of Life Sciences:
Prague, Czech Republic; pp. 7–16.
Leterme, B.; Blanc, P.; Jacques, D. A reactive transport model for mercury fate in soil—application to different
anthropogenic pollution sources. Environ. Sci. Pollut. R. 2014, 21, 12279–12293. [CrossRef] [PubMed]
Leterme, B.; Jacques, D. A reactive transport model for mercury fate in contaminated soil-sensitivity analysis.
Environ. Sci. Pollut. R. 2015, 22, 16830–16842. [CrossRef] [PubMed]
Bessinger, B.A.; Marks, C.D. Treatment of Mercury-Contaminated Soils With Activated Carbon: A Laboratory,
Field, and Modelling Study. Remediat. J. 2010, 21, 115–135. [CrossRef]
Thaysen, E.M.; Jacques, D.; Jessen, S.; Andersen, C.E.; Laloy, E.; Ambus, P.; Postma, D.; Jakobsen, I. Inorganiccarbon fluxes across the vadose zone of planted and unplanted soil mesocosms. Biogeosciences 2014, 11,7179–7192. [CrossRef]
Thaysen, E.M.; Jessen, S.; Postma, D.; Jakobsen, R.; Jacques, D.; Ambus, P.; Laloy, E.; Jakobsen, I. Effects oflime and concrete waste of vadose zone carbon cycling. Vadose Zone J. 2014, 13, 11. [CrossRef]
Zhou, D.; Thiele-Bruhn, S.; Arenz-Leufen, M.G.; Jacques, D.; Lichtner, P.; Engelhardt, I. Impact ofmanure-related DOM on sulfonamide transport in arable soils. J. Contam. Hydrol. 2016, 192, 118–128.
[CrossRef] [PubMed]
Makselon, J.; Zhou, D.; Engelhardt, I.; Jacques, D.; Klumpp, E. Experimental and Numerical Investigationsof Silver Nanoparticle Transport under Variable Flow and Ionic Strength in Soil. Environ. Sci. Technol. 2017,51, 2096–2104. [CrossRef] [PubMed]
US Congress. Chemicals Used in Hydraulic Fracturing; United States House of Representatives—Committee
on Energy and Commerce: Washington, DC, USA, 2011; p. 29.

Water 2017, 9, 385

77.

78.
79.

80.
81.

82.

83.
84.
85.
86.
87.

88.
89.

90.
91.

92.
93.

94.

95.
96.
97.

32 of 33

Kahrilas, G.A.; Blotevogel, J.; Stewart, P.S.; Borch, T. Biocides in Hydraulic Fracturing Fluids: A Critical
Review of Their Usage, Mobility, Degradation, and Toxicity. Environ. Sci. Technol. 2014, 49, 16–32. [CrossRef][PubMed]
QGC. Approval Conditions 49 to 52: Stage 2 CSG Water Monitoring and Management Plan; 2nd Revision to 23
April 2012 Submission; 13.0 Well Stimulation; QGC: Brisbane, Australia, 2012.
Madsen, T.; Buchardt Boyd, H.; Nylén, D.; Pedersen, A.R.; Petersen, G.I.; Simonsen, F. Environmental and
Health Assessment of Substances in Household Detergents and Cosmetic Detergent Products; Environmental ProjectNo. 615; Danish Environmenal Protection Agency: Copenhagen, Denmark, 2001.
US EPA. OPP Pesticide Ecotoxicity Database. Available online: http://www.ipmcenters.org/ecotox/(accessed on 23 June 2016).
US EPA. Registration Eligibility Decision (RED) Bronopol; United States Environmental Protection Agency,
Office of Prevention, Pesticides and Toxic Substances (EPA-738-R-95-034), U.S. Government Printing Office:Washington, DC, USA, 1995.
Swenberg, J.A.; Kerns, W.D.; Mitchell, R.I.; Gralla, E.J.; Pavkov, K.L. Induction of squamous cell carcinomasof the rat nasal cavity by inhalation exposure to formaldehyde vapor. Cancer Res. 1980, 40, 3398–3402.
[PubMed]Dunnett, P.C.; Telling, G.M. Study of the fate of bronopol and the effects of antioxidants on N-nitrosamine
formation in shampoos and skin creams. Int. J. Cosmet. Sci. 1984, 6, 241–247. [CrossRef] [PubMed]
Challis, B.C.; Yousaf, T.I. Facile formation of N-nitrosamines from bromonitromethane and secondary amines.
J. Chem. Soc. Chem. Commun. 1990, 22, 1598–1599. [CrossRef]
Loeppky, R.N. Nitrosamine and N-Nitroso Compound Chemistry and Biochemistry. In ACS Symposium
Series; American Chemical Society: Washington, DC, USA, 1994.
Douglass, M.L.; Kabacoff, B.L.; Anderson, G.A.; Chent, M.C. The chemistry of nitrosamine formation,
inhibition and destruction. J. Cosmet. Sci. 1978, 29, 581–606.
Cui, N.; Zhang, X.; Xie, Q.; Want, S.; Chen, J.; Huang, L.; Qiao, X.; Li, X.; Cai, X. Toxicity profile of labile
preservative bronopol in water: The role of more persistent and toxic transformation products. Environ. Pollut.
2011, 159, 609–615. [CrossRef] [PubMed]
Scow, K.M.; Johnson, C.R. Effect of sorption on biodegradation of soil pollutants. Adv. Agron. 1997, 58, 1–56.
[CrossRef]
Chiou, C.T. Roles of Organic Matter, Minerals, and Moisture in Sorption of Non-Ionic Compounds andPesticides by Soil. In Humic Substances in Soil and Crop Sciences; Selected Readings; MacCarthy, P., Clapp, C.E.,
Malcolm, R.L., Bloom, P.R., Eds.; SSSA, ASA: Madison, WI, USA, 1991; pp. 111–160.
Delle Site, A.D. Factors Affecting Sorption of Organic Compounds in Natural Sorbent Water Systems andSorption Coefficients for Selected Pollutants. A Review. J. Phys. Chem. Ref. Data 2001, 30, 187–436. [CrossRef]
Ringrose-Voase, A.J.; Young, R.R.; Paydar, Z.; Huth, N.I.; Bernardi, A.L.; Cresswell, H.P.; Keating, B.A.;
Scott, J.F.; Stauffacher, M.; Banks, R.G.; et al. Deep Drainage under Different Land Uses in the Liverpool
Plains Catchment; Report 3, Agricultural Resource Management Report Series; NSW Agriculture: Canberra,
Australia, 2003.
Bennett, S. Key Spatiotemporal Factors Determining Uncertainty of Deep Drainage in a Semi-Arid Area.
Ph.D. Thesis, The University of New South Wales, New South Wales, Australia, 2012.
ANZECC. Australian and New Zealand Guidelines for Fresh and Marine Water Quality; Volume 1: The
Guidelines (Chapters 1–7); Australian and New Zealand Environment and Conservation Council: Canberra,
Australia, 2000.
Stearns, M.; Tindall, J.; Cronin, G.; Friedel, M.; Bergquist, E. Effects of coal-bed methane discharge waters on
the vegetation and soil ecosystem in Powder River Basin, Wyoming. Water Air Soil Pollut. 2005, 168, 33–57.
[CrossRef]
Beletse, Y.G.; Annandale, J.G.; Steyn, J.M.; Hll, I.; Aken, M.E. Can crops be irrigated with sodium bicarbonate
rich CBM deep aquifer water? Theoretical and field evaluation. Ecol. Eng. 2008, 33, 26–36. [CrossRef]
Nghiem, L.D.; Ren, T.; Aziz, N.; Porter, I.; Regmi, G. Treatment of coal seam gas co-produced water forbeneficial use in Australia: A review of best practices. Desal. Water Treat. 2011, 32, 316–323. [CrossRef]
Mace, J.E.; Amrhein, C. Leaching and reclamation of a soil irrigated with moderate SAR waters. Soil Sci. Soc.
Am. J. 2001, 65, 199–204. [CrossRef]

Water 2017, 9, 385

33 of 33

98.

Vance, G.F.; King, L.A.; Ganjegunte, G.K. Soil and plant responses from land application of saline-sodic
waters: Implications of management. J. Environ. Qual. 2008, 37, S139–S148. [CrossRef] [PubMed]
99. Crosbie, R.; Morrow, D.; Cresswell, R.; Leaney, F.; Lamontagne, S.; Lefournour, M. New Insights to the Chemical
and Isotopic Composition of Rainfall across Australia; Water for a Healthy Country Flagship Report Series;
CSIRO: Canberra, Australia, 2012; ISSN 1835-095X.
100. Ringrose-Voase, A.J. Water Balance and Deep Drainage under Irrigated Cotton. In WaterPak; Dugdale, H.,
Harris, G., Neilsen, J., Richards, D., Roth, G., Williams, D., Eds.; Cotton Research and Development Corporation:
New South Wales, Australia, 2004; pp. 17–28.
101. Bureau of Meteorology (BOM). Recent Evapotranspiration. 2014. Available online: http://www.bom.gov.
au/watl/eto/ (accessed on 2 June 2015).
102. Appelo, C.A.J.; Postma, D. Geochemistry, Groundwater, and Pollution, 2nd ed.; Balkema, A.A., Ed.; CRC Press:
Rotterdam, The Netherlands, 2005.
103. McNeal, B.L. Prediction of the effect of mixed-salt solutions on soil hydraulic conductivity. Soil Sci. Soc.
Am. Proc. 1968, 32, 190–193. [CrossRef]
104. McNeal, B.L. Soil Salts and Their Effects on Water Movement. In Drainage for Agriculture; Agronomy No. 17;
van Schilfgaarde, J., Ed.; American Society of Agronomy: Madison, WI, USA, 1974.
105. Biggs, A.; Witheyman, S.L.; Williams, K.M.; Cupples, N.; de Voil, C.A.; Power, R.E.; Stone, B.J. Assessing theSalinity Impacts of Coal Seam Gas Water on landscapes and Surface Streams. In Final Report of Activity 3 of
the Healthy HeadWaters Coal Seam Gas Water Feasibility Study; Department of Natural Resources and Mines:
Toowoomba, Australia, 2012.
106. De Caritat, P.; Lech, M.E. Thomson Region Geochmical Survey, Northwestern New South Wales; CRC LEME Open
File Report 145; CRC LEME: Bentley, WA, USA, 2007.
© 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open accessarticle distributed under the terms and conditions of the Creative Commons Attribution(CC BY) license (http://creativecommons.org/licenses/by/4.0/).