Journal of Hydrology 596 (2021) 126053

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Adaptation and validation of the ParSWMS numerical code for simulationof water flow and solute transport in soilless greenhouse substrates
Mohammad R. Gohardoust a, Jirka Šimůnek b, Horst Hardelauf c, Markus Tuller a, *a

Department of Environmental Science, The University of Arizona, Tucson, AZ, USA
Department of Environmental Sciences, University of California, Riverside, CA, USAc
Institute of Bio-and Geosciences, IBG-3 Agrosphere, Jülich, Germanyb

A R T I C L E I N F O

A B S T R A C T

This manuscript was handled by Corrado Cor­
radini, Editor-in-Chief, with the assistance ofPhilip Brunner, Associate Editor

Numerical simulation of three-dimensional water flow and solute transport in containerized variably saturated
soilless substrates with complex hydraulic properties and boundary conditions necessitates high-resolution dis­cretization of the spatial and temporal domains, which commonly leads to several million nodes requiring nu­
merical evaluation. Even today’s computing prowess of workstations is not adequate to tackle such problems
within a reasonable timeframe, especially when numerous realizations are required to optimize the geometry,
substrate properties, and irrigation and fertigation management of soilless plant growth modules. Hence, theparallelization of the numerical code and utilization of high performance computing (HPC) are essential. Here,
we adapted and applied the ParSWMS parallelized code that is amenable to solving the 3D Richards equation forwater flow and the convection-dispersion equation for solute transport subject to linear solute adsorption. The
code was modified to allow for nonlinear equilibrium solute adsorption with new boundary conditions andapplied to simulate water flow and nitrogen and phosphorus transport in containerized soilless substrates. Multisolute transport simulations with the modified Linux ParSWMS code were first performed on a workstation andreferenced to the Windows-based HYDRUS (2D/3D) numerical code. After confirming the agreement between
the modified ParSWMS code and HYDRUS (2D/3D), various preconditioners and iterative solvers were evaluatedto find the computationally most efficient combinations. The performance of the modified ParSWMS code and its
stability were compared to HYDRUS (2D/3D) simulations for three soilless substrates consisting of horticultural
perlite, volcanic tuff, and a volcanic tuff/coconut coir mixture. Considering the solute mass balance error as astability measure, ParSWMS outperformed HYDRUS (2D/3D). Moreover, simulations with the modified
ParSWMS code were about 22% faster than simulations with HYDRUS (2D/3D) on the workstation. Tests of themodified ParSWMS on two HPC clusters with 28 and 94 cores revealed a potential computational speedup of
94% relative to the HYDRUS (2D/3D) simulations performed on the workstation.

Keywords:
Soilless culture
Water flow and solute transport
Numerical simulations
High performance computing
HYDRUS (2D/3D)
ParSWMS

1. Introduction
Because of the urgent need to secure and sustain the food and watersupply for an ever-growing human population, especially in underde­veloped arid and semiarid regions of the world, and an increasing de­mand for out-of-season fruits, vegetables, and ornamentals in theindustrial world, there is a momentous incentive to shift from soil-based
crop production to more resource-efficient containerized soilless pro­duction systems (Raviv et al., 2019). Soilless substrates exhibit majoradvantages over soils. Besides the alleviated risk for spreading soilborne
pathogens, their physical and hydraulic properties (i.e., plant water

availability after irrigation and the aeration characteristics) are superior
to those of arable soils (Savvas, 2003; Savvas and Gruda, 2018).
Furthermore, the nutrient availability to plant roots can be better
controlled in soilless substrates, which commonly leads to healthier
plants and higher yields than in soil-based production (Raviv et al.,2019). It is even possible to tailor (i.e., engineer) substrates that exhibit
specific physical and hydraulic properties beneficial for the growth ofdifferent crops. Though the same physical principles apply to bothsoilless substrates and soils, their physical and hydraulic properties arevastly different.
In addition, there are fundamental differences with regard to* Corresponding author.
E-mail address: mtuller@arizona.edu (M. Tuller).
https://doi.org/10.1016/j.jhydrol.2021.126053
Received 14 October 2020; Received in revised form 31 December 2020; Accepted 28 January 2021
Available online 10 February 2021
0022-1694/© 2021 Elsevier B.V. All rights reserved.

M.R. Gohardoust et al.

Journal of Hydrology 596 (2021) 126053

dynamic water, air, and nutrient distribution processes and root growthand development between spatially confined growth containers andunconfined field soils. While an impermeable container bottom (withdrainage holes) restricts water flow and nutrient transport in growthcontainers, water drains and redistributes to much deeper layers inagricultural soils unless natural impediments exist. This leads to vastlydifferent infiltration and redistribution dynamics, requiring more
intensive management of soilless systems (Gohardoust et al., 2020).
Taking into consideration that currently most of the soilless culture
plant growth experiments are based on trial and error, threedimensional (3D) numerical flow and transport simulations can pro­vide a scientifically sound basis for a priori elimination of substrates
with unfavorable growth potential, as well as guidance for growthmodule design and irrigation and fertigation management, thereby
preventing costly mistrials.
Today, numerical computer codes that solve the governing 3D waterflow and solute transport equations are indispensable for advancing theunderstanding of complex porous media processes. However, intricate
hydraulic properties and boundary conditions demand high-resolution
discretization of the spatial and temporal domains, often leading toseveral million nodes requiring numerical evaluation (Šimůnek et al.,2016). Even the current computing prowess of workstations is notadequate to tackle such simulations within a reasonable timeframe,
especially when numerous realizations are desired to optimize the ge­
ometry, soilless substrate properties, and irrigation and fertigation
management of plant growth modules. Hence, parallelization of thenumerical codes and utilization of high-performance computing (HPC)
clusters are vital.
Hardelauf et al. (2007) were among the first to develop a parallelized
code for 3D simulation of water flow and solute transport that theytermed ParSWMS. The code is based on the SWMS-3D model of Šimůnek
et al. (1995), which numerically solves the Richards equation for satu­
rated/unsaturated water flow and the convection–dispersion equation(CDE) for solute transport using a Galerkin-type linear finite element
scheme. The resulting nonlinear partial differential equation (PDE) forwater flow and the linear PDE for the CDE (i.e., assuming linear
adsorption isotherms) are solved in SWMS-3D with preconditioned
conjugate gradient and Orthomin (Vinsome, 1976) methods, respec­tively. In ParSWMS, Hardelauf et al. (2007) employed ParMETIS, an
open source MPI-based library, to distribute the simulation to multiple
processing nodes, and the Portable, Extensive Toolkit for Scientific
Computation (PETSc) library (https://www.mcs.anl.gov/petsc/) forpreconditioning and solving the resulting system of linear equations.
Driven by the need to execute a vast number of 3D flow and transport
simulations to aid with the design and management of containerized
plant growth modules, we adopted ParSWMS and modified the code to:(1) permit the simultaneous application of time-variable flux and freedrainage boundary conditions (BCs) and allow for the time-variable flux
BC to be treated as an atmospheric BC; and (2) enable multi-solutetransport simulations with nonlinear equilibrium adsorption. These
modifications are essential for the simulation of surface drip irrigationand to account for the nonlinear adsorption behavior of phosphorus andammonium.
Using HYDRUS (2D/3D) version 3.02.0560 (Šimůnek et al., 2018) asa reference, we first meticulously verified the modified ParSWMS code
for a simple, yet realistic, soilless growth module setup. Thereafter, we
focused on computational efficiency and evaluated a vast number of
preconditioner/solver combinations for both water flow and solutetransport based on a real growth module used in production scale
greenhouse experiments. Subsequently, we benchmarked the soilless
substrate flow and transport simulation results against HYDRUS (2D/3D) on a workstation with both Windows and Linux OS installed. As a
final step, we executed the modified ParSWMS code on two HPC clusters
to demonstrate the significant computational speedup relative to the
HYDRUS (2D/3D) simulations performed on the workstation. Please
note that for brevity we refer to HYDRUS instead of HYDRUS (2D/3D) in

the following.
2. Modifications to ParSWMS
2.1. Boundary conditions
In the ParSWMS water flow simulation each node is assigned an
integer code ranging from − 6 to 6. While positive values are allocated tonodes with prescribed pressure heads, zero and negative values repre­
sent nodes with prescribed volumetric water fluxes. In the original
ParSWMS code (Hardelauf et al., 2007) both time-variable flux and freedrainage boundary nodes are assigned –3. Hence, these BCs cannot be
used simultaneously. To resolve this issue, the free drainage boundary
condition code was set to − 6, similar as in HYDRUS, which was notallocated in the original ParSWMS. The discharge rate at a free drainage
node n is determined as Q(n) = − width(n) × K(h), where width(n) is thesurface area associated with node n, and K(h) is the hydraulic conduc­tivity as a function of the pressure head h.
To account for stage-I (i.e., potential evaporation rate controlled byatmospheric demand) and stage-II (i.e., falling evaporation rate limitedby the ability of the substrate to transmit water to the surface) evapo­
ration, it is necessary to treat the variable flux BC imposed for simulating
surface drip irrigation as an atmospheric BC. The atmospheric BC that is
implemented in both HYDRUS and ParSWMS is expressed as (Neuman
et al., 1974):) ⃒⎧ ⃒⃒ ( ∂h⃒⃒K K A + K A ni ⃒ ≤ E
⎪
⎪
ij
iz
⃒
⎪⃒
⎪
∂
xj
⎨
(1)
⎪
and
⎪
⎪
⎪
⎩
hCritA ≤ h ≤ hCritS
where K is the hydraulic conductivity, KAij are the components of theanisotropy tensor K A , xi (i = 1,2,3) are the spatial coordinates, ni are thecomponents of the outward unit vector normal to the atmosphericboundary, E is the maximum potential evaporation rate, and hCritS and
hCritA are two limiting values for the surface pressure head h. The
Einstein summation convention is applied to Eq. (1), which implies
summation over all possible values of the index in algebraic terms withrepeated indices. While the hCritS value specifies the maximum allowed
pressure head at the soil surface (usually 0.0), the hCritA value repre­
sents the minimum allowed surface pressure head defined based on
equilibrium conditions between soil water and atmospheric vapor. The
modified ParSWMS code now automatically switches between Dirichlet
and Neumann boundary conditions for nodes with a variable flux BC if
one of these limiting points is reached.
2.2. Multi-solute transport with nonlinear equilibrium adsorption
To allow realistic simulations of phosphorus and ammonium trans­
port processes in containerized soilless substrates, an option for
nonlinear equilibrium solute adsorption was added to the modified
ParSWMS code. The underlying CDE in ParSWMS is given as:()∂θc ∂ρs∂∂c∂qi c+ μw θc + μs ρs + γw θ + γ s ρ − Sc
(2)
θDij
−
+
=
∂xi
∂t
∂t
∂xj
∂xi
where θ is the volumetric substrate water content [L3 L− 3], c is the so­lution concentration [M L− 3], s is the concentration adsorbed to the solid
phase [M M− 1], ρ is the bulk density of the substrate [M L− 3], Dij is thedispersion coefficient tensor [L2 T− 1], qi is the i-th component of thevolumetric flux [LT− 1], μw and μs are first-order rate constants for sol­utes in the liquid and solid phases [T− 1], respectively, γ w and γs are zeroorder rate constants for the liquid [M L− 3 T− 1] and solid [T− 1] phases,
respectively, S is the sink term in the water flow equation [T− 1], and c isthe concentration of the sink term [M L− 3] with potential i andj values of2

M.R. Gohardoust et al.

Journal of Hydrology 596 (2021) 126053

1, 2, 3. The equilibrium adsorption isotherm in its general form can be
written as:
s=

ks cβ
1 + ηcβnected to a pressure manifold with a high-resolution pressure/vacuum
regulator and gauge. Initially, saturated samples were sequentiallydesaturated via application of increasing pressures after each equili­
bration phase. The saturated hydraulic conductivities (Ksat ) were
measured with an automated constant head permeameter. Details about
the SWC and Ksat measurements are provided in Gohardoust et al.
(2020).
After thoroughly testing several SWC models comprising the van
Genuchten (1980) model, a modified version of the van Genuchten
model (Vogel and Cislerova, 1988; Vogel et al., 2000), and the Brooks
and Corey (1966) model, we found that the latter was best suited for
mitigation of numerical instability issues caused by the extreme
nonlinearity of the hydraulic conductivity function near saturation for
all considered soilless substrates (Fig. 1). The Brooks and Corey (BC)SWC model is given as:⎧( )λ⎪⎨ θ(h) = θ + (θ − θ ) hb
h < hb
r
s
r
h
(4)
⎪
⎩
θ(h) = θsh ≥ hb(3)

where ks [L3β M− β ] is the distribution coefficient, and β [-] and η[L3β M− β ] are coefficients of the nonlinear Freundlich and Langmuiradsorption isotherm models, respectively. For linear adsorption β = 1and η = 0 with dimensions of [L3 M− 1]), for nonlinear Langmuiradsorption β = 1, and for nonlinear Freundlich adsorption η = 0.
Substituting Eq. (3) into Eq. (2) yields a nonlinear differential
equation that can be solved with nonlinear iterative schemes. Picard
iteration was chosen for the nonlinear CDE in the same fashion as for theRichards equation. Note that numerical solvers and preconditioners
specific for nonsymmetric matrices are required for the discretized CDE(see Section 3 below). To increase the computational efficiency, simul­
taneous multi-solute transport simulations (i.e., phosphorus andammonium) were enabled in the modified ParSWMS code – again,
similar as in the current version of HYDRUS.

( )− 2− 3λ
h
hb

(5)

3. Hydraulic and chemical soilless substrate properties

K(h) = Ksat

The selection of the three soilless substrates that were considered forthe simulations was guided by production scale greenhouse tomato
growth experiments at the Volcani Center in Israel. The substrates
included perlite, volcanic tuff, and a 70/30 vol-% volcanic tuff/coconut
coir mixture.
Horticultural perlite (Fig. 1a) is a natural amorphous volcanic glassthat is formed through the hydration of obsidian. Perlite is commonly
heated to 1000˚C, which causes structural water to evaporate and thevolume to expand to about 4 to 20 times of its original size whenrehydrated (Bar-Tal et al., 2019). The sieved perlite aggregates are lightweight, inert, and pathogen free (Noland et al., 1992).
Tuff (Fig. 1b) is a pyroclastic volcanic material with high porosity
and surface area. The mineral composition and weathering stage, in
conjunction with physical alterations (i.e., grinding and sieving),
determine its physicochemical properties (Silber et al., 1994). Tuff
possesses a high buffering capacity and may adsorb and release nutrients
throughout the growing season, especially phosphorus (Silber et al.,
1999; Silber and Raviv, 1996). It exhibits bulk densities between 0.8 and1.5 g cm− 3 and a total porosity ranging from 60 to 80%.
Coconut coir consists of short and medium-length fibers of themesocarp of Cocos nucifera L. that are left from various industrial ap­
plications. It exhibits remarkable physical and chemical characteristics
such as high water holding capacity, good drainage and aeration prop­
erties, and high cation exchange capacity (Evans et al., 1996; Abad et al.,2005). It is commonly used as a surrogate for peat moss and mixed withmineral substrates such as tuff (Fig. 1c).
The substrate water characteristics (SWC) were measured with
Tempe cells (Soilmoisture Equipment Corp., Santa Barbara, CA) con­

where h is the pressure head, hb is the air entry pressure (i.e., pressure
threshold related to the onset of drainage of the largest pore in thesystem – a transition from fully to partially saturated conditions), θ is thewater content expressed as a function of h, θs is the saturated watercontent, θr is the residual water content, λ is an empirical shapeparameter, and K is the hydraulic conductivity expressed as a function ofh.
The parameters of the BC model were determined via least-square
fitting to the measured θ − h data pairs (Fig. 2), with SWC parametersand Ksat listed in Table 1.
The phosphorus and ammonium adsorption isotherms were
measured with a combination of calorimetric spectrometry and induc­tively coupled plasma mass spectrometry. Details are provided in
Gohardoust et al. (2020). The derived isotherm parameters are listed inTable 2.
4. Verification of the modified ParSWMS code
Using HYRUS version 3.02.0560 (Šimůnek et al., 2018) as a refer­ence, we meticulously verified the modified ParSWMS code using asimple domain representative of a containerized growth module (Fig. 3).
We note that HYDRUS has been extensively tested for specific flow ortransport problems based on analytical solutions and measured data (e.
g., Cook et al., 2006; Kandelous and Šimůnek, 2010; Karlsson et al.,2015; Lassabatere et al., 2014; Luo and Sophocleous, 2010; Neumann
et al., 2011; Vanderborght et al., 2005), and code-to-code validation is ascientifically accepted approach (e.g., Greenwald, 2010; Hardelauf
et al., 2007; Kačur and Minár, 2013; Orgogozo et al., 2014).

Fig. 1. Soilless substrates considered for flow and transport simulations. (a) perlite, (b) tuff, and (c) 70/30 vol-% volcanic tuff/coconut coir mixture.
3

M.R. Gohardoust et al.

Journal of Hydrology 596 (2021) 126053

Fig. 2. The Brooks and Corey (1964) SWC model fitted to the measured data
for tuff (left), the 70/30 vol–% tuff/coconut coir mixture (center), andperlite (right).
Table 1
Parameters of the Brooks & Corey SWC model and measured Ksat .
Substrateθs(cm3 cm− 3)

θr(cm3 cm− 3)

hb
(cm)

λ
(–)

Ksat
(cm h− 1)

Perlite
Tuff
Tuff/Coconut Coir

0.818
0.483
0.549

0.001
0.014
0.014

0.458
0.382
1.102

0.149
0.247
0.219

305.1
304.2
110.7

Fig. 3. Flow domain geometry with water flow (top) and solute transport(bottom) BCs. The Cauchy-1 BC is used for the solute concentrations at the top
atmospheric and bottom free drainage boundaries. The Cauchy-2 BC is used forthe solute concentrations at the drip irrigation boundaries.

0.032, 0.063, and 0.068 cm2 h− 1, respectively. Similar values have beenreported in Buffle et al. (2007) and Hashitani and Tanaka (1983). The
water content tolerance was set to 2 × 10− 5 cm3 cm− 3, and the absoluteand relative concentration tolerances were set to 10− 6 mmol cm− 3 and10− 3, respectively. The flow and transport domains were discretizedwith finite elements of 1.4-cm in the horizontal direction and 0.7–cm inthe vertical direction. Mesh refinements for the top 1-cm layer (i.e.,
using half of the finite element size) led to a total of 25,264 spatial
nodes. The temporal discretization was dynamically calculated
throughout the simulation based on the specified initial time step of
2.4E–7 hr. The total simulation time was ten days.
Input files for the modified ParSWMS code were generated with the“Export to ParSWMS” function in HYDRUS. To extract and analyze
ParSWMS simulation results, dedicated MATLAB® – Version R2019b(MathWorks, Natick, MA, USA) scripts were created. The simulations
were performed on a workstation with two 2.40 GHz Intel® Xeon® E52630 v3 processors. The hard drive was partitioned with Windows 10
Pro (Version 1909, OS Build 18363.1016) as well as Ubuntu 18.04.5 LTS
installed to run the HYDRUS and modified ParSWMS codes,
respectively.
The modified ParSWMS and HYDRUS water flow and solute trans­
port simulation results were compared at the end of the simulations
based on the normalized root mean squared error (NRMSE) given as:√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
N∑
1(CParSWMS i − CHYDRUS i )N

Table 2
Nonlinear adsorption isotherm parameters for phosphorus and ammonium.
Substrate

Nonlinear Langmuir
Phosphorus

Perlite
Tuff
Tuff/Coconut Coir

Perlite
Tuff
Tuff/Coconut Coir

Ammoniumks(cm3 g− 1)

η(cm3 mg− 1)

ks(cm3 g− 1)

η

17.71
17.86
24.66

984
66
102

147.19
58.43
42.94

3376
135
83

Nonlinear Freundlichksβ(-)(cm3β g− β )ks(cm3β g− β )

β
(-)

9.12
39.19
23.10

27.66
85.23
72.34

0.16
0.36
0.43

0.28
0.41
0.62

(cm3 mg− 1)

To mimic irrigation and fertigation with drip emitters, water, phos­
phorus, ammonium, and nitrate were applied to the surface within twocircles with variable flux BCs. When no water or solutes were applied,the variable flux BC was treated as an atmospheric BC to conform withthe remaining evaporating substrate surface. The 70/30 vol-% tuff/co­
conut coir mixture (see Table 1 for hydraulic parameters), which wasdry (i.e., − 5000 cm pressure head) at the start of the simulation, wasused for validation. While nonlinear Langmuir adsorption was consid­
ered for ammonium (NH4 + ) and phosphorus (H2 PO4 − ) (see adsorption
isotherm parameters in Table 2), nitrate (NO3 − ) did not interact (i.e.,adsorb) with the solid phase. The applied concentrations for H2 PO4 − ,NH4 + , and NO3 − were 20 mg l− 1, 20 mg l− 1, and 80 mg l− 1, respectively.
Water was applied every 8 h at a rate of 1.0 l h− 1 for 4 min with onefourth allowed to evaporate. The longitudinal and transverse dis­persivities were assumed to be 2.0 and 0.2 cm, respectively. The diffu­
sion coefficients for H2 PO4 − , NH4 + , and NO3 − in water were assumed as

NRMSE =i=1

CHYDRUSmax − CHYDRUS min(6)

where C is the specific output quantity to be compared, and N is thenumber of spatial nodes.
The comparison of the modified ParSWMS and HYDRUS volumetric−water content (θ) and H2 PO4 − , NH+
4 , and NO3 concentration simulation
results after 10 days are depicted in Fig. 4. The coefficients of determi­−
nation for θ and H2 PO4 − , NH+
4 , and NO3 concentrations exceeded 0.999
with associated NRMSEs of less than 1%, which indicates excellent
agreement between HYDRUS and the modified ParSWMS code. This isfurther solidified by the very low mass balance errors of less than 0.3%
4

M.R. Gohardoust et al.

Journal of Hydrology 596 (2021) 126053−
Fig. 4. Comparison of the modified ParSWMS and HYDRUS volumetric water content (θ) (a) and H2 PO4 − (b), NH+
4 (c), and NO3 (d) concentration simulation results
for the 25,264 spatial nodes after 10 days.

for water and less than 2.5% for the solutes (Table 3). Notably, the massbalance errors of the modified ParSWMS code are smaller than the errorsof HYDRUS.

extended computation times, especially for large three-dimensional
problems (Kwack et al., 2016; Duff and Scott, 2004; Benzi, 2002). On
the other hand, sparse iterative solvers such as preconditioned Krylov
subspace solvers can most efficiently handle massive linear equation
systems (Dongarra and Sullivan, 2000). To resolve Eq. (7) iteratively,
the solution is approximated via constructing a sequence of xn , startingwith an initial guess x0 , by moving in an affine subspace x0 + K ⊂RN . In
order to create suitable subspaces K , methods containing efficient op­
erations such as matrix–vector products (e.g., Krylov subspaces) are ofgreat value (Del Corso et al., 2015). Assuming the residual of the nthiteration rn = b − Axn , Krylov subspace projection methods can beseparated into three main categories (Golub and van der Vorst, 1997):(1) Ritz-Galerkin methods in which xn is constructed such that the re­sidual is orthogonal to the current subspace – examples are conjugate
gradients (CG), the full orthogonalization method (FOM) (Saad, 1981),
and generalized conjugate gradients (GENCG) (Eisenstat et al., 1983);(2) minimum residual methods in which xn minimizes the second normof the residual over K n – the generalized minimal residual (GMRES)(Saad and Schultz, 1986), the minimum residual (MINRES), ORTHOMIN(Vinsome, 1976), and ORTHODIR (Young and Jea, 1980) methods fall

5. Evaluation of preconditioners and solvers
To simulate water flow and solute transport, separate sets of linear
equations need to be solved. To assure computational efficiency, it isimperative to apply the optimal preconditioner/solver combinations,which differ for flow and transport. The discretization of the associatedpartial differential equations via finite element or finite difference
methods yields linear systems of the form:(7)

Ax = bwith A ∈ RN×N is the coefficient matrix, b ∈ RN is the right-hand side
vector, and x is the vector of unknowns.
Depending on the extent of the problem, direct or iterative solvers
may be employed to solve Eq. (7). To enhance the numerical robustness
and stability of ill-posed problems with very high condition numbers,
sparse direct solvers are preferred despite their high memory usage andTable 3Cumulative water fluxes and mass balance information.

ParSWMSHYDRUS

Irrigation(cm3)

EATM *
(cm3)

EVF †
(cm3)

FD‡(cm3)

4001.864001.88

902.8
902.7

49.13
48.88

912.5
912.9

*

Water evaporation from areas with an atmospheric BC.
Water evaporation from areas with a variable flux BC when not irrigated.
‡
Water loss from the free drainage boundaries.
¶
Change in water storage.
†

5

ΔSw ¶
(cm3)

Mass Balance Error (%)
Water

H2 PO4 −

NH4 +

NO3 −

2147.412147.41

0.249
0.250

1.806
2.421

1.401
2.439

0.451
0.514

M.R. Gohardoust et al.

Journal of Hydrology 596 (2021) 126053

within this category; and (3) Petrov-Galerkin methods in which xn isdetermined so that the residual is orthogonal to some other n-dimen­
sional subspace – the biconjugate gradient (Bi-CG) (Fletcher, 1975) andquasi-minimal residual (QMR) methods (Freund and Nachtigal, 1991)
are part of this category. In addition, there are hybrids of the above
categories such as the conjugate gradient squared (CGS) (Sonneveld,1989), biconjugate gradient stabilized (Bi-CGSTAB) (van der Vorst,1992), transpose-free quasi-minimal residual (TFQMR) (Freund, 1993),
and flexible generalized minimal residual (FGMRES) (Saad, 1993)
methods.
In brief, preconditioning is applied to convert a given problem into aform that is more amenable for numerical solution (Herbst et al., 2008).
This means that Eq. (7) is transformed to the form:
M −1 1 AM −2 1 y = M −1 1 b

y = M2 x

surface drip emitters and the solution drains through circular openings
in the bottom of the concrete containers (Fig. 5b). The 70/30 vol-% tuff/
coconut coir mixture with its hydraulic properties listed in Table 1 wasused as the substrate for the test case. The applied water flow and solutetransport BCs are consistent with the flow domain depicted in Fig. 3. For
water flow, this means an atmospheric BC at the top with ten circularinclusions with variable flux BCs to mimic surface drip irrigation (i.e.,
two angle arrow drippers per tomato plant) and seven openings at thebottom to allow for free drainage. For root water uptake, a simple root
distribution model (Vrugt et al., 2001a, 2001b) was considered for eachof the five tomato plants (Šimůnek et al., 2018). While uptake conforms
with the Feddes et al. (1978) water stress response function, the po­tential uptake rate was constrained at 13% of the evapotranspiration
rate. Nutrients were applied together with irrigation water (see Section 4for concentrations). Drip irrigation occurred every two hours for 2.75
min at a rate of 1.0 L per hour per angle arrow dripper (i.e., a total of 5.5
L were applied per day).
To consider heterogeneity, independent scaling factors for the pres­sure head and hydraulic conductivity function were introduced. The
scaling factors were generated via sequential Gaussian simulation with
the VISIM program (Hansen and Mosegaard, 2008) in the open source
mGstat MATLAB® toolbox (Hansen, 2020) with a variance of 0.05, amean value of 1.0, and correlation lengths of 4 and 2 cm in the hori­zontal and vertical directions, respectively. The spatial distribution of
the hydraulic conductivity scaling factors within the flow domain is
depicted in Fig. 6.
The flow domain was discretized with finite elements of 1.0-cm inthe horizontal direction and 0.5–cm in the vertical direction. Mesh re­finements for the top 1-cm layer, and in areas with variable flux and freedrainage BCs, led to a total of 408,848 nodes and 2,315,763 3D meshelements.
To solve Eq. (7) with its symmetric sparse coefficient matrix forwater flow, a total of 15 preconditioners and 13 solvers were considered(Table 4). All preconditioners and solvers, but two, were extracted fromthe Portable, Extensible Toolkit for Scientific Computation (PETSc) library(https://www.mcs.anl.gov/petsc/). Two preconditioners were incor­
porated from the High Performance Preconditioners (HYPRE) library(http://www.llnl.gov/CASC/hypre/) via the PETSc interface. The
sparsity pattern of the coefficient matrix associated with the water flow
simulation is depicted in Fig. 7. Because of the large number ofconsidered preconditioner/solver combinations, the total time for each(8)

where M1 and M2 are nonsingular matrices with inverses that can be
rapidly and accurately solved. This leads to faster convergence and re­quires much less memory (Simoncini and Szyld, 2007).
To find optimal preconditioner/solver combinations with regard toboth robustness and computational efficiency, we consider a realistic
flow and transport simulation applicable for soilless culture growthmodules. Note that the optimal preconditioner/solver combinations
differ for water flow and solute transport simulations as they exhibit
symmetric and asymmetric coefficient matrices, respectively. Hence,
they are evaluated separately.
5.1. Preconditioners and solvers for water flow
Computational efficiency is essential when simulating flow andtransport processes in complex porous systems, where intricate hy­
draulic properties and boundary conditions demand high-resolution
discretization of the spatial and temporal domains, often leading toseveral million nodes requiring numerical evaluation. This holds true forsoilless plant growth modules, especially when numerous realizations
are desired to optimize container geometry, soilless substrate properties,
and irrigation and fertigation management.
To evaluate potential preconditioner/solver combinations for water
flow, we considered a typical growth module for greenhouse tomato
production at the Agricultural Research Organization (ARO) Volcani
Center in Israel (Fig. 5a), where water and nutrients are applied via

Fig. 5. Greenhouse tomato growth trial at the ARO Volcani Center in Israel (a). Rendering of a concrete growth module used in the greenhouse trial (b).
6

M.R. Gohardoust et al.

Journal of Hydrology 596 (2021) 126053

Fig. 6. Spatial distribution of the scaling factor (αK ) for the hydraulic con­ductivity function within the flow domain.
Table 4
Preconditioners and solvers considered for the water flow simulations.
Preconditioner

Abbr.

SolverAbbr.

JACOBI (diagonalscaling)
Block JACOBI*JACOBI

Conjugate GradientCG

BJACOBICGS

Additive Schwarz(Restrict)
Additive Schwarz(Interpolate)
Additive Schwarz(Basic)Additive Schwarz(None)
Shared Blocks ASM(Restrict)
Shared Blocks ASM(Interpolate)
Shared Blocks ASM(Basic)Shared Blocks ASM(None)
Successive Over
Relaxation†
Symmetric SOR(Eisenstat)
Classical Algebraic
Multigrid‡
Sparse Approximate
Inverse‡
Geometric Algebraic
Multigrid¶

ASM-R

Conjugate Gradient
SquaredFlexible Conjugate
Gradient
Pipelined FCG

GMRES

ratio of utilized cores and the number of cores used for the base case (i.e.
12):

ASM-N

Generalized Minimal
Residual§
Deflated GMRES

DGMRESTp =GASM-R

Loose GMRESLGMRESGASM-I

Biconjugate GradientBICG

GASM-B

BCGS

SOR

Biconjugate GradientStabilized
Improved Stabilized
BICG
Enhanced BICG

BCGSL

EISEN

Minimum Residual

MINRESBOOMERAMG

Chebyshev Iterative
MethodCHEBYSHEV

ASM-I
ASM-B

GASM-N

PARASAIL§

Fig. 7. Sparsity pattern of the symmetric coefficient matrix for water flow
simulations.

FCG
PIPEFCG

pbase
Tbase
p

(9)

where Tp is the simulation time when p processors are utilized, and pbase
and Tbase are the number of processors used for the base case and theassociated simulation time, respectively. The inset in Fig. 8 shows all
successful combinations. As apparent, the choice of the preconditioner/solver combination significantly impacts computational efficiency, as
the determined wall clock times vary by more than one order ofmagnitude. Note that some of the simulations with the GAMG pre­conditioner did not converge when 94 cores were utilized. The reported
wall clock times represent the median of 5 separate simulation runs. The
simulation speedup follows the ideal speedup line up to 48 cores,
beyond which the performance starts to degrade. This is attributable tothe increase of the communication to computation ratio (Hammondet al., 2014). For the executed simulations, the average number of nodes
for each core was about 8500 and 4400 when utilizing 48 and 94 cores,
respectively. This suggests that for optimal scalability, the minimum
number of nodes per processor should fall within this range.
Because the order of the most efficient combinations is dependent on
the number of cores (NC), the wall clock times were normalized by thesmallest value for each NC group (i.e., 12, 24, 48, and 94) and thenaveraged. The five fastest combinations and their associated normalized
times are listed in Table 5. Although the DGMRES solver applied inconjunction with the BJACOBI preconditioner performed better thanother widely used solvers such as CG, the differences are minor between
the five computationally most efficient combinations (Table 5). This
leads to the conclusion that BJACOBI, which performs the incomplete
lower–upper factorization in each block (with zero level fill), is thepreconditioner of choice for our water flow simulations. This finding in
general concurs with Herbst et al. (2008), who investigated the perfor­
mance of five preconditioners in combination with the CG solver for
simulation of one- and three-dimensional water flow (i.e., Richards
equation) in porous media. In contrast to their results, the algebraic
multigrid preconditioner (BOOMERAMG) did not work well for our

IBCGS

GAMG

*
Inner preconditioner for each block is ILU(0), i.e., incomplete lower–upper
factorization with zero level fill.
†
Equivalent to the block Jacobi with SOR on each block.
‡
From the HYPRE library. BOOMERAMG with “HMIS” coarsen type, “ext + I”
interpolation type, zero number of levels of aggressive, and 0.5 threshold value
for being strongly connected.
¶
GAMG with type ‘aggregate’ and one smoothing step.
§
With restart parameter 60.

simulation was limited to 10 h. The simulations were performed on the
University of Arizona Puma HPC cluster with system specifications
provided in Table S1. The number of cores (NC) utilized were 12, 24, 48,
and 94, respectively.
To compare the computational efficiency of various preconditioner/solver combinations, the wall clock time, which includes preconditioner
setup and application time, as well as the time for solver application,
was used as a metric. Fig. 8 depicts the ten fastest preconditioner/solver
combinations together with the ideal speedup line that is defined as the7M.R. Gohardoust et al.

Journal of Hydrology 596 (2021) 126053

Fig. 8. Wall clock times for the ten most efficient preconditioner/solver combinations for simulation of water flow. The five most efficient combinations are rep­resented by solid lines and symbols, the remaining five by dashed lines. The inset depicts all successful preconditioner/solver combinations.

Table 5The most efficient preconditioner/solver combinations for the water flow
simulations.
Rank

Preconditioner/Solver Combination

Normalized Time

1
2
3
4
5

BJACOBI/DGMRESBJACOBI/GMRESBJACOBI/LGMRESBJACOBI/CG
SOR/DGMRES

1.01281.02761.04501.07041.1407

Table 6Preconditioners and solvers considered for the solute transport simulations.

simulations, which is most likely because of the more complex water
distributions within our flow domain due to the implementation of
surface drip irrigation, root water uptake, and free drainage from thebottom boundary.
5.2. Preconditioners and solvers for solute transport
The test case for solute transport is identical to the test case in Section
5.1 with regard to water flow. In addition, H2 PO4 − , NH4 + , and NO3 −were applied with the irrigation water at concentrations of 20 mg l− 1,
20 mg l− 1, and 80 mg l− 1, respectively. While nonlinear Langmuiradsorption was considered for ammonium (NH4 + ) and phosphorus(H2 PO4 − ), nitrate (NO3 − ) did not interact with the solid phase. Fifteen
preconditioners and 12 solvers were evaluated (Table 6).
The ten computationally most efficient preconditioner/solver com­binations for solute transport are depicted in Fig. 9 together with theideal speedup line. The inset shows all successful combinations, whichvary over two orders of magnitude. As for water flow, the simulation
speedup follows the ideal speedup line up to 48 cores and then theperformance starts degrading. Combinations with the GAMG precondi­
tioner were the least efficient. Although the BOOMERAMG precondi­
tioner performed better for solute transport, it is not among the topranked. These results differ from Sbai and Larabi (2021), who simu­
lated field scale chemical transport in groundwater. The deviations canbe attributed to the consideration of adsorbing solutes in our simulations
as well as the generally better performance of the algebraic multigrid
preconditioner for large transport domains (Sbai and Larabi, 2020). The
most efficient combinations consisted of preconditioners of the succes­
sive over relaxation category (i.e., SOR and EISENSTAT) and variationsof the GMRES solver. It should be noted that the parallel SOR and
EISENSTAT preconditioners are equivalent to BJACOBI with SOR on

Preconditioner

Abbr.

SolverJACOBI (diagonalscaling)
Block JACOBI*
Additive Schwarz(Restrict)
Additive Schwarz(Interpolate)
Additive Schwarz(Basic)Additive Schwarz(None)
Shared Blocks ASM(Restrict)
Shared Blocks ASM(Interpolate)
Shared Blocks ASM(Basic)Shared Blocks ASM(None)
Successive Over
Relaxation†
Symmetric SOR(Eisenstat)
Classical Algebraic
Multigrid‡
Sparse Approximate
Inverse‡
Geometric Algebraic
Multigrid¶

Abbr.

JACOBI

Transpose Free QMR

§

TFQMR

BJACOBI
ASM-R

TCQMR
GMRES

ASM-I

A Variant of QMRGeneralized Minimal
Residual||
Flexible GMRES||

ASM-B

Deflated GMRES

DGMRES

ASM-N

Pipelined GMRES

PGMRESGASM-R

Pipelined Flexible
GMRES
LGMRES||

PIPEGMRESGASM-IGASM-BGASM-NSOR
EISENBOOMERAMG

Conjugate Gradient
Squared
Stabilized Biconjugate
Gradient
Enhanced BICG
Improved Stabilized
BICG||

FGMRES

LGMRES
CGS
BCGS
BCGSL
IBCGS

PARASAILGAMG*

Inner preconditioner for each block is ILU(0), i.e., incomplete lower–upper
factorization with zero level fill.
†
Equivalent to BJACOBI with SOR on each block.
‡
From the HYPRE library. The BOOMERAMG with “HMIS” coarsen type, “ext
+ i“ interpolation type, zero number of levels of aggressive, and 0.5 thresholdvalue for being strongly connected.
¶
GAMG with type ‘aggregate’ and zero smoothing steps.
§
QMR: Quasi Minimal Residual method.
||
With restart parameter 60.

each block and are therefore not true parallel preconditioners. As
evident from Fig. 9 and Table 7, the SOR/FGMRES combination is best
suited for our application.

8M.R. Gohardoust et al.

Journal of Hydrology 596 (2021) 126053

Fig. 9. Wall clock times for the ten most efficient preconditioner/solver combinations for simulation of solute transport. The five most efficient combinations are
represented by solid lines and symbols, the remaining five by dashed lines. The inset depicts all successful preconditioner/solver combinations.

Table 1). The simulation duration was 14 days (336 hrs). The flow
domain was discretized with finite elements of 1.90-cm in the horizontal
direction and 0.95–cm in the vertical direction. Mesh refinements for the
top 1-cm layer and in areas with variable flux and free drainage BCs, led
to a total of 132,570 nodes and 726,264 3D mesh elements. Note that thehydraulic characteristics are calculated directly from the hydraulicfunctions by setting the lower and upper limits of the internal interpo­
lation tables to zero in both HYDRUS and the modified ParSWMS code.
The BJACOBI/DGMRES and SOR/FGMRES preconditioner/solver com­binations were used for ParSWMS water flow and solute transport
simulations, respectively.
Fig. 10 shows the temporal evolution of mass balance errors forwater flow and H2 PO4 − , NH4 + , and NO3 − transport. The water massbalance errors are below 0.15% for all substrates for both the HYDRUS
and ParSWMS simulations, which we consider quite accurate. The massbalance errors for solute transport are higher for both HYDRUS and themodified ParSWMS. An extreme error of about 80% for NO3 − is evident
for the HYDRUS simulation for tuff, which can be attributed to simu­
lation instability issues starting at around 48 hrs. This is likely due tovery localized increases in NO3 − concentrations at the atmospheric BC
because of evaporation (Fig. 11, bottom right). In comparison, themaximum mass balance error for NO3 − and tuff simulated with themodified ParSWMS is 0.54%. Note that the instability issues experienced
with HYDRUS can be potentially resolved by increasing the density of
the 3D mesh close to the surface at the cost of computational efficiency.
Notably, the ParSWMS solute mass balance errors are significantly lower
than the HYDRUS errors.
Fig. 12 shows a comparison of computational efficiencies of the
HYDRUS and modified ParSWMS codes. Besides the simulations per­formed on the workstation, we ran ParSWMS on the Ocelote HPC (28cores) and the recently installed Puma HPC (94 cores) clusters (forspecifications see Table S1) to demonstrate the massive increase inscalability when HPC environments are utilized. The simulation speedup
expressed in terms of normalized wall clock times (see values on top of
the bars in Fig. 12) of the modified ParSWMS code is about 22% fasterthan that of HYDRUS on the workstation, and about 94% faster on the
Puma HPC cluster with 94 cores. Note that while an increase in thenumber of processing cores generally speeds up the simulations, beyonda certain core threshold, the computational efficiency starts to decline
due to an increase of the communication to computation ratio.

Table 7The most efficient preconditioner/solver combinations for the solute transport.
Ranking

Preconditioner/Solver Combination

Normalized Time

1
2
3
4
5

SOR/FGMRES
EISENSTAT/DGMRES
EISENSTAT/GMRES
EISENSTAT/LGMRESJACOBI/FGMRES

1.00561.01861.0312
1.04211.0810

6. Computational efficiency and stability of HYDRUS and themodified ParSWMS
As already indicated in Section 5, computational efficiency and sta­bility are essential when simulating flow and transport processes insoilless plant growth modules, especially when a vast number of re­alizations is required. Hence, it is of interest to compare the Windows OS
based HYDRUS with the Linux OS based modified ParSWMS. Because
workstations and HPC clusters are not intercomparable, such compari­son is only viable on a workstation, with both Windows (i.e., Windows
10 Pro Version 1909, OS Build 18363.1016) and Linux (i.e., Ubuntu
18.04.5 LTS) installed on a partitioned hard drive. A workstation withtwo 2.40 GHz Intel® Xeon® E5-2630 v3 Dual 8-Core processors and 32
GB RAM was used for the comparison.
HYDRUS benefits from a powerful GUI with numerous advanced
features such as the ability to create or import complex domain geom­etries or to specify various spatially variable properties such as mate­rials, initial conditions, boundary conditions, and domain properties
either directly on the finite element mesh or on geometric objects in­dependent of the mesh (Šimůnek et al., 2016). To utilize parallel
computing and the computing prowess of the latest multicore worksta­
tions, HYDRUS relies on the HYdrus PARallelized (HYPAR) module andthe Microsoft Parallel Patterns Library to numerically solve Eq. (7) forwater flow and solute transport with the CG and ORTHOMIN solvers,
respectively. Equation (7) is preconditioned with the incomplete lower
upper factorization method. Both the solvers and the preconditioner are
from the ORTHOFEM library (Mendoza et al., 1991).
The test case for the stability and efficiency comparison conforms
with the descriptions in Sections 5.1 and 5.2, except that the SWC andhydraulic conductivity functions were not scaled for the considered tuff,
tuff/coconut coir, and perlite substrates (see hydraulic properties in9M.R. Gohardoust et al.

Journal of Hydrology 596 (2021) 126053

Fig. 10. Temporal mass balance error evolutions for water flow (left) and H2 PO4 − , NH4 + , and NO3 − transport (right) simulated with both HYDRUS and the modified
ParSWMS code.

Fig. 11. Water content and NO3 − concentration distributions at the end of the HYDRUS and modified ParSWMS simulations. Note the localized very high NO3 −concentrations depicted in the bottom right graph for HYDRUS.

7. Conclusions

multi-solute transport simulations, but also the application of boundary
conditions to realistically represent typical soilless culture systems. The
modified ParSWMS code was thoroughly tested using HYDRUS as areference. To optimize the computational efficiency of the modified
ParSWMS, numerous preconditioner/numerical solver combinations
were tested for both the water flow and solute transport equations.

With the goal of performing a vast number of flow and transport
simulations on HPC clusters to aid the design and management of soil­
less culture growth modules, we modified the open source 3D ParSWMS
code to not only enable nonlinear equilibrium solute adsorption and
10

M.R. Gohardoust et al.

Journal of Hydrology 596 (2021) 126053

CRediT authorship contribution statement
Mohammad R. Gohardoust: Conceptualization, Methodology,
Software, Formal analysis, Investigation, Writing - original draft,
Writing - review & editing. Jirka Šimůnek: Conceptualization, Meth­odology, Writing - review & editing. Horst Hardelauf: Conceptualiza­
tion, Software, Writing - review & editing. Markus Tuller:
Conceptualization, Methodology, Formal analysis, Writing - originaldraft, Writing - review & editing, Supervision, Funding acquisition.
Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Acknowledgments
This research was funded by The United States - Israel Binational
Agricultural Research and Development Fund (BARD), grant number
US-4764-14 R, and by the United States Department of Agriculture(USDA) National Institute of Food and Agriculture (NIFA), Hatch/MultiState project number ARZT-1370600-R21-189. We thank Dr. Asher BarTal from the ARO Volcani Center in Israel for providing details about thedesign and management of the tomato growth experiment that was
simulated in this study.

Fig. 12. Comparison of computational efficiencies of the HYDRUS and themodified ParSWMS codes for perlite (left), tuff (center), and tuff/coconut coir(right). The values on top of the bars represent the clock times normalized with
the HYDRUS clock times.

Appendix A. Supplementary data

While the BJACOBI/DGMRES preconditioner/solver combination was
the most efficient for water flow, SOR/FGMRES worked best for solutetransport. In general, GMRES-type and CG solvers applied in conjunc­tion with the BJACOBI preconditioner were very efficient for water flow
simulations. The most efficient combinations for solute transport con­sisted of preconditioners of the successive over relaxation category (i.e.,
SOR and EISENSTAT) and variations of the GMRES solver. To compare
the stability and computational efficiency of modified ParSWMS with
HYDRUS, a real growth module from a production scale greenhouse
experiment at the ARO Volcani Center in Israel was simulated consid­
ering three soilless substrates and H2 PO4 − , NH4 + , and NO3 − transport.
The simulations were performed on a workstation and two HPC clusters.
The results revealed that the modified ParSWMS was about 22% more
efficient than HYDRUSwhen the simulations on the workstation werecompared. Simulations on the HPC clusters were up to 94% more effi­
cient when 94 cores were utilized. If the full computational power of the
HPC clusters (i.e., no competing projects and up to 23,616 available
cores) was to be utilized, simulation times would again be drasticallyreduced to a very small fraction of the HYDRUS wall clock time. While
we experienced some stability issues for tuff and NO3 − with HYDRUS,
leading to an unrealistically high mass balance error, all modified
ParSWMS simulations were stable on both the workstation and HPC
clusters. In general, all water and solute mass balance errors generated
with the modified ParSWMS where bellow the mass balance errors of
HYDRUS.
The availability of the modified and tested ParSWMS code opens new
avenues for streamlining 3D water flow and solute transport simulations
for complex porous media via the utilization of HPC environments.
Simulations that took days to complete can now be resolved within amatter of a few hours, allowing for a vast number of realizations within ashort period of time.

Supplementary data to this article can be found online at https://doi.
org/10.1016/j.jhydrol.2021.126053.
References
Abad, M., Fornes, F., Carrión, C., Noguera, V., Noguera, P., Maquieira, Á ngel,
Puchades, R., 2005. Physical properties of various coconut coir dusts compared topeat. HortScience 40 (7), 2138–2144.
Bar-Tal, A., Saha, U.K., Raviv, M., Tuller, M., 2019. Inorganic and Synthetic Organic
Components of Soilless Culture and Potting Mixtures. In: Raviv, M., Lieth, H.J., BarTal, A. (Eds.), Soilless Culture – Theory and Practice, 2nd Ed. Academic Press(Elsevier), London, pp. 264–267.
Benzi, M., 2002. Preconditioning techniques for large linear systems: a survey.
J. Comput. Phys. 182 (2), 418–477. https://doi.org/10.1006/jcph.2002.7176.
Brooks, R.H., Corey, A.T., 1966. Properties of porous media affecting fluid flow. J. Irrig.
Drain. Div., Am. Soc Civ. Eng. 92 (2), 61–88.
Buffle, J., Zhang, Z., Startchev, K., 2007. Metal Flux and Dynamic Speciation at (Bio)interfaces. Part I: Critical Evaluation and Compilation of Physicochemical
Parameters for Complexes with Simple Ligands and Fulvic/Humic Substances.
Environ. Sci. Technol. 41, 7609–7620. https://doi.org/10.1021/es070702p.
Cook, F.J., Fitch, P., Thorburn, P.J., Charlesworth, P.B., Bristow, K.L., 2006. Modelling
trickle irrigation: comparison of analytical and numerical models for estimation ofwetting front position with time. Environ. Model. Softw. 21 (9), 1353–1359. https://doi.org/10.1016/j.envsoft.2005.04.018.
Del Corso, G.M., Menchi, O., Romani, F., 2015. Krylov Subspace Methods for Solving
Linear Systems. Technical Report. University of Pisa, Department of Computer
Science, Pisa, Italy.
Dongarra, J., Sullivan, F., 2000. Guest editors’ introduction: the top 10 algorithms.
Comp. Sci. Eng. 2 (1), 22–23. https://doi.org/10.1109/MCISE.2000.814652.
Duff, I.S., Scott, J.A., 2004. A parallel direct solver for large sparse highly unsymmetriclinear systems. ACM T. Math. Software 30 (2), 95–117. https://doi.org/10.1145/
992200.992201.
Eisenstat, S.C., Elman, H.C., Schultz, M.H., 1983. Variational iterative methods fornonsymmetric systems of linear equations. SIAM J. Numer. Anal. 20 (2), 345–357.
Evans, M.R., Konduru, S., Stamps, R.H., 1996. Source variation in physical and chemical
properties of coconut coir dust. HortScience 31 (6), 965–967.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field Water Use and Crop
Yield. John Wiley & Sons, New York.
Fletcher, R., 1975. Conjugate Gradient Methods for Indefinite Systems. In: Watson, G.A.
(Ed.), Numerical Analysis. Springer Verlag, Berlin, pp. 73–89.
Freund, R.W., Nachtigal, N.M., 1991. QMR: a quasi-minimal residual method for nonhermitian linear systems. Numer. Math 60 (1), 315–339.
Freund, R.W., 1993. A transpose-free quasi-minimal residual algorithm for nonhermitian linear systems. SIAM J. Sci. Comput. 14 (2), 470–482. https://doi.org/
10.1137/0914029.
Gohardoust, M.R., Bar-Tal, A., Effati, M., Tuller, M., 2020. Characterization ofphysicochemical and hydraulic properties of organic and mineral soilless culture

8. Availability of software
All computer codes and Matlab scripts created for this project areavailable from the corresponding author upon reasonable request.

11

M.R. Gohardoust et al.

Journal of Hydrology 596 (2021) 126053
Raviv, M., Lieth, H.J., Bar-Tal, A., 2019. Significance of Soilless Culture in Agriculture.
In: Raviv, M., Lieth, H.J., Bar-Tal, A. (Eds.), Soilless Culture – Theory and Practice,
2nd Ed. Academic Press (Elsevier), London, pp. 3–14.
Saad, Y., Schultz, M.H., 1986. GMRES: A Generalized Minimal Residual Algorithm for
Solving Nonsymmetric Linear Systems. SIAM J. Sci. Stat. Comput. 7 (3), 856–869.
Saad, Y., 1981. Krylov subspace methods for solving large unsymmetric linear systems.
Math. Comput. 37 (155), 105.
Saad, Y., 1993. A flexible inner-outer preconditioned GMRES algorithm. SIAM J. Sci.
Comput. 14 (2), 461–469. https://doi.org/10.1137/0914028.
Savvas, D., 2003. Hydroponics: a modern technology supporting the application of
integrated crop management in greenhouse. J. Food Agri. Environ. 1, 80–86.
Savvas, D., Gruda, N., 2018. Application of soilless culture technologies in the modern
greenhouse industry – a review. Eur. J. Hortic. Sci. 83 (5), 280–293.
Sbai, M.A., Larabi, A., 2021. On solving groundwater flow and transport models withalgebraic multigrid preconditioning. Groundwater 59 (1), 100–108. https://doi.org/
10.1111/gwat.v59.110.1111/gwat.13016.
Silber, A., Raviv, M., 1996. Effects on chemical surface properties of tuff by growing rose
plants. Plant Soil 186 (2), 353–360.
Silber, A., Bar-Yosef, B., Chen, Y., 1999. pH-dependent kinetics of tuff dissolution.
Geoderma 93 (1–2), 125–140. https://doi.org/10.1016/S0016-7061(99)00048-8.
Silber, A., Bar-Yosef, B., Singer, A., Chen, Y., 1994. Mineralogical and chemicalcomposition of three tuffs from Northern Israel. Geoderma 63 (2), 123–144. https://doi.org/10.1016/0016-7061(94)90002-7.
Simoncini, V., Szyld, D.B., 2007. Recent computational developments in Krylov subspacemethods for linear systems. Numer. Linear Algebra Appl. 14 (1), 1–59.
Šimůnek, J., Huang, K., & van Genuchten, M.Th. (1995). The SWMS-3D Code forSimulating Water Flow and Solute Transport in Three-Dimensional VariablySaturated Media. Riverside: USDA ARS U.S. Salinity Laboratory, (Research Report
No. 139).
Šimůnek, J., van Genuchten, M.Th., & Šejna, M. (2016). Recent Developments andApplications of the HYDRUS Computer Software Packages. Vadose Zone J., 15(7).
https://doi.org/10.2136/vzj2016.04.0033.
Šimůnek, J., Šejna, M., van Genuchten, M.T., 2018. New features of version 3 of the
HYDRUS (2D/3D) computer software package. J. Hydrol. Hydromech. 66 (2),
133–142. https://doi.org/10.1515/johh-2017-0050.
Sonneveld, P., 1989. CGS, A Fast Lanczos-Type Solver for Nonsymmetric Linear Systems.
SIAM J. Sci. and Stat. Comput. 10 (1), 36–52. https://doi.org/10.1137/0910004.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898.
van der Vorst, H.A., 1992. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG
for the solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 13 (2),631–644. https://doi.org/10.1137/0913035.
Vanderborght, J., Kasteel, R., Herbst, M., Javaux, M., Thiery, D., Vanclooster, M.,
Mouvet, C., Vereecken, H., 2005. A set of analytical benchmarks to test numericalmodels of flow and transport in soils. Vadose Zone J. 4 (1), 206–221. https://doi.
org/10.2113/4.1.206.
Vinsome, P.K.W. (1976, February). Orthomin, an Iterative Method for Solving Sparse
Sets of Simultaneous Linear Equations. Conference paper, Society of Petroleum
Engineers (SPE) Symposium on Numerical Simulation of Reservoir Performance, Los
Angeles, CA. https://doi.org/10.2118/5729-MS.
Vogel, T., Cislerova, M., 1988. On the reliability of unsaturated hydraulic conductivity
calculated from the moisture retention curve. Transp. Porous Media 3 (1), 1–15.
Vogel, T., van Genuchten, M.T., Císlerová, M., 2000. Effect of the Shape of the Soil
Hydraulic Functions Near Saturation on Variably-Saturated Flow Predictions. Adv.
Water Resour. 24 (2), 133–144. https://doi.org/10.1016/S0309-1708(00)00037-3.
Vrugt, J.A., Hopmans, J.W., Šimunek, J., 2001a. Calibration of a Two-Dimensional Root
Water Uptake Model. Soil Sci. Soc. Am. J. 65 (4), 1027–1037.
Vrugt, J.A., van Wijk, M.T., Hopmans, J.W., Šimunek, J., 2001b. One-, Two-, and ThreeDimensional Root Water Uptake Functions for Transient Modeling. Water Resour.
Res. 37 (10), 2457–2470.
Young, D.M., Jea, K.C., 1980. Generalized Conjugate-Gradient Acceleration of
Nonsymmetrizable Iterative Methods. Linear Algebra Appl. 34, 159–194.

substrates and mixtures. Agronomy 10 (9), 1403. https://doi.org/10.3390/
agronomy10091403.
Golub, G.H., van der Vorst, H.A., 1997. Closer to the solutions: iterative linear solvers. In:
Duff, I.S., Watson, G.A. (Eds.), The State of the Art in Numerical Analysis. Clarendon
Press, Oxford, pp. 63–92.
Greenwald, M., 2010. Verification and validation for magnetic fusion. Phys. Plasmas 17(5), 058101. https://doi.org/10.1063/1.3298884.
Hammond, G.E., Lichtner, P.C., Mills, R.T., 2014. Evaluating the performance of parallelsubsurface simulators: an illustrative example with PFLOTRAN. Water Resour. Res.
50 (1), 208–228. https://doi.org/10.1002/wrcr.v50.110.1002/2012WR013483.
Hansen, T.M., Mosegaard, K., 2008. VISIM: sequential simulation for linear inverse
problems. Comput. Geosci. 34 (1), 53–76. https://doi.org/10.1016/j.
cageo.2007.02.003.
Hansen, T.M. (2020). mGstat: A Geostatistical Matlab toolbox. http://mgstat.
sourceforge.net/ Accessed 4 September 2020.
Hardelauf, H., Javaux, M., Herbst, M., Gottschalk, S., Kasteel, R., Vanderborght, J.,
Vereecken, H., 2007. PARSWMS: a parallelized model for simulating threedimensional water flow and solute transport in variably saturated soils. Vadose Zone
J. 6 (2), 255–259. https://doi.org/10.2136/vzj2006.0156.
Hashitani, T., Tanaka, K., 1983. ChemInform Abstract: measurements of self-diffusion
coefficients of the nitrate ion in aqueous solutions of potassium nitrate and calcium
nitrate. Chem. Informationsd. 14 (46) https://doi.org/10.1002/chin.
v14.4610.1002/chin.198346019.
Herbst, M., Gottschalk, S., Reißel, M., Hardelauf, H., Kasteel, R., Javaux, M.,
Vanderborght, J., Vereecken, H., 2008. On preconditioning for a parallel solution of
the Richards equation. Comput. Geosci. 34 (12), 1958–1963. https://doi.org/
10.1016/j.cageo.2008.02.020.
Kačur, J., Minár, J., 2013. A benchmark solution for infiltration and adsorption ofpolluted water into unsaturated-saturated porous media. Transp. Porous Media 97(2), 223–239. https://doi.org/10.1007/s11242-012-0119-5.
Karlsson, S.C., Langergraber, G., Pell, M., Dalahmeh, S., Vinnerås, B., Jönsson, H., 2015.
Simulation and verification of hydraulic properties and organic matter degradationin sand filters for greywater treatment. Water Sci. Technol. 71, 426–433. https://doi.
org/10.2166/wst.2015.003.
Kandelous, M.M., Šimůnek, J., 2010. Comparison of numerical, analytical, and empiricalmodels to estimate wetting patterns for surface and subsurface drip irrigation. Irrig.
Sci. 28 (5), 435–444. https://doi.org/10.1007/s00271-009-0205-9.
Kwack, JH., Bauer, G.H., & Koric, S. (2016, May). Performance Test of Parallel Linear
Equation Solvers on Blue Waters – Cray XE6/XK7 system. Conference paper,
Proceedings of the Cray User Group Meeting (CUG 2016), London, England.
Lassabatere, L., Yilmaz, D., Peyrard, X., Peyneau, P.E., Lenoir, T., Šimůnek, J., AnguloJaramillo, R. (2014). New Analytical Model for Cumulative Infiltration into DualPermeability Soils. Vadose Zone J., 13. https://doi.org/10.2136/vzj2013.10.0181.
Luo, Y.i., Sophocleous, M., 2010. Seasonal groundwater contribution to crop-water use
assessed with lysimeter observations and model simulations. J. Hydrol. 389 (3-4),
325–335. https://doi.org/10.1016/j.jhydrol.2010.06.011.
Mendoza, C.A., Therrien, R., Sudicky, E.A., 1991. ORTHOFEM User’s Guide, Version
1.02. Waterloo Centre for Groundwater Research, University of Waterloo, Waterloo,
Ontario, Canada.
Noland, D.A., Spomer, L.A., Williams, D.J., 1992. Evaluation of pumice as a perlite
substitute for container soil physical amendment. Commun. Soil Sci. Plan. 23(13–14), 1533–1547. https://doi.org/10.1080/00103629209368685.
Neuman, S. P., Feddes, R. A., & Bresler, E. (1974). Finite element simulation of flow insaturated-unsaturated soils considering water uptake by plants, Third Annual
Report, Project No. A10-SWC-77, Hydraulic Engineering Lab., Technion, Haifa,
Israel.
Neumann, L.E., Šimůnek, J., Cook, F.J., 2011. Implementation of quadratic upstream
interpolation schemes for solute transport into HYDRUS-1D. Environ. Model. Softw.
26 (11), 1298–1308. https://doi.org/10.1016/j.envsoft.2011.05.010.
Orgogozo, L., Renon, N., Soulaine, C., Hénon, F., Tomer, S.K., Labat, D., Pokrovsky, O.S.,
Sekhar, M., Ababou, R., Quintard, M., 2014. An open source massively parallel
solver for Richards equation: mechanistic modelling of water fluxes at the watershedscale. Comput. Phys. Commun. 185 (12), 3358–3371. https://doi.org/10.1016/j.
cpc.2014.08.004.

12