Agricultural Water Management 189 (2017) 52–69

Contents lists available at ScienceDirect

Agricultural Water Management
journal homepage: www.elsevier.com/locate/agwat

Determining water quality requirements of coal seam gas produced
water for sustainable irrigation
Dirk Mallants a,∗ , Jirka Šimůnek b , Saeed Torkzaban aab

CSIRO Land and Water, Waite Road – Gate 4, Glen Osmond, SA 5064, Australia
University of California, Riverside, CA 92521, USAa r t i c l e

i n f oArticle history:
Received 2 March 2017
Received in revised form 18 April 2017
Accepted 20 April 2017
Keywords:
Soil management
Salinity risk
Coupled processes
Major ion chemistry
HYDRUSa b s t r a c t
Coal seam gas production in Australia generates large volumes of produced water that is generally highin total dissolved solids and has a high sodium absorption ratio (SAR) which may affect soil structure,
hydraulic conductivity, and crop production if used untreated for irrigation. By coupling major ion soil
chemistry and unsaturated ﬂow and plant water uptake, this study incorporates effects of salt concentrations on soil hydraulic properties and on root water uptake for soils irrigated with produced waterfeaturing different water qualities. Simulations provided detailed results regarding chemical indicators ofsoil and plant health, i.e. SAR, EC and sodium concentrations. Results from a base scenario indicated thatthe use of untreated produced water for irrigation would cause SAR and EC values to signiﬁcantly exceed
the soil quality guide values in Australia and New Zealand (ANZECC). The simulations provided further
useful insights in the type of coupled processes that might occur, and what the potential impacts could
be on soil hydrology and crop growth. Calculations showed that the use of untreated produced water
resulted in a decrease in soil hydraulic conductivity due to clay swelling causing water stagnation, additional plant-water stress and a reduction in plant transpiration. In case the produced water was mixed
with surface water in a 1:3 ratio prior to irrigation, the calculated soil SAR values were much lower andgenerally acceptable for sandy to sandy-loam soil. The use of reverse osmosis treated produced water
yielded an acceptable salinity proﬁle not exceeding guide values for SAR and EC; the plant water stress
was limited as there was no additional salinity stress associated with the low level of salts. Results further illustrated that accounting for coupled geochemical, hydrological and plant water uptake processes
resulted in more accurate water balance calculations compared to an approach where such interactions
were not implemented. Coupling unsaturated ﬂow modelling with major ion chemistry solute transport using HYDRUS provides quantitative evidence to determine suitable water quality requirements forsustainable irrigation using coal seam gas produced water.
© 2017 Commonwealth Scientiﬁc and Industrial Research Organisation. Published by Elsevier B.V.
This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/).

1. Introduction
Production of coal seam gas in Australia requires depressurisation of the coal seam layers by extracting large volumesof groundwater. Reducing the hydrostatic pressure on the gas in
the pores of coal seams allows gas to ﬂow through natural andenhanced fracture networks to the gas production well (Moore,2012). Water extraction rates in the main coal seam gas areas inAustralia vary over time and geographically (QWC, 2012). Basedon a 50-year time frame that is generally assumed for forecast-∗ Corresponding author.

ing water production (DNRM, 2012), average water production
rates range from 0.028 GL/year in the Clarence Moreton basin (NewSouth Wales) (RPS, 2011) to 75–98 GL/year in the Surat Basin(Queensland) (QWC, 2012). By 2050, most Australian coal seam gas
areas will have reached the end of their production stage (DNRM,
2012).
Treated coal seam gas water may be utilised for beneﬁcial uses,
including agriculture (irrigation and stock), urban uses (town watersupply) and industrial uses (construction and processing). Assessment of the beneﬁcial use of treated coal seam gas water for theCentral Condamine Alluvium (Queensland) indicates the region has
the capacity to deliver 854 GL (gigalitres), or around 35% of thehistoric depletion, by 2050 to irrigators and the Chinchilla Weir
water supply scheme (DNRM, 2013). Other irrigation projects in

E-mail address: Dirk.Mallants@csiro.au (D. Mallants).
http://dx.doi.org/10.1016/j.agwat.2017.04.011
0378-3774/© 2017 Commonwealth Scientiﬁc and Industrial Research Organisation. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

Australia where produced water is used include large-scale forestry
in Queensland’s Bowen Basin and legume plantations (RPS, 2011).
In the Bowen Basin, up to 8 ML/day of reverse osmosis (RO) treatedproduced water was used to irrigate 234 ha of pasture crop andan 800,000 tree timber plantation (Santos, 2009, 2010). Also in
Queensland, 300 ha of oilseed-bearing legume tree (Pongamia pinnata) plantation are irrigated with RO treated produced water(Parsons, 2010).
Produced water is generally unsuitable for direct surface discharge or irrigation without any treatment or amendment (Stearns
et al., 2005; Beletse et al., 2008; Nghiem et al., 2011). While sodium
ions cause soil particles to disperse, particularly if the soils containmontmorillonite clays, most ions increase the aggregation of soil
particles (Nghiem et al., 2011). Irrigation water with a high SAR1
can lead to a decrease in inﬁltration and deterioration of the soil
structure as the dispersed clay minerals, once dry, cause soils to
become dense, cloddy and structureless, destroying natural particle aggregation. SAR values greater than 13 pose a risk to the soil
ecosystem (Stearns et al., 2005), and even SAR values between 5 and8 have been shown to cause irreversible plugging of soil pores andswelling (Mace and Amrhein, 2001). Traditional treatments to mitigate saline-sodic irrigation water can be used, such as for example
the addition of gypsum and elemental sulphur (Vance et al., 2008;
Šimůnek et al., 2006).
Studies have been undertaken to assess the feasibility of using
sodium-rich produced water for salt-tolerant crop production(Johnston et al., 2008; Vance et al., 2008; Beletse et al., 2008).
Vance et al. (2008) used saline-sodic coal seam gas produced waterwith an SAR between 17 and 57 for irrigating grasslands and hayﬁelds in the Powder River Basin, Wyoming. For the use of produced
water to be effective, the method further required applicationsof gypsum and elemental S to provide calcium ions and an acidiﬁed soil environment to promote calcite dissolution. Further ﬁeld
studies in Wyoming with bioenergy feedstock species indicatedthat produced water can be used for short-period (2 years) irrigation of such crops (Burkhardt et al., 2015). However, the study
also concluded that prolonged use of untreated produced water
for irrigation would likely have deleterious long-term effects on
the soil and plants unless the water was treated or diluted with
good-quality water.
Water produced during the depressurisation phase of coal seamgas mining in Waterberg district, South Africa, is highly salineand dominated by sodium bicarbonate. With careful management,
Beletse et al. (2008) determined that coal seam gas irrigation waterwith SAR values of 85 could be used to grow certain crops, but
additions of gypsum and organic matter to the soil were necessary
to counteract inﬁltration problems that arose due to the excessive
sodium that had accumulated in the soil.
Investigations of soil and vegetation recovery in British
Columbia document the natural recovery of salt-affected plots(Leskiw et al., 2012), with natural attenuation from above normal rainfall in the ﬁrst few years after the salt deposition being
effective in removing salts from the root zone and subsoil. These
results, however, were obtained for a humid climate; it is expectedthat natural attenuation under arid or semi-arid climates would bemuch less effective. Bright and Addison (2002a,b) developed standardised guidelines to assist spill response and soil remediation innorthern British Colombia at sites where salt-containing produced
water is released as part of oil and gas exploration and extraction activities. Generic soil quality standards for human health,aquatic life, and soil ecological functioning, as manifested throughsoil invertebrate and plant responses, were developed for salt ions.

1
Sodium adsorption ratio (SAR) = [Na+ ]/{([Ca2+ ] + [Mg2+ ])/2}1/2 where [Na+ ],[Ca2+ ], and [Mg2+ ] are in meq/L.

53

Whilst the impairment of metabolic functioning of soil microbes
is recognised as playing a major role in nutrient cycling and otherprocesses important to terrestrial ecosystems, there is insufﬁcient
data available to deﬁne a threshold for salt ions; instead, a microbial functional impairment standard is used that is developed fromnutrient and energy cycling data.
In Australia, the Queensland Department of Environment andHeritage Protection (DEHP) recommends the minimum standards
for using produced water for irrigation purposes. The guidelines
indicate that the electrical conductivity (EC) should be less than
950 ␮S/cm or 0.95 dS/m (95th percentile over a one-year period),
and the SAR should be less than 6 for heavy soils and less than
12 for light soils (95th percentile over a one-year period) (DEHP,
2014). Another potential issue when produced water is used for
irrigation is the relatively high boron content. If not removed, boron
can be damaging to some ﬂora. The boron ANZECC guide value forlong-term irrigation (up to 100 years) is 0.5 mg/L.
In Australia and New Zealand, the minimum quality of irrigationwater is deﬁned in The Australian and New Zealand Guidelines forFresh and Marine Water Quality (ANZECC, 2000). Plant salt tolerances and recommended trigger values for metals and metalloids in
irrigation water are also provided in ANZECC (2000). For instance,for different soil types, permissible levels of chloride, sodium, andSAR have been deﬁned (Table 1).
To support the sustainable management of soil and landscapes
under irrigated agriculture, Biggs et al. (2012) developed a framework to assess the salinity risk associated with the use of coal seamgas water for irrigation in the Queensland Murray-Darling Basin.
Biggs et al. (2012) identiﬁed the description of the unsaturated zone
as the key knowledge and data gap when conducting salinity risk
assessments. In particular, studies are needed about (i) critical soil
attributes (porosity, initial soil water content, substrate hydraulicconductivity) to depths below the “root zone”, (ii) the extent of
the root zone for different crops, and (iii) the quantiﬁcation oflateral ﬂow processes. Biggs et al. (2012) further highlighted theneed to better deﬁne critical thresholds for agronomic and degradation purposes (e.g. acceptable root zone salinity) across landscapes,
improve process understanding such as the fate of excess water(lateral ﬂow or drainage), and deeper soil/regolith processes andproperties.
The current study endeavours to address those knowledge
gaps that relate to better process understanding (e.g., critical soil
properties, critical thresholds) using the unique features of the
HYDRUS-1D simulator that are relevant for salinity risk assessment (Jacques et al., 2013a; Šimůnek et al., 2006, 2016). To this end,
this study simulates coupled processes of variably saturated water
ﬂow, plant water uptake and coupled transport of multiple major
ions in soils irrigated with produced water of different water qualities. Through this process coupling, we explicitly include criticalsoil processes required for salinity risk assessment associated withcoal seam gas produced water in the analysis. Simulations of the
movement of multiple ions in a vegetated soil proﬁle under irrigation are based on the major ion chemistry module UnsatChem(Šimůnek and Suarez, 1994) implemented in the ﬁnite element
code HYDRUS-1D (Šimůnek et al., 2008).

2. Materials and methods
2.1. Soil hydrological model
Simulation of variably saturated ﬂow in soil requires a mathematical relationship between (i) the soil water content () and thesoil pressure head (h), i.e. the soil water retention curve (h), and (ii)either the water content or the pressure head and the unsaturated
hydraulic conductivity K() and K(h), respectively. We applied the

54

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

Table 1
Permissible levels of irrigation water quality parameters to preserve soil and plant health.
Parameter

Maximum (DEHP, 2014)

Electrical conductivity (␮S/cm)k
SAR
Bicarbonate (mg/L)
Fluoride (mg/L)
pH
Chloride (mg/L)
Sodium (mg/L)
Heavy metals and metalloids (trigger values in
mg/L; short-term use up to 20 years)<950i<6i for heavy soils; <12i for light soils
100j2
6–8.5iANZECC (2000)<650>20a ; 20–8b ; 13–6c ; 11–5d
No trigger value recommended
1
6–9<175e ; 175–350f ; 350–700g ; > 700h
<115e ; 115–230f ; 230–460g ; >460h
Al (20); As (2); Cd (0.05); Cr (1); Co (0.1); Cu (5); Fe (10); Li (2.5); Pb (5); Mn (10);Hg (0.002); Mb (0.05); Ni (2); Zn (5)aSand-sandy loam.
Loam-silty loam.
c
Clay loam.
dLight clay.
eSensitive crop.
fModerately sensitive crop.
g
Moderately tolerant crop.
hTolerant crop.
i
95th percentile over a one-year period.
j
DERM (2010).
k1000 ␮S/cm = 1 dS/m.
b

analytical model of van Genuchten, i.e., nonhysteretic, to describe
the soil water retention curve, (h) (van Genuchten, 1980), since it
permits a relatively good description of (h) for many soils using
only a limited number of parameters. The van Genuchten soil moisture retention characteristic is deﬁned as:
(h) = r +s − r(1 + |ah|n )(1)mwhere  r is the residual water content [cm3 /cm3 ],  s is the saturated water content [cm3 /cm3 ], and ˛ [1/m], n [–] and m (=1 − 1/n)[–] are the retention curve shape parameters. The van Genuchten(1980) (h) equation has the additional advantage that when coupled with the unsaturated hydraulic conductivity K(h) model of
Mualem (1976), it produces the following closed-form expression:
1/m m

K(h) = Ks Sel [1 − (1 − Se

) ]

2

(2)

where Ks is saturated hydraulic conductivity [cm/day],Se = ( −  r )/( s −  r ) is the effective saturation [–], m = 1 − 1/n,
and n > 1. The pore-connectivity parameter (l) was estimated by
Mualem (1976) to be about 0.5 as an average for many soils. In thisstudy, we use Eq. (2) with the l parameter ﬁxed at 0.5.
2.2. Soil geochemical model
Migration of chemicals in soil is primarily by advection
and dispersion (ignoring gaseous transport) when advection isgreater than the rate of diffusion. For the transport of contaminants during transient water ﬂow, HYDRUS implements theadvection–dispersion equation (Šimůnek et al., 2008). In this study,
a dispersivity of 0.1 m was taken, or one tenth of the total travel
distance from top to bottom of the 1-m deep soil proﬁle.
The solute transport module of HYDRUS-1D as described above
is limited to single ions, or ions subject to relatively simple
ﬁrst-order consecutive decay or transformation reactions (e.g.,nitriﬁcation–denitriﬁcation chains, or radionuclide decay chains).
As an extension to this approach, HYDRUS-1D implements a major
ion chemistry module based on the UnsatChem model (Šimůnek
and Suarez, 1994, 1997; Šimůnek et al., 2008). The UnsatChem module of HYDRUS-1D enables quantitative predictions of processes
involving major ions, such as simulations of the effects of salinity
on plant growth and estimating the amount of water and amendments required to reclaim soil proﬁles to desired levels of salinity
and ESP (exchangeable sodium percentage) (Šimůnek et al., 2006).

The UnsatChem module (Šimůnek and Suarez, 1994) considers the transport of major ions and carbon dioxide in soils. The
major variables of the chemical system are Ca, Mg, Na, K, SO4 ,Cl, NO3 , H4 SiO4 , alkalinity, and CO2 (Supplementary Material
Table S1). The model accounts for equilibrium chemical reactions
between these components such as complexation, cation exchange,and precipitation–dissolution. For the precipitation–dissolution of
calcite and the dissolution of dolomite, either equilibrium or multicomponent kinetic expressions are used, including both forward
and backward reactions. Other precipitation–dissolution reactions considered involve gypsum (CaSO4 ·2H2 O), hydromagnesite(Mg5 (CO3 )4 (OH)2 ·4H2 O), nesquehonite (MgCO3 ·3H2 O) (Mg), andsepiolite (Mg2 Si3 O7.5 (OH)·3H2 O). Since the ionic strength of soil
solutions can vary considerably in time and space and often reachhigh values, both the modiﬁed Debye–Hückel and Pitzer expressions were incorporated into the model, thus providing options forcalculating single-ion activities (Šimůnek et al., 2013).
Partitioning of dissolved major ions between the solid andsolution phases is described with the Gapon equation (Whiteand Zelazny, 1986) provided in the major ion chemistry module.
This requires the deﬁnition of the Gapon Exchange Constants forexchange of calcium and magnesium, calcium and potassium, andcalcium and sodium. The Gapon Exchange Constants are takenfrom Jacques et al. (2012) and are interpreted as follows: for theCa2+ /Mg2+ exchange (KCa/Mg = 1.2), Ca2+ is the preferred ion relative to Mg2+ on the exchange complex. For the Ca2+ /Na+ exchange(KCa/Na = 2.9), Ca2+ is the preferred ion relative to Na+ . Finally, for
the Ca2+ /K+ exchange (KCa/K = 0.2), K+ is the preferred ion relative
to Ca2+ .
UnsatChem also considers the effects of solution compositionon the unsaturated soil hydraulic properties. The accumulation ofmonovalent cations, such as sodium and potassium, often leads toclay dispersion, swelling, ﬂocculation and overall poor soil physicomechanical properties. These processes have an adverse effect on
the soil hydraulic properties including hydraulic conductivity, inﬁltration rates and soil water retention as a result of swelling and clay
dispersion. These negative effects are usually explained based on
the diffuse double layer theory.
The effect of solution chemistry on the hydraulic conductivity
in the major ion chemistry module is calculated as:
K(h, pH, SAR, C0 ) = r(pH, SAR, C0 )K(h)(3)

where SAR is the sodium adsorption ratio, C0 is the total salt concentration of the ambient solution (meq/L), and r is a scaling factor

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

which represents the effect of the solution composition on the ﬁnal
hydraulic conductivity [–], and which is related to pH, SAR, andsalinity (Šimůnek et al., 2008). Conceptually, this effect is related
to clay swelling as a result of changes in ion composition here
expressed via the parameters SAR and C0 (McNeal, 1968, 1974).
Although soil water retention curves may also be affected by SAR
and electrolyte concentrations (e.g., Lenhard and Brooks, 1986), thisis not being considered in the current simulations.
We note that other effects of solution chemistry on hydraulicproperties such as mineral precipitation resulting in a decreaseof pore space and hence a decrease in hydraulic conductivity
may be equally important to determine a soil’s water balance
under irrigation in arid and semi-arid regions. Examples of amore general approach to couple major ion chemistry with porous
media hydraulic and solute transport parameters are available fromJacques et al. (2013a,b).
2.3. Soil hydrological data
The soil proﬁle used in this study was assumed to be 1 m
deep, with three soil layers: 0–0.3 m, 0.3–0.7 m, and 0.7–1.0 m. The
van Genuchten soil hydraulic properties were taken from Bennett(2012). These soil hydraulic properties were calculated on the basisof the pedotransfer functions of Minasny and McBratney (2002)
using soil data points within the Cox’s Creek catchment, NewSouth Wales. For three soil horizons, the particle size distribution and the bulk density from the SALIS (New South Wales Soil
And Land Information System) dataset were used as predictor variables in the pedotransfer functions. From six soil types (Tenosols,
Dermosols, Vertosols, Chromosols, Sodosols, Kandosols), Vertosols
were most prominent in 84 out of 143 soil proﬁles (59%). Vertosol
hydraulic properties were therefore selected here for the simulations (Supplementary Material Table S2). The mean soil water
retention curves together with their uncertainty bands (mean ± 2
standard deviations) are shown in Supplementary Material Figure
S1.
In order to reach the (pseudo) steady-state concentrations proﬁle, the following hypothetical irrigation scheme was used. The
potential ET is equal to 1 cm/day (based on the average from January
2014 in Narrabri, BOM, 2014), and is divided into potential evaporation (2 mm/day) and potential transpiration (8 mm/day). Over a
5-day period 6 cm of irrigation was applied on the ﬁrst day followedby 4 days of not irrigating; this resulted in a total of 73 irrigation
days out of a simulation period of 365 days. For each 5-day period,
this yielded an irrigation excess of 20% above the potential ET foreach day. The annual total irrigation of 438 cm is much higher thanthe combined rainfall and irrigation in irrigation areas of the Northern Darling Basin for crops such as cotton (up to nearly 140 cm) orsorghum (up to 110 cm) (Silburn and Montgomery, 2004). However, irrigation of pasture, especially for grazing (beef and lamb),
uses more water: for example, nearly 150 cm was recorded in aﬂood irrigated pasture in the South East of South Australia (Horizon,
2006). The application rate used in this study, therefore, corresponds to the equivalent applied water of two to three growing
seasons. Our purpose was not to have the most realistic irrigation regime, rather, we applied a simpliﬁed water application ratethat would simplify the analysis of the complex coupled chemicalhydrological processes.
Zero precipitation was considered for the modelling; while thisis a simpliﬁcation of a soil’s hydrological balance and does notaccount for effects of rainfall salinity on the salt mass balance,
it does allow to focus on effects of irrigation only when coupled chemical-hydrological-biological processes are considered.
The root zone was assumed to be 0.50 m deep and root distribution was assumed to decrease linearly with depth. The initialpressure head at the beginning of each simulation was assumed

55

to be -1 m throughout the soil proﬁle. The top boundary condition
was speciﬁed as an atmospheric boundary with surface runoff, i.e.
whenever rainfall intensity exceeds the soil’s inﬁltration capacity
ponding does not occur and all excess water is considered runoff
and removed from the soil surface. At the bottom of the proﬁle, afree drainage BC was imposed. The soil water and chemical massbalances were calculated for one year.
2.4. Crop data
Two cases were considered to simulate stress effects: (i) water
stress without salinity stress and (ii) both water and salinity
stresses on crop production. First, the effect of water stress was
separately analysed using the Feddes’ model (Feddes et al., 1978).
A conceptual representation of the suction-dependent plant-water
stress response function is shown in Supplementary Material Figure S2 and relevant parameters based on Wesseling et al. (1991) in
Table S3.
When salinity stress is considered in HYDRUS, a selection must
be made whether the effect of salinity stress is additive or multiplicative to water stress. Here the multiplicative Threshold Modelaccording to Maas (1990) was selected. The Threshold Model has
two parameters, the threshold parameter and the slope parameter(Table S3). The ﬁrst parameter represents the value of the minimum osmotic head (the salinity threshold) above which root water
uptake occurs without any reduction. The second parameter is theslope of the curve determining the fractional root water uptake
decline per unit increase in salinity below the threshold. Root water
uptake reduction parameters for alfalfa were used to investigate
combined water and salinity stress as no salinity stress parameters
are available for pasture (all other root water uptake parameters
are those for pasture).
The volume of water removed from a unit volume of soil per unit
time due to plant water uptake, S(h) [day−1 ], is deﬁned as (Šimůnek
et al., 2008)
S(h) = ˛(h)b(z)Tp(4)

where ˛(h) is the plant-water stress response function (between
0 and 1, see Supplementary Material Figure S2), b(z) is thenormalised water uptake distribution [cm−1 ] (here assumed linearly decreasing with depth), and Tp is potential transpiration[cm/day]. With Tp = 0.8 cm/day, and b(z = 0) = 0.04 cm−1 (whenrooting depth = 50 cm; 2 × 1/50), S(h) at the soil surface equals
0.032 day−1 .
2.5. Chemical composition of irrigation water
Three different chemical compositions of irrigation water wereassumed (Table 2); the ﬁrst corresponds to untreated produced
water, the second is a 3:1 surface-water amended produced water,
and the third is RO-treated produced water. The ﬁrst (untreated)
and third (treated) compositions represent end-members that canbe mixed to various degrees to obtain ﬁt-for-purpose irrigationwater. Unlike the amended and treated produced waters thatare actually being ﬁeld-tested and/or being used in commercial
irrigation schemes, the untreated produced water would almost
certainly not receive any permission for use as irrigation waterowing to its high potassium and SAR values (Table 2). It is included
here only to demonstrate effects of several interacting processes,
such as major ion composition triggering a reduction in soil
hydraulic conductivity and subsequent modiﬁcation of water redistribution in the soil proﬁle and how root water uptake responds
to such changes. It is worth noting that the treated water is stillmuch more saline than rainfall. In other words, even irrigation withtreated produced water is likely to have some effect on the soil’s
chemical balance.

56

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

Table 2Composition of untreated, amended, and treated produced waters from Australian coal basins. Treated/rain water = ratio of treated produced water to rain water composition.
Alkalinity

Cl

SO4

Ca

K

Mg

Na

SAR

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

1406
61.1
3064
53.7

77
Total EC = 5705
100%

Amended produced water (3:1 ratio surface water:produced water)
183
164
mg/L
4.5
4.6
meq/L
ECi a (␮S/cm)
134
353
8.5
22.5
%ECi a

337
7.0
561
35.7

10.30
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

11
Total EC = 1573
100%

Treated produced water
mg/L
meq/L
ECi a (␮S/cm)%ECi aTreated/rain water

1
0.02
1.67
0.6
3

10.33
0.52
30.7
10.5
43

0.43
0.011
0.81
0.3
4

7.6
0.62
33
11.3
95

28.7
1.25
62.5
21.5
60

119

73.5

106

50.1

Untreated produced water
mg/L
meq/L
ECi a (␮S/cm)%ECi a

25
0.61
18.2
6.3
417

67
1.89
144.3
49.6
103

Electrical conductivity for ions as trace concentration in water at 25 ◦ C
44.5
76.35
160
S cm2 /mola

1.6
Total EC = 291
100%

The electrical conductivity contribution of each ion (%ECi ) to the total EC of the solution is calculated as (Appelo and Postma, 2004):

Table 3
Initial concentrations of exchangeable cations (Source: Ringrose-Voase, 2004) andinitial soil solution composition used to calculate initial solution in equilibrium withthe cation exchange complex. Data represents composition of rainwater in WaggaWagga (Source: Crosbie et al., 2012). Gapon Exchange Constants (Source: Jacqueset al., 2012).
Soil layer

CEC (meq/kg soil)

1 (0–0.3 m)
2 (0.3–0.7 m)
3 (0.7–1.0 m)

275
278
274

The initial solution composition in the soil proﬁle (Table 3) was
deﬁned using the solution composition of rainwater measured in
Wagga-Wagga (Crosbie et al., 2012). Note that this solution will be
automatically brought into equilibrium with the cation exchange
complex.

Gapon Exchange Constants
K [Ca/Mg]

K [Ca/Na]K [Ca/K]

1.2
1.2
1.2

2.9
2.9
2.9

0.2
0.2
0.2

Exchangeable cations for layer 1
Ions

Ca

Mg

Na

K

mmol/kg
meq/kg

86.3
173

44
88

4.6
4.6

9.4
9.4

Initial solution composition (meq/L)
Ca

Mg

Na

K

Alkalinity

SO4

Cl

0.012

0.00660.021

0.00280.00098

0.00690.018

It is important to note that from all the ions considered in thesimulation, only potassium (K) has an ANZECC guide value for irrigation (ANZECC, 2000). For long-term use (up to 100 years), theguide value for K is 0.05 mg/L. For short term use (up to 20 years),
the guide value is 0.8 mg/L. The calculated SAR for the untreatedproduced water is 77 (Table 2); this is considerably higher than the
DEHP maximum guide value of 12 for light soils (95th percentile
over a one-year period) (DEHP, 2014).
2.6. Initial soil chemical composition
The initial concentrations of soil exchangeable cations (Table 3)
were taken from a lysimeter study by Ringrose-Voase (2004).
The lysimeter is located in paddock C1 at the Australian Cotton
Research Institute (ACRI) near Narrabri in northern New South
Wales (30◦ 11.53 South, 149◦ 36.31 East). The soil is a Haplic, selfmulching, Grey Vertosol. The bulk density for different soil depths
was estimated to range between 1.32 and 1.45 g/cm3 (RingroseVoase, 2004). Because the total cation exchange capacity is similar
for all three soil layers, the cation exchange capacity and the concentration of exchangeable cations of the ﬁrst layer was assigned
to all three layers.

2.7. Simulation scenarios
Scenario 1: A base scenario against which other scenarios canbe compared. The base scenario involves soil water balance andgeneral solute transport calculations without (i) effects of major
ion chemistry on hydraulic properties, (ii) effects of salinity stresson root water uptake, and (iii) effects of competing adsorption
between major ions on ion migration in soil. Root water uptake
without solute stress is calculated for pasture and alfalfa. The irrigation water is solute-free and only contains a single tracer, assuming
a hypothetical unit tracer concentration of 1 mg/L, used to evaluate
solute migration based on a 6/5 cm/day irrigation rate (averaged
across one day of irrigation at 6 cm/day followed by 4 days of zero
irrigation). The standard Hydrus “General Solute Transport” module was used to model salinity as a non-reactive tracer with a unit
input concentration.
Scenario 2: Water balance calculations account for root water
uptake with a multiplicative solute stress. Major ion chemistry
solute transport is invoked excluding the effects of a reduction
in hydraulic conductivity. This scenario uses untreated produced
water at a 6 cm/day irrigation rate (for chemical composition see
Table 3). This scenario uses the Hydrus “Major Ion Chemistry”(UnsatChem) module to calculate concentrations of major ions (Ca,Mg, Na, K, alkalinity, Cl, SO4 ) present in irrigation water.
Scenario 3: The main difference with Scenario 2 is that major
ion chemistry solute transport now includes calculation of reduction in hydraulic conductivity due to soil solution chemistry (basedon Eq. (3)). Untreated produced water is used at a lower 4 cm/day
irrigation rate followed by 4 days of no irrigation to avoid water
ponding due to reduced hydraulic conductivity (calculations with
6 cm/day irrigation rate produced ponding, results not shown).
In addition to the processes considered in Scenario 2, the reduction of hydraulic conductivity due to the solution composition is
accounted for (McNeal, 1968, 1974). Scenario 4: Same as Scenario
3, but now surface water-amended produced water (at a 3:1 surface water:produced water ratio) is used at a 6 cm/day irrigation
rate (for chemical composition see Table 2).

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

57

Fig. 1. (A) Simulated solute concentrations following inﬁltration of solute-containing irrigation water (Scenario 1). Comparison of models including ET (solid lines), including
E only (evaporation, dashed lines) or excluding ET (evapotranspiration, symbols). (B) Comparison of simulated soil water contents for models including ET, excluding T, or
excluding ET.

Scenario 5: Same as Scenario 4, but here treated produced water
is used at a 6 cm/day irrigation rate (for chemical composition see
Table 2).

3. Results and discussion
3.1. Base scenario simulation
Simulated water content and pressure head for the base scenario (Supplementary Material Figure S3) time series show typical
patterns linked to alternating periods of soil wetting during andfollowing irrigation and drying periods in the intermittent periods between irrigation events. At deeper depths in the soil proﬁle,the variations become smaller as the water redistribution becomesincreasingly buffered as it moves through the soil.
The solute accumulation in the soil proﬁle at different depths isshown in Fig. 1A. Three different simulation results are shown toillustrate the effect evaporation and transpiration processes haveon solute concentration.
First, consider the solute concentration proﬁles at increasing
times since the start of the irrigation for simulations that include
evaporation and transpiration (Fig. 1A). Proﬁles are obtained at
the end of each irrigation-free four-day period, just before irrigation starts (see output time steps on Supplementary Material
Figure S2). During this four-day period, as a result of evaporation, water content at the soil surface decreases (Fig. 1B), creating
an upward pressure head gradient, which produces upward waterﬂow towards the soil surface. Water ﬂowing towards the soil surface carries with it solute ions, which accumulate at the soil surface
during evaporation. Note that during evaporation, water leaves thesoil surface while leaving solute ions behind. As a result of decreasing water contents and upward water ﬂow towards the surface(bringing additional solute ions), the concentration of solute ionsincreases signiﬁcantly during this time period (visible in the top
0.1 m of the soil proﬁle, Fig. 1A). Transpiration will further reduce
the water content across the root zone (Fig. 1B), which further
increases solute concentrations.

As a whole, the solute proﬁles reﬂect the equilibrium between
inﬁltration with irrigation water, water and solute redistribution
due to transpiration and evaporation, and downward leaching(Fig. 1A ). The prominent effects of transpiration and evaporation on the solute proﬁle become further apparent when theseprocesses are suppressed (Fig. 1A). For instance, when transpiration T is suppressed during the simulations a much smaller soluteconcentrations across the entire soil proﬁle become apparent. If
subsequently also soil evaporation E is suppressed, the effect ismainly a further concentration reduction near the soil surface while
deeper in the soil the proﬁles that exclude T and ET (soil evaporation and plant transpiration) are similar at 30 days and only slightly
different at 365 days. These effects are mainly caused by differencesin soil water contents between the different models (with or without T and ET): the highest water content is observed when bothtranspiration and evaporation are suppressed, whereas the lowest
water contents exist when those processes are included (Fig. 1B).
As noted, the effect of suppressing transpiration is much greaterthan suppressing evaporation (daily potential transpiration rate is
0.8 cm/day while daily potential evaporation rate is 0.2 cm/day). At
the top of the soil proﬁle, the water content without ET is 1.7 times
larger than when ET is included. This ratio is close to the ratio of
2.2 obtained by dividing solute concentration with ET and without
ET at the same depth.
Cumulative water ﬂuxes across model boundaries (soil surface
and bottom of the model domain) were calculated for pasture andalfalfa (note that only the water stress response function is applied
for both crops). As can be observed from Fig. 2A and Table 4, there
is no noticeable difference between the two crops. For both crops,
the potential (292 cm) and actual (281 cm) annual cumulative transpiration ﬂuxes are very similar, indicating there was little water
stress during the simulations. The cumulative potential and actual
surface ﬂuxes (considering inﬁltration and evaporation) are both
negative (downward) and identical (−365 cm), meaning potential
and actual evaporation are equal (73 cm). The surface ﬂux is basedon an annual total irrigation of −438 cm (negative in the water balance calculation). The cumulative actual bottom ﬂux (i.e. drainage)

58

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

Fig. 2. Cumulative potential (black) and actual (purple) surface ﬂuxes (inﬁltration minus evaporation), cumulative potential (red) and actual (blue) transpiration ﬂuxes, andcumulative actual (green) bottom ﬂux. Cumulative potential and actual surface ﬂuxes are identical across all scenarios. (A) Scenario 1 (pasture and alfalfa). (B) Scenario 2(pasture). (C) Scenario 3 (pasture). (D) Scenario 4 (pasture). (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of thearticle.)Table 4
Water balance components for a 365 day simulation period. Negative inﬁltration values means inﬂow; negative net drainage values means outﬂow.
Scenario

Cumulative inﬁltration(cm)

Cumulative actual
transpiration (cm)Cumulative actual
evaporation (cm)Cumulative net
drainage (cm)

1: Pasture – tracer only–ADE
1: Alfalfa – tracer only–ADE
2: Pasture – UPW–MIC–WSS
3: Pasture – UPW–MIC–WSS–KRM
4: Pasture – SAPW–MIC–WSS–KRM
5: Pasture – TPW–MIC–WSS–KRM

−438a
−438a
−438a
−292b
−438a
−438a

281 [96%]281 [96%]240 [82%]
188 [64%]281 [96%]
280 [96%]

73 [100%]
73 [100%]
73 [100%]
69 [95%]73 [100%]
73 [100%]

−84
−83
−123
−29
−84
−84

ADE: standard advection–dispersion model; UPW: untreated produced water; MIC: major ion chemistry combined with ADE; WSS: water and solute stress model; KRM:
hydraulic conductivity reduction model; SAPW: surface water amended produced water; TPW: treated produced water.
a0.06 cm/day once every 5 days.
b
0.04 cm/day once every 5 days.

is negative, i.e. −84 cm for pasture and −83 cm for alfalfa. Thisdrainage ﬂux is expected as the irrigation rate exceeds the potential ET by a factor of 6 during day one of the irrigation cycle (6 cm
irrigation with 1 cm potential ET).

3.2. Major ion chemistry
We now consider major ion chemistry to simulate behaviour ofmultiple solutes added to the soil proﬁle through the use of pro-

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

59

Fig. 3. Solute tracer concentration proﬁles for Scenarios 2 (A), 3 (B), 4 (C) and 5 (D). Dashed lines show corresponding solute tracer concentration proﬁles for Scenario 1.

duced water; Scenario 2 with untreated water is discussed ﬁrst. This
scenario provides insight into the distribution of individual ions(both in the aqueous phase and adsorbed on the solid phase), theEC and the SAR in soil. The simulations also provide concentrationsof any minerals such as calcite, dolomite or gypsum, should they
form (see Supplementary Material Table S1 for a full list of minerals
simulated). Scenario 2 further links the soil chemistry (i.e. salinity) to the salinity stress of plants (in addition to the water stressinvoked in Scenario 1), allowing for calculations of the reduction
in water uptake by vegetation due to differences in water quality.
The reduction in water uptake, in turn, has an effect on the waterand solute distribution in the soil proﬁle, and thus on the drainage
ﬂux of water and solutes across the bottom of the proﬁle. Unlike
Scenarios 3, 4, and 5, the effects of salinity on hydraulic properties
were not considered in Scenario 2.
Due to the combined effects of water and salinity stresses on root
water uptake, which will be explained in the following paragraph,
actual transpiration for Scenario 2 (240 cm) is considerably lower
than potential transpiration (292 cm) (Fig. 2 and Table 4). Lowertranspiration results in more water remaining in the soil (Fig. 4A),

and more water draining through the bottom of the proﬁle (123 cm
for Scenario 2 versus 84 cm for Scenario 1, Table 4).
The tracer proﬁle developed over one-year simulation for Scenario 2 is quite different from Scenario 1 (Fig. 3A). If only water
stress is accounted for in the absence of major ion chemistry (Scenario 1), the tracer concentration gradually builds up over the entire
soil proﬁle, with a sharp concentration gradient near the soil surface(Fig. 3A). Also note that the output is obtained just before irrigation starts (at the end of a period with irrigation at day 1 followedby 4 days of ET), when water contents in the top layer are lowest(and consequently concentration are the highest). When major ion
chemistry and the salinity stress are accounted for (Scenario 2), thesalinity stress results in a signiﬁcant reduction in the root water
uptake and hence a reduction in actual transpiration, i.e. from 281(Scenario 1) to 240 cm (Fig. 2B). This, in turn, mainly changes thesolute concentration distributions whereas water content is less
impacted as exempliﬁed in Figs. 3A and 4A , respectively: the lower
concentrations in the bottom part of the soil proﬁle are due to thehigher drainage ﬂux (−123 cm for Scenario 2 compared to −84 cm
for Scenario 1) causing dilution of the tracer.

60

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

Fig. 4. Water content proﬁles for Scenarios 2 (A), 3 (B), 4 (C) and 5 (D) and their comparison with the water content proﬁle for Scenario 1 at 365 days.

Water content at a particular depth in the soil proﬁle seems
quite constant during the simulation (except for Scenario 3, Fig. 4);
however, this is only so because we output water content proﬁles
at the same time of the irrigation cycle (i.e. just before irrigation
starts, see Supplementary Material Figure S3). On the other hand,
there is a gradual increase in water content with depth as a resultof the plant water uptake in the top 0.5 m (rooting depth is from 0
to 0.5 m).
The simulated SAR and EC evolutions illustrate high SAR valuesof over 100 (Fig. 5A and B) and EC values exceeding 10 dS/m by theend of the simulation period (Fig. 6A and B). The rate at which SAR
moves through soil is determined by the migration of Na, Ca, and
Mg. Because Na, Ca, and Mg migration is retarded compared to thetracer (Fig. 7A), SAR migration (Fig. 5A) will be retarded too compared to the tracer. The precipitation of Ca in the top of the proﬁlefurther inﬂuences the value of SAR: a decrease in dissolved Ca as aresult of calcite production (Fig. 8) increases SAR. Note that pH andbicarbonate concentrations are internally calculated even thoughinputting such data is not required. The SAR ANZECC trigger valueof 20 (sandy-sandy loam soil, see Table 1) is exceeded relatively
quickly, i.e. down to a depth of 0.4 m after 150 days.

The electrical conductivity (Fig. 6) progresses much faster than
SAR because it is mainly determined by the non-reactive and fastermoving chloride; the EC reaches half its maximum concentration(6 dS/m) at the bottom of the proﬁle after 150 days (Fig. 6A), while
after 150 days the SAR value is at half its maximum (SAR = 60)
between 0.3–0.4 m depth (Fig. 5A). The EC ANZECC trigger valueof 0.65 dS/m (see Table 1) is exceeded very quickly after the startof irrigation, i.e. from the surface down to 0.7 m after 30 days.
From the six minerals considered in the Major Ion Chemistry(calcite, dolomite, gypsum, nesquehonite, hydromagnesite, sepeolite), only calcite is formed. The concentration of calcite gradually
increases with depth as more irrigation water enters the soil
(Fig. 8A). The effect of calcite precipitation is clearly reﬂected in
the concentration of dissolved calcium ions (Fig. 7A): close to thesurface, dissolved calcium is nearly absent as a result of precipitation of calcite. Using a calcite mineral density 2.72 g/cm3 , total
porosity (based on saturated water content) and soil bulk densityfrom Table S2 (Supplementary Material), the calculated mineralvolume of precipitated calcite after 365 days was shown to haveincreased to approximately 0.015 cm3 /cm3 at the soil surface. In
other words, the mineral phase occupies nearly 3% of total soil

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

61

Fig. 5. Simulated SAR for Scenario 2 (A), Scenario 3 (B), Scenario 4 (C) and Scenario 5 (D). ANZECC (2000) guide value for SAR is 5 for light clay and 20 for sand-sandy loamsoil. Note different scales in different ﬁgures.

porosity (0.515 cm3 /cm3 in soil layer 1). While this is a relatively
small fraction of the available pore space, this volume is expected
to increase if irrigation continues for many years. These long-term
effects of salt formation on soil hydrology and plant response needfurther corroboration.
There are clear differences in the way the major cations Ca, Mg,
Na and K migrate through the soil (Fig. 7). Calcium is mainly affectedby precipitation in the top soil; ahead of the precipitation front
calcium is only weakly retarded compared to chloride. Magnesium
displays accumulation at the soil surface due to a decrease in soil
water owing to ET; the similarity between magnesium and chloride
indicates the former is little retarded. Sodium displays the largest
retardation as is expected from the larger Gapon cation exchange
constant (see Table 3). The migration of potassium is somewhat
similar to Mg, although it exhibits nearly constant values in the top
0.3 m of the soil proﬁle.

3.3. Major ion chemistry effects on soil hydrology
Scenario 3 was run using the Major Ion Chemistry module allowing calculation of concentrations of major ions associated withirrigation of untreated produced water, while now also considering
reduction of the hydraulic conductivity due to solution composition (based on the McNeal model, Eq. (3)). During the simulations
and based on an irrigation rate of 6 cm/day once in ﬁve days, thehydraulic conductivity was reduced automatically due to the salinity effect up to the point that water could not get into the proﬁlewithout causing ponding. Therefore, we reduced the irrigation rate
to 4 cm/day once in ﬁve days to avoid ponding conditions. Under
these new soil hydrological conditions, the following observations
were made regarding the soil water balance and major ion chemistry. Changes in SAR and electrolyte concentration will also affectclay swelling, which in turn may affect the soil water retention
curve. Because the focus of this study was on reduction of the

62

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

Fig. 6. Simulated electrical conductivity for Scenarios 2 (A), 3 (B), 4 (C) and 5 (D). ANZECC (2000) EC trigger value of 0.65 dS/m. Note different scales in different ﬁgures.

hydraulic conductivity, effects on water retention curve were notconsidered.
The cumulative actual transpiration becomes severely reduced,
i.e. down to 188 cm for Scenario 3 compared to 240 cm for Scenario 2 (Fig. 2C and Table 4). The decrease in root water uptake
between 0.2 and 0.4 m soil depth (Fig. 9B) is not due to the salinity stress caused by a too high EC. This decrease occurs at the very
end of the simulation, when the SAR front reaches the second soil
layer, causing the hydraulic conductivity to decrease, which in turncauses water logging in this soil layer. The pressure head associated with this water-logged layer exceeds the −25 cm parameter
P1 (Fig. 9A), which caused the water stress response function ˛(h)to decrease below its optimal rate (i.e. less than one): this, in turn,caused a reduction in plant water uptake rate S(h) (a deviation fromits potential value (straight line at t = 0 days, with S(h) = 0.032 day−1
at z = 0 in Fig. 9B) and a reduction in actual transpiration.
Compared to the two previous scenarios, the SAR moves downwards about two times slower compared to Scenario 2, while themagnitude is more or less similar (Fig. 5B). The electrical conductivity displays a different distribution across the proﬁle: higher EC
values exist at 0.4 m depth while at greater depths EC values aremuch smaller (Fig. 6B).

3.4. Irrigation with amended produced water
Scenario 4 considers the use of amended irrigation water (basedon a 3:1 surface water:produced water ratio) with the chemical composition deﬁned in Table 2, while all other processes andparameters are identical to those of Scenario 3. The cumulative
potential and actual surface and bottom ﬂuxes shown in Table 4
reveal that the use of lower salinity irrigation water results in a signiﬁcantly reduced salinity stress: the actual transpiration is nearly
equal to its maximum value (281 cm, as in Scenario 1), 41 cm greaterthan in Scenario 2 (untreated produced water) and 93 cm greaterthan in Scenario 3 (untreated produced water with considerationof the hydraulic conductivity reduction model). As a result, thedrainage across the bottom of the soil proﬁle is 39 cm lower thanin Scenario 2 (for the same amount of irrigation water). Comparedto Scenario 2, this scenario results in a more efﬁcient use of theapplied irrigation water.
The simulated SAR proﬁles shown in Fig. 5C demonstrate that
the SAR ANZECC trigger value (SAR = 20) for sand-sandy loam isexceeded after 365 days, but only in a small section of the proﬁle.
The SAR ANZECC trigger value (SAR = 5) for light clay is exceeded in
the top 0.3 m after 150 days, whereas after 365 days it is exceededin the top 0.5 m.

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

63

Fig. 7. Simulated concentrations of dissolved major ions Ca (A), Na (B), Mg (C), K (D), and Cl (A, dashed line) for Scenario 2). ANZECC (2000) sodium trigger value of 20 meq/L(B). Dashed lines are for chloride.

The EC ANZECC trigger value of 0.65 dS/m is exceeded across theentire soil proﬁle after 90 days. The relatively high EC is due to the(still relatively high) EC of the irrigation water (1.57 dS/m). In termsof contributions of individual ions in the amended irrigation water
to the EC, the most important contributors to EC are sulfate (a 36%
contribution to total EC), sodium (a 26% contribution to total EC)
and chloride (a 23% contribution to total EC) (Table 2).
Optimisation of irrigation water quality can be done by reducing
concentrations of those ions that contribute most to the total EC.
In the above example, in decreasing order of importance, that issulfate, sodium, and chloride, where the contribution of each ion
to the total EC is calculated using the procedure given in Table 2.
3.5. Irrigation with treated produced water
In Scenario 5 irrigation water has the chemical composition ofwater treated by reverse osmosis, as deﬁned in Table 2, while allother processes and parameters are identical to those of Scenarios

3 and 4. The water balance is very similar to Scenario 4, with cumulative actual transpiration equal to 280 cm and drainage ﬂux equal
to −84 cm (Fig. 2D and Table 4). Time series of soil water content
and pressure head at the 10 cm depth display wetting-drying cycles
as a result of the applied irrigation regime and water redistribution
due to evapotranspiration and drainage (Fig. 10). Every irrigation
event adds solutes to the soil which partly migrate downwards
while another part accumulates in the top of the soil proﬁle. The
pressure head data indicates that the maximum pressure (about−15 cm) is slightly smaller than the pressure head below which root
water uptake is optimal (P1 = −25 cm) (Fig. 10). On the other hand,
the minimum pressure head (about −272 cm) exceeds the limitingpressure head below which root water uptake is no longer optimal(P2 = −200 cm). In other words, the water uptake by the grass vegetation considered here occurs at a maximum rate for about 3 outof 4 days, whereas on the 4th day it is suboptimal (see schematic
root water uptake function in Supplementary Material Figure S2).

64

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

Fig. 8. Calcite concentration for Scenario 2 (A). Calcite mineral volume after 365 days of irrigation (B).

Fig. 9. Pressure head (A) and root water uptake (B) proﬁles for Scenario 3. The pressure head exceeds the P1 = −25 cm threshold after 300 days (A).

A more optimal irrigation regime would involve starting irrigation
at day 4 to avoid the soil from becoming too dry.
The simulated SAR proﬁles now remain below the lowest SARANZECC trigger value (SAR = 5) for light clay (Fig. 5D). This result
is as expected, given the low SAR value (i.e. 1.6) of the irrigationwater and the absence of any initial salt build-up in the soil proﬁle(the initial SAR value is 1.25 throughout the entire soil proﬁle) orwithout other salt inputs into the soil proﬁle (e.g. from saline shallow groundwater or lateral inﬂow from water discharging at the
break-of-slope, Biggs et al., 2012).
The EC ANZECC trigger value of 0.65 dS/m is generally notexceeded across the entire soil proﬁle, except at the very top
few centimetres. Although the EC of the irrigation water is low(0.29 dS/m), EC in the top of the soil proﬁle increases to a maximum value of 0.73 dS/m due to evapotranspiration leading to dry

soil conditions, triggering salts to concentrate near the soil surface.
In terms of contributions of individual ions in the treated irrigationwater to the EC, the most important contributors to EC are chloride(a 49.6% contribution to total EC), sodium (a 21.5% contribution tototal EC) and magnesium (a 11.3% contribution to total EC) (Table 2).
Depth distributions of aqueous concentrations of major cations
Ca, Mg, Na and K in the soil proﬁle are shown in Fig. 11. Calcium,magnesium, and sodium display similar patterns; an accumulation
at the soil surface due to a decrease in soil water owing to ET and
gradual accumulation at deeper depths due to downward leaching. Sodium displays the largest retardation as is expected fromthe Gapon cation exchange constant (see Table 3). The migrationof potassium is somewhat different in the top 0.3 m of the soil proﬁle, without a distinctive accumulation as observed with the otherD. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

65

Fig. 10. Calculated soil water contents and pressure heads at the 10 cm depth (Scenario 5; pasture). Every ﬁve days a 24-h irrigation event takes place. P1 = pressure head
below which root water uptake is optimal and P2 = limiting pressure head below which root water uptake is no longer optimal.

ions. The sodium ANZECC trigger value of 20 meq/L (460 mg/L, see
Table 1) is not exceeded anywhere in the soil proﬁle.
Fig. 12 shows concentrations of sorbed/exchangeable cations
on the cation exchange sites. While the concentration front for
exchangeable Na moves relatively quickly down the soil proﬁle,the front reaching half the total soil depth after 365 days, otherconcentration fronts (i.e., for Ca and Mg) move much slower. Thisis mainly due to large amounts of Na being added with irrigationwater and relatively small buffering capacity with respect to Na(note that the initial concentration for exchangeable Na was only
about 2.6 meq/kg). On the other hand, the initial concentrations of
exchangeable Ca and Mg were about 175 and 89 meq/kg, resulting in a signiﬁcantly larger buffering capacity for these two cations
compared to Na.
From the six minerals considered in the Major Ion Chemistry (calcite, dolomite, gypsum, nesquehonite, hydromagnesite,
sepeolite), only calcite is formed in Scenario 5. However, the concentration of calcite is very small as it reaches only 0.5 meq/kg(results not shown), compared to several hundred meq/kg for Scenario 2.
3.6. Summary of water and chemical mass balance calculations
For each scenario, water balance components were calculated
based on a 365-day simulation period. Cumulative inﬁltration,
actual transpiration (root water uptake), evaporation and net
drainage (outﬂow–inﬂow) values are shown for each of ﬁve scenarios evaluated (Table 4). The cumulative actual transpiration andevaporation are each also expressed as a percentage of their respective cumulative potential values (292 and 73 cm, respectively).
For all scenarios except Scenario 3 cumulative inﬁltration equals
the imposed value of 438 cm based on a 6 cm/day irrigation rate
once every 5 days. Scenario 3 had a reduced irrigation rate of
4 cm/day once every 5 days resulting in a lower cumulative inﬁltration of 292 cm. The lower irrigation rate was required to avoid
ponding conditions as a result of the reduction in hydraulic conductivity when untreated produced water was applied.
The cumulative actual transpiration is a useful metric to evaluate combined effects of soil chemical and soil physical processeson the soil water balance and on plant water stress. In all scenarios

the cumulative actual transpiration is less than its potential value;
for Scenarios 1, 4, and 5 transpiration reaches 281 cm or 96% of
the potential value. The observation that the plant experiences avery small water stress is most likely due to the sub-optimal irrigation regime which generates soil pressure heads slightly outside
the optimal range, especially near the more negative soil pressure
heads where the limiting pressure head P2 becomes exceeded (seeFig. 10). Larger differences between actual and potential transpiration exist for Scenario 2 (82% of potential value) and 3 (64% of
potential value). Scenario 2 applied the combined water and solutestress model - in combination with the major ion chemistry module - to control root water uptake: this resulted in an additional
stress on plant water uptake, mainly caused by a salinity stress(irrigation water EC = 5700 ␮S/cm). Smaller actual transpiration for
Scenario 3 (188 cm) is due in part to (i) the increased water stressas a result of a drier soil (4 cm/day rather than 6 cm/day irrigation rate, every 5 days), (ii) sup-optimal pressure head conditions
owing to a decrease in the hydraulic conductivity, generating nearsaturated conditions within the root zone, and (iii) to a salinity
stress (irrigation water EC = 5700 ␮S/cm).
Cumulative actual evaporation reached its maximum value of
73 cm (based on a 0.2 cm/day potential evaporation rate) in all but
one scenario (i.e., Scenario 3 yielded 69 cm or 95% of the potential value owing to the reduced inﬁltration, generating drier soil
moisture conditions and hence reduced evaporation).
Cumulative net drainage (the difference between cumulative
outﬂow and inﬂow) is identical (and negative, deﬁned here as outﬂow through the bottom of the soil proﬁle) for Scenarios 1, 4, and 5.
Higher drainage for Scenario 2 (−123 cm) is due to the lower transpiration, leaving more water in the soil proﬁle that is available forgravity drainage. Much lower drainage for Scenario 3 (−29 cm) is
the result of the smaller irrigation rate. Note that the decrease inirrigation by ∼33% results in a decreased in drainage of ∼65%. Calculated net drainage values from Table 4 are considerably higher
than calculated values for furrow irrigation on cracking clay soils(e.g., 10–20 cm/year) as reported by (Biggs et al., 2012). Our higher
drainage values are the result of the higher irrigation rate applied
constantly for a full year.
A comparison of Scenarios 1 and 2 is interesting from anotherperspective: it illustrates the effect of including or excluding certain

66

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

Fig. 11. Simulated concentrations of dissolved major ions Ca (A), Na (B), Mg (C), K (D), and Cl (A), Scenario 5. Dashed lines are for chloride.

coupled soil processes on soil water balance components. Scenario
1 has a standard advection–dispersion solute transport model without major ion chemistry; the irrigation water is solute-free andonly contains a tracer. Scenario 2, on the other hand, accounts not
only for major ion chemistry but also calculates salinity effects onroot water uptake. Scenario 2 is an example where a link is madebetween soil geochemistry, soil hydrology, and plant water uptake.
Accounting for such coupled processes resulted in (i) a decrease of
the actual transpiration to 82% of its potential value, and (ii) an
increase in net drainage by 46% compared to Scenario 1. This isa clear demonstration that soil water balance modelling involving major ion chemistry should account for the above interacting
processes. Note that the importance of such interacting processes
decreases as the applied irrigation water has a less extreme chemical composition (such as the amended and treated produced water).
Mass ﬂuxes of major cations and anions were calculated forupper and lower model boundary (Table 5). This is important for
estimating the potential effects of salt leaching from upper soil
layers to deeper soil layers as leaching may negatively affect theconcentration of adsorbed sodium and the soil hydraulic conduc-tivity of deeper soil layers. To better protect groundwater resourcesfrom pollution with salts and nutrients, it is further importantto estimate the potential mass ﬂux of salts leaching to shallow
groundwater. For all scenarios, the mass ﬂux for calcium out of
the soil is larger than the mass of calcium added via irrigation.
This results from calcium being replaced on the adsorption complex by mono-charged sodium cations whose added mass via
the irrigation water is two orders of magnitude larger than thatof calcium. As the sodium concentration in the irrigation water
decreases when untreated produced water is replaced by surfacewater amended produced water (Scenario 4) or treated produced
water (Scenario 5), the calcium mass ﬂux through the bottom of
the soil proﬁle decreases due to decreased competition for adsorption sites between calcium and sodium. For magnesium, the mass
ﬂux through the bottom of the soil proﬁle is larger than that added
via irrigation only when untreated produced water is used (Scenarios 2 and 3). Similarly as for Ca, Mg is displaced from the solid
phase by the large concentrations of Na. For all other scenarios,
magnesium accumulates in the soil (input ﬂux is larger than output ﬂux via drainage). As a result of the displacement of Ca and Mg

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

67

Fig. 12. Surface species concentrations for Ca (A), Na (B), Mg (C) and K (D) versus soil depth during irrigation cycles. The total adsorbed concentration for Ca + Na + Mg + K isconstant (275 meq/kg for soil layer 1, 278 meq/kg for soil layer 2, and 274 meq/kg for soil layer 3) (Scenario 5).

Table 5Chemical mass ﬂuxes (g/m2 ) through upper and lower model boundary for a 365 day simulation period. Negative ﬂuxes for the upper model boundary are inputs to the soil;
negative values for the lower model boundary refer to ﬂuxes leaving the soil. For each element and each scenario, the ﬁrst row is a ﬂux through the upper boundary (input);the second row is a ﬂux through the lower boundary (output).
Scenario

Ca

Mg

Na

K

HCO3

SO4

Cl

2: Pasture –UPW–MIC–WSS−4.56E+01−6.76E+02−3.93E+01−7.91E+01−6.18E+03−6.90E+02−2.41E+01−7.35E+01−5.64E+03−1.44E+01−9.99E+01−7.92E+01−2.61E+03−2.07E+03

3: Pasture –UPW–MIC–WSS–KRM−3.02E+01−4.42E+02−2.59E+01−4.99E+01−4.07E+03−1.69E+02−1.59E+01−3.19E+01−3.72E+03−2.92E+00−6.58E+01−3.80E+01−1.72E+03−9.93E+02

4: Pasture –SAPW–MIC–WSS–KRM−4.52E+01−2.06E+02−2.86E+01−2.07E+01−8.09E+02−5.52E+01−1.53E+02−2.55E+01−6.04E+02−1.01E+01−1.48E+03−9.13E+02−7.20E+02−4.43E+02

5: Pasture –TPW–MIC–WSS–KRM−4.54E+01−5.70E+01−3.33E+01−6.27E+00−1.26E+02−2.80E+01−1.90E+00−1.63E+01−8.24E+01−1.17E+01−4.39E+00−2.74E+00−2.95E+02−1.75E+02

UPW: untreated produced water; MIC: major ion chemistry; WSS: water and solute stress model; KRM: hydraulic conductivity reduction model; SAPW: surface water
amended produced water; TPW: treated produced water.

68

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69

by Na, sodium accumulates on the exchanger resulting in a tenfold smaller sodium ﬂux at the bottom of the soil proﬁle compared
to its input ﬂux (Scenario 2). The behaviour of potassium is similar to that of magnesium: as a smaller mass of potassium goes
into the soil under Scenarios 4 and 5, the mass ﬂux of potassium
leaving the soil proﬁle decreases. When comparing the effect ofusing different water types on salt leaching, the use of treated produced water clearly generates the overall smallest mass ﬂuxes (forall ions except HCO3 ) from soil to groundwater. This is not just theresult of a smaller input of solutes via the treated produced water,
but also because decreased sodium concentrations displace smaller
amounts of major ions previously present on the soil’s exchange
complex.

4. Conclusions
An in-depth analysis has been undertaken to provide information for potential optimisation of irrigation water quality that will
yield minimal long-term effects of dissolved ions on soil and plant
health. Three different irrigation water qualities were evaluated,
where the water source was considered to be produced coal seamgas water: untreated or treated produced water or produced water
amended with surface water. The analysis also involved the use of
two different conceptualisation of water quality evolution in thesoil. The ﬁrst conceptualisation served as a baseline and did only
consider the coupled variably-saturated water ﬂow and standard
advection–dispersion transport models implemented in HYDRUS1D. The second conceptualisation considered the HYDRUS-1D codeand its major ion chemistry module UnsatChem to simulate thesimultaneous movement of multiple major ions (Na, K, Ca, Mg, Cl,
SO4 , alkalinity) present in irrigation water in a vegetated soil proﬁle(pasture).
Simulations with different irrigation water qualities provided
detailed results regarding chemical indicators of soil and plant
health, i.e. SAR, EC and sodium concentrations. By comparing such
indicators in the soil proﬁle with water quality guideline values(here the Australian and New Zealand ANZECC values), an assessment was made of the suitability of the applied produced water
for long-term irrigation. The assessment evaluated the build-up ofsalts in the soil proﬁle and its effect on the plant salinity stress,
a condition which reduces the capacity of plants to uptake waterand results in yield reduction. Finally, the simulations also allowed
to test under which conditions soil hydraulic properties, in particular, the hydraulic conductivity, are negatively impacted by salt
concentrations owing to clay dispersion. Changes in the soil water
retention characteristics were not modelled.
One of the three irrigation scenarios considered the use of
untreated produced water (EC = 5705 ␮S/cm and SAR = 77). Because
irrigation with produced water is strictly regulated, the purpose ofthis scenario was to provide useful insights into the type of coupled chemical-hydrological processes that might occur, and what
the potential impacts would be on soil and plant health shouldsuch water end up in the soil. For example, by using saline irrigation water, we demonstrated that accumulation of monovalent ionscauses clay swelling, resulting in a decrease in hydraulic conductivity of a particular soil layer. This, in turn, caused water stagnation
above this layer (the pressure head exceeds the threshold for optimum water uptake by plant roots), leading to an additional plantwater stress and reduced transpiration (188 cm versus 240 cm
without the water stress, or more than 20% less). Due to lower
transpiration, more water resides in the soil proﬁle after irrigation, causing drainage to increase by 46% (−123 cm compared to−84 cm).
In case when produced water was amended with surface
water (at a ratio of 3:1 of surface water:produced water, with

EC = 1573 ␮S/cm and SAR = 11), the calculated soil SAR values weremuch lower and generally acceptable for sandy to sandy-loam soil
(i.e. less than 20). Permissible SAR values for light clay (i.e. lessthan 5) were still exceeded in most of the soil proﬁle. Results also
demonstrated the formation of calcite in the top of the soil proﬁle owing to solute concentration effects under dry soil moisture
conditions. This case did not result in an additional stress on plantwater uptake as a result of the acceptable salinity levels.
A ﬁnal case used treated produced water with a relatively low
concentration of major ions (EC = 291 ␮S/cm, SAR = 1.6). The permissible SAR values for light clay were not exceeded, and all otherions remained at tolerable levels. Under these conditions, actual
transpiration was not reduced by an additional salinity stress,
meaning that near optimal soil hydrological and chemical conditions existed for plant growth.
Our analysis further illustrates that accounting for coupled geochemical, hydrological and plant water uptake processes resultedin (i) a decrease of the actual transpiration to 82% of its potential
value, and (ii) an increase in net drainage by 46% when compared toa water balance calculation where such interacting processes weredisregarded. This is a clear demonstration that soil water balance
modelling involving major ion chemistry should account for theabove interacting processes.
The major ion chemistry module UnsatChem implemented in
HYDRUS-1D is a very powerful tool to assess in a very detailed way
potential effects of different irrigation water qualities on soil andplant health. This offers a cost-effective way to optimise the qualityof irrigation water derived from coal seam gas water or other saline
waters and thus to ensure soils are managed in a sustainable manner. To further underpin such assessments, ﬁeld trials providing
detailed soil chemical and physical data need to be considered totest the validity of the simulations and to improve process parameterisation. One of the processes that needs improvement andtesting is the change in hydraulic properties as a result of mineral precipitation/dissolution causing partial or full pore blockage.
Based on a simulation period of a single year, calcite precipitation
was found to occupy up to 3% of the soil’s pore space in the top of
the soil proﬁle. The long-term evolution of salinity and its effect onsoil hydraulic properties also needs further corroboration.
Acknowledgements
The authors would like to acknowledge Sarah Bennett (formerly
University of NSW) for providing detailed soil hydraulic property
data on the Cox Creak Catchment (NSW). We equally acknowledge
the Land and Water Business Unit for funding through the strategic appropriation research project “Next generation methods andcapability for multi-scale cumulative impact and management”.
Appendix A. Supplementary data
Supplementary data associated with this article can be found,in the online version, at http://dx.doi.org/10.1016/j.agwat.2017.04.
011.
References
Appelo, T., Postma, D., 2004. Geochemistry, Groundwater, and Pollution, second ed.
ANZECC, 2000. Australian and New Zealand Guidelines for Fresh and Marine Water
Quality. Volume 1: The Guidelines (Chapters 1–7).
Beletse, Y.G., Annandale, J.G., Steyn, J.M., Hll, I., Aken, M.E., 2008. Can crops be
irrigated with sodium bicarbonate rich CBM deep aquifer water? Theoreticaland ﬁeld evaluation. Ecol. Eng. 33, 26–36.
Bennett, S., (Ph.D. thesis) 2012. Key Spatiotemporal Factors Determining
Uncertainty of Deep Drainage in a Semi-arid Area. The University of New South
Wales.
Biggs, A., Witheyman, S.L., Williams, K.M., Cupples, N., de Voil, C.A., Power, R.E.,
Stone, B.J., August 2012. Assessing the Salinity Impacts of Coal Seam Gas Water

D. Mallants et al. / Agricultural Water Management 189 (2017) 52–69
on Landscapes and Surface Streams. Final Report of Activity 3 of the Healthy
HeadWaters Coal Seam Gas Water Feasibility Study. Department of Natural
Resources and Mines, Toowoomba.
BOM, 2014. Bureau of Meteorology, Recent Evapotranspiration, http://www.bom.
gov.au/watl/eto/.
Bright, D.A., Addison, J., 2002a. Derivation of Matrix Standards for Salt Under the
British Columbia Contaminated Sites Regulation. Royal Roads University.
Bright, D.A., Addison, J., 2002b. Derivation of Matrix Standards for Salt Under the
British Columbia Contaminated Sites Regulation. Addendum A: Technical
Options Analysis for a Soil Chloride Standard for Drinking Water Protection(Aesthetic). Applied Research Division, Royal Roads University, Victoria, BC.
Burkhardt, A., Gawde, A., Cantrell, C.L., Baxter, H.L., Joyce, B.L., Stewart Jr., C.L.,
Zheljazkov, V.D., 2015. Effects of produced water on soil characteristics, plantbiomass, and secondary metabolites. J. Environ. Qual. 44, 1938–1947, http://dx.doi.org/10.2134/jeq2015.06.0299.
Crosbie, R., Morrow, D., Cresswell, R., Leaney, F., Lamontagne, S., Lefournour, M.,
2012. New Insights to the Chemical and Isotopic Composition of Rainfall
Across Australia. Water for a Healthy Country Flagship Report Series ISSN:
1835-095X. CSIRO.
DERM, 2010. Guideline: Approval of Coal Seam Gas Water for Beneﬁcial Use.
Department of Environment and Resource Management, Queensland,
Australia.
DEHP, 2014. General Beneﬁcial Use Approval—Irrigation of Associated Water(Including Coal Seam Gas Water). Department of Environment and Heritage
Protection, Queensland, Australia.
DNRM, 2012. Forecasting Coal Seam Gas Water Production in Queensland’s Suratand Southern Bowen Basins. Report Prepared by Klohn Crippen Berger:
Healthy HeadWaters Coal Seam Gas Water Feasibility Study. Department ofEnvironment and Resource Management, Queensland, Australia.
DNRM, 2013. Assessment of Alternative Use Options for Coal Seam Gas Water
Proposed for Central Condamine Alluvium Recharge Schemes. Healthy
HeadWaters Coal Seam Gas Water Feasibility Study. Department ofEnvironment and Resource Management, Queensland, Australia.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field Water Use andCrop Yield. John Wiley & Sons, New York, NY.
Horizon, 2006. Irrigated Pasture Systems Comparative Project. Stage 3 Final Report.
Jacques, D., Smith, C., Šimůnek, J., Smiles, D., 2012. Inverse optimization of
hydraulic, solute transport, and cation exchange parameters using HP1 and
UCODE to simulate cation exchange. J. Contam. Hydrol. 142–143, 109–125.
Jacques, D., Šimůnek, J., Mallants, D., van Genuchten, M.Th., 2013a. The HPx
reactive transport models: summary of recent developments and applications.
In: 4th International Conference “Hydrus Software Applications to Subsurface
Flow and Contaminant Transport Problems”. Czech University, pp. 7–16.
Jacques, D., Perko, J., Seetharam, S., Mallants, D., Govaerts, D., 2013b. Modelling
long-term evolution of cementitious materials used in waste disposal. In:
IAEA-TECDOC-1701 Behaviours of Cementitious Materials in Long Term
Storage and Disposal., pp. 26.
Johnston, C.R., Vance, G.F., Ganjegunte, G.K., 2008. Irrigation with coalbed natural
gas co-produced water. Agric. Water Manag. 95, 1243–1252.
Lenhard, R.J., Brooks, R.H., 1986. Effects of clay–solution interactions on water
retention. Am. Soc. Civil Eng. J. Irrig. Drain. Eng. 112, 28–38.
Leskiw, L.A., Sedor, R.B., Welsh, C.M., Zeleke, T.B., 2012. Soil and vegetation
recovery after a well blowout and salt water release in northeastern British
Columbia. Can. J. Soil Sci. 92, 179–190.
Maas, E.V., 1990. Crop salt tolerance. In: Tanji, K.K. (Ed.), Agricultural Salinity
Assessment and Management. ASCE Manuals and Reports on Engineering
Practice, No. 71. , New York.
Mace, J.E., Amrhein, C., 2001. Leaching and reclamation of a soil irrigated withmoderate SAR waters. Soil Sci. Soc. Am. J. 65, 199–204.
McNeal, B.L., 1968. Prediction of the effect of mixed-salt solutions on soil hydraulicconductivity. Soil Sci. Soc. Am. Proc. 32, 190–193.

69

McNeal, B.L., 1974. Soil salts and their effects on water movement. In: van
Schilfgaarde, J. (Ed.), Drainage for Agriculture, Agronomy No. 17. Am. Soc. Agr.,
Madison, WI.
Minasny, B., McBratney, A., 2002. The neuro-m method for ﬁtting neural network
parametric pedotransfer functions. Soil Sci. Soc. Am. J. 66, 352–361.
Moore, T.A., 2012. Coalbed methane: a review. Int. J. Coal Geol. 101, 36–81.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity ofunsaturated porous media. Water Resourc. Res. 12, 513–522.
Nghiem, L.D., Ren, T., Aziz, N., Porter, I., Regmi, G., 2011. Treatment of coal seam gas
co-produced water for beneﬁcial use in Australia: a review of best practices.
Desalin. Water Treat. 32, 316–323.
Parsons, W., 2010. Chapter 12 in Australia Paciﬁc LNG Project Environmental
Impact Statement. Adaptive Associated Water Management, vol. 2. Report toOrigin and Conoco Phillips (Asia Paciﬁc LNG Joint Venture), Brisbane.
QWC, July 2012. Underground Water Impact Report for the Surat Cumulative
Management Area. Queensland Water Commission.
Ringrose-Voase, A.J., 2004. Water Balance and Deep Drainage Under Irrigated
Cotton. ‘WaterPak’. In: Dugdale, H., Harris, G., Neilsen, J., Richards, D., Roth, G.,
Williams, D. (Eds.). Cotton Research and Development Corporation, Narrabri,
New South Wales, pp. 17–28.
RPS, 2011. Onshore co-produced water: extent and management. In: LTD, R.A.P.
(Ed.), Waterlines Report Series No. 54. National Water Commission, Canberra.
Santos, 2009. Sustainability Report 2009: The Energy for Sustainability, Australia.
Santos, 2010. Upstream-Fairview Project Area Environmental Management Plan;
Appendix B Fairview Project Area Coal Seam Gas Water Management Plan.
Silburn, D.M., Montgomery, J., 2004. Deep drainage under irrigated cotton inAustralia: a review. In: WATER-PAK a Guide for Irrigation Management in
Cotton. Cotton Research and Development Corporation and Australian Cotton
Cooperative Research Centre, Narrabri, pp. 29–40.
Šimůnek, J., Suarez, D.L., 1994. Two-dimensional transport model for variably
saturated porous media with major ion chemistry. Water Resour. Res. 30,
1115–1133.
Šimůnek, J., Suarez, D.L., 1997. Sodic soil reclamation using multicomponent
transport modelling. ASCE J. Irrig. Drain. Eng. 123, 367–376.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2008. Development and applicationsof the HYDRUS and STANMOD software packages and related codes. Vadose
Zone J. (Special Issue “Vadose Zone Modeling”) 7, 587–600, http://dx.doi.org/
10.2136/VZJ2007.0077.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.Th., 2013. The
HYDRUS-1D Software Package for Simulating the Movement of Water, Heat,and Multiple Solutes in Variably Saturated Media, Version 4.17, HYDRUS
Software Series 3. Department of Environmental Sciences, University ofCalifornia Riverside, Riverside, CA, USA, pp. 343.
Šimůnek, J., Jacques, D., van Genuchten, M.T., Mallants, D., 2006. Multicomponent
geochemical transport modelling using Hydrus-1D and HP1. J. Am. Water
Resour. Assoc. 42, 1537–1547.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J. 15(7), 25, http://dx.doi.org/10.2136/vzj2016.04.0033.
Stearns, M., Tindall, J., Cronin, G., Friedel, M., Bergquist, E., 2005. Effects of coal-bed
methane discharge waters on the vegetation and soil ecosystem in PowderRiver Basin, Wyoming. Water Air Soil Pollut. 168, 33–57.
Vance, G.F., King, L.A., Ganjegunte, G.K., 2008. Soil and plant responses from land
application of saline-sodic waters: implications of management. J. Environ.
Qual. 37, S139–S148.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
Wesseling, J.G., Elbers, J.A., Kabat, P., van den Broek, B.J., 1991. SWATRE:
Instructions for Input, Internal Note. Winand Staring Centre, Wageningen, The
Netherlands.
White, N., Zelazny, L.W., 1986. Charge properties in soil colloids. In: Sparks, D.L.
(Ed.), Soil Physical Chemistry. CRC Press, Boca Raton, FL.