Published online February 15, 2018

Special Section: The Root
Zone: Soil Physics and Beyond

Implementation and Application ofa Root Growth Module in HYDRUS

Anne Hartmann,* Jiří Šimůnek, Moses Kwame Aidoo,
Sabine J. Seidel, and Naftali Lazarovitch

Core Ideas• A HYDRUS module was implemented
to model root growth as a functionof environmental stresses.
• Simulations were compared withexperimental data for temperaturestress effects on root growth.
• Temperature-dependent root
growth parameters were fitted toexperimental data using DREAM.
• A sensitivity analysis of model
parameters revealed key parameters of the modeling approach.

A. Hartmann, Research & Development, UGT Environmental Measure­­ment
Devices GmbH, Eberswalder Str. 58,
Muencheberg
15374,
Germany;
J. Šimůnek, Dep. of Environmental
Sciences, Univ. of California, Riverside,
CA 92521; M.K. Aidoo and N. Lazarovitch, French Associates Institutefor Agriculture and Biotechnology of
Drylands, Jacob Blaustein Institutes for
Desert Research, Ben-Gurion Univ. of
the Negev, Sede Boqer Campus 84990,
Israel; S.J. Seidel, Institute of Crop
Science and Resource Conservation,
Univ. of Bonn, Katzenburgweg 5, 53115
Bonn, Germany. *Corresponding author(hartmann.anne.de@gmail.com).

Received 6 Feb. 2017.
Accepted 13 June 2017.
Citation: Hartmann, A., J. Šimůnek,M.K. Aidoo, S.J. Seidel, and N. Lazarovitch. 2017. Implementation andapplication of a root growth modulein HYDRUS. Vadose Zone J. 17:170040.
doi:10.2136/vzj2017.02.0040© Soil Science Society of America.
This is an open access article distributed
under the CC BY-NC-ND license(http://creativecommons.org/licenses/
by-nc-nd/4.0/).

A root growth module was adapted and implemented into the HYDRUS
software packages to model root growth as a function of different environmental stresses. The model assumes that various environmental factors,
as well as soil hydraulic properties, can influence root development undersuboptimal conditions. The implementation of growth and stress functions in
the HYDRUS software opens the opportunity to derive parameters of thesefunctions from laboratory or field experimental data using inverse modeling.
One of the most important environmental factors influencing root growth issoil temperature. The effects of temperature in the root growth module was
the first part of the newly developed HYDRUS add-on to be validated bycomparing modeling results with measured rooting depths in an aeroponic
experimental system with bell pepper (Capsicum annuum L.). The experiment was conducted at root zone temperatures of 7, 17, and 27°C. Inverse
optimization was used to estimate a single set of parameters that was foundto well reproduce measured time series of rooting depths for all temperature
treatments. A sensitivity analysis showed that parameters such as the maximum rooting depth and cardinal temperatures had only a small impact on
the model output and can thus be specified using values from the literature
without significantly increasing prediction uncertainties. On the other hand,parameters that define the growth rate or the shape of the temperaturestress function had a high influence. The root growth module that considers
temperature stress only slightly increased the complexity of the standard
HYDRUS models.

Water uptake by root systems can greatly affect water flow through the soil
(Hao et al., 2005; Yu et al., 2007). The spatial pattern of root water uptake is determinedby the spatial distribution of the root system, the knowledge of which is essential for predicting the spatial distribution of water contents and water fluxes in soils. The spatial
distribution of roots and their growth are sensitive to various physical, chemical, and biological factors, as well as to soil hydraulic properties that influence the availability of waterand oxygen for plants. It is thus important to describe root growth under the influence ofvarious environmental factors to accurately simulate agricultural systems.
Various attempts have been made in the past to develop root growth models that accountfor the influence of various environmental factors such as temperature (Stone et al., 1983;
Stone and Taylor, 1983; Williams et al., 1989; Jones et al., 1991; Kaspar and Bland, 1992;
Bingham and Wu, 2011), aeration (Williams et al., 1989; Jones et al., 1991; Asseng et al.,1997), soil water availability (Jones et al., 1991; Clausnitzer and Hopmans, 1994; Sommaet al., 1998; Tsutsumi et al., 2003; Leitner et al., 2010), and soil strength (Dexter, 1986;
Williams et al., 1989; Jones et al., 1991; Grant, 1993; Asseng et al., 1997). Additional
studies evaluating the impact of various environmental factors on root growth werelisted by Wu et al. (2005). Existing models are either complex, three-dimensional rootarchitecture models (Clausnitzer and Hopmans, 1994; Somma et al., 1998; Leitner etal., 2010; Bingham and Wu, 2011; Couvreur et al., 2012) or simpler root growth models
that are implemented within more complex models such as EPIC (Williams et al., 1989)
or CERES (Robertson et al., 1993). The model by Jones et al. (1991) considers theinfluence of various environmental factors on root growth in terms of the root length

Vadose Zone Journal | Advancing Critical Zone Scienceper unit depth and as a result was found to be highly suitable for
implementation into HYDRUS.
One of the most important environmental factors influencing root
growth is soil temperature (Brouwer and de Wit, 1968; Aidoo etal., 2016). Extreme (both high and low) root zone temperatures
may influence root growth and root development and affect boththe fundamental functions of root nutrient and water uptake andtolerance of elemental deficiency and toxicity (Wu and Cheng,2014). The ability of plants to absorb and accumulate both waterand mineral nutrients from the soil is related to their capacity todevelop an extensive root system (Taiz and Zeiger, 2002) to reachlarge parts of the soil profile (Garnett et al., 2009). In addition,it has been reported that salinity and mycorrhiza, coupled withextreme root zone temperatures, affect several morphologicaland anatomical changes such as a decrease in root biomass, theroot specific surface area, the cortex width, and an altered xylem
vessel density in different species (Valenzuela-Estrada et al., 2008;
Rewald et al., 2012). These changes have been associated withextreme temperatures, freezing, and cold tolerance (Fowler et al.,1981). An overview of studies evaluating the impact of soil temperature on root growth for several crops such as maize (Zea mays L.),
rapeseed (Brassica napus L.), wheat (Triticum aestivum L.), barley(Hordeum vulgare L.), rye (Secale cereale L.), oat (Avena sativa L.),and soybean [Glycine max (L.) Merr.] was provided by Kaspar andBland (1992) and Smith et al. (2005). These studies have indicatedthat non-optimal root zone temperatures affect root growth androot development, resulting in a reduction in the rooting depth
and the formation of smaller root systems. Especially with regard
to climate change scenario analyses, there is a demand for improved
modeling predictions of root growth under heat stress or heat stresscombined with further environmental stresses (Seidel et al., 2016).
The HYDRUS software packages (Šimůnek et al., 2016) are widelyused numerical models to simulate water flow and heat and solutetransport in one-, two-, and three-dimensional, variably saturated
porous media. The governing water flow equation that is solvedby the HYDRUS models incorporates a macroscopic sink term toaccount for root water uptake (Feddes et al., 1978), which may be
reduced from its positive values due to salinity and drought stresses.
However, the macroscopic root water uptake model in the standard
versions of HYDRUS (Šimůnek et al., 2008, 2016) does not consider the feedback of environmental conditions in the root zone
on root growth, and root growth is fully defined using predefined
input parameters (Šimůnek and Hopmans, 2009).
The objectives of this study thus were (i) to develop an add-on
root growth module for the HYDRUS software packages (bothHYDRUS-1D and HYDRUS-2D) to model root growth as afunction of different environmental stresses such as temperature,
aeration, and chemical soil conditions, and (ii) to use this newly
developed add-on module to evaluate the effects of temperature on
root growth using experimental data and to carry out a sensitivity

VZJ | Advancing Critical Zone Science

analysis. The study was conducted as follows. First, the modeling approach developed by Jones et al. (1991) was adapted andimplemented into the HYDRUS software packages. The originalmodeling approach by Jones et al. (1991) is a layered root growthmodel that needs to be coupled to a crop growth model, which provides it with the daily allocation of dry matter to the root system,and that furthermore needs information about the soil profile
characteristics influencing root growth. To overcome the need for
input variables provided by a crop model, the modeling approachis instead combined with time-dependent root growth functions.
Second, the capability of the newly developed root growth module
to account for the effects of temperature on the vertical root penetration was evaluated using experimental data. The evaluation
was performed using measured rooting depths with bell pepper
in an aeroponic experimental system, which guarantees that thetemperature-dependent root growth approach was evaluated using
experimental root growth data generated only under the influenceof temperature while limiting the effects of other factors. Third, aglobal sensitivity analysis using the Sobol¢ method (Sobol¢, 1993)
was conducted to identify the key parameters of the root growthmodule and to indicate which parameters can be fitted using experimental data. The parameters to which the model is not sensitive
cannot be estimated using experimental data and may be fixedusing literature values. Finally, the global optimization algorithm
DREAM (Vrugt, 2016) was used to fit the modeling results to theexperimental data. The ability of the modeling approach to reproduce the data was evaluated using the sample standard deviation.

66Materials and Methods

We first describe the root growth model of Jones et al. (1991) and
its implementation into the HYDRUS software packages. Then,
we describe the aeroponic root growth experiment, which wasused to validate the implemented root growth and temperaturestress models. Finally, we describe the statistical approaches used
to analyze the collected experimental data, such as the sensitivity
analysis and parameter optimization.

Modeling Root Growth as a Functionof Environmental Stresses

Root Growth Model of Jones et al. (1991)
Jones et al. (1991) proposed a root growth model that simulates
daily root growth in a layered soil as a function of different environmental conditions that are characterized by growth stressfactors. Root growth and the development of root length density
depend on environmental stress factors that range from 0 (nogrowth) to 1 (no stress). The model simulates (i) daily increases in
the rooting depth, (ii) the length/weight ratio of new roots, (iii)root proliferation within different soil layers, and (iv) root senescence under the influence of various environmental factors. Joneset al. (1991) divided the stress factors into static and dynamic stressfactors. Stress factors due to Al toxicity, Ca deficiency, and coarse
fragments are considered to be static. Stress factors due to extreme
p. 2 of 16

temperatures, soil strength, and poor aeration are considered tobe dynamic.
To evaluate a potential daily increase in the rooting depth, Jones etal. (1991) proposed using either a function of the growth stage ora function of thermal time. Both functions require an additional
input that should be provided by an output of a crop growth model.
An actual increase in the rooting depth is then obtained by reducing the potential increase due to various stress factors (for moredetails, see the equations below).
The model also evaluates the depth-dependent root distribution
within the soil profile based on the root length/mass ratio, which
is used to describe a potential increase in root mass in a particularlayer and which additionally depends on the root length density.
To calculate an increase in the root mass in a soil layer, the model
requires daily values of the dry matter allocation to the root system
as an additional input, which has to be provided, in a way similar tothe daily increases in rooting depth, by a crop growth model. The
depth distribution of the root length density can also be influenced
by root senescence, which is considered to be driven not only bytime but also by low soil water contents or poor aeration (Joneset al., 1991).

Implementation into HYDRUS
The HYDRUS programs numerically solve the Richards equation
for saturated–unsaturated water flow and convection–dispersion
type equations for heat and solute transport (Šimůnek et al., 2008,2016). Extraction of water from the soil via the root system due totranspiration is expressed via a sink term in the Richards equation.
The spatial variation in root length density needs to be specified toevaluate the spatial variation in the intensity of root water uptake.
The standard code of HYDRUS-1D offers the possibility of considering a time-dependent alteration of the root length density
distribution. The time-dependent rooting depth can be either
given as an input in a tabulated form or calculated using the
Verhulst–Pearl logistic growth function. The time-dependent
rooting depth is then used together with the Hoffman and van
Genuchten (1983) function to calculate the spatial distributionof the root length density. The standard version of HYDRUS-1D
does not consider the feedback between root growth and conditions in the soil; root growth is fully defined using input
parameters (Šimůnek and Hopmans, 2009).
The standard code of HYDRUS-2D considers the root system to
be static with time and does not allow for the dynamic development of the root system. To allow for the dynamic development of
the root system and implementation of the Jones et al. (1991) root
growth model, a two-dimensional, root length density distributionfunction (Vrugt et al., 2001a, 2001b) was combined with either
the Verhulst–Pearl logistic growth function or a tabulated inputof rooting depths (Hartmann and Šimůnek, 2015), in a similar

VZJ | Advancing Critical Zone Science

way to that in HYDRUS-1D. HYDRUS-2D can consider either
a two-dimensional or radially symmetrical root length density
distribution function and requires an additional parameter todescribe the maximum horizontal or radial distance, respectively,of the root system.
The modeling approach proposed by Jones et al. (1991) was
adapted and implemented into both HYDRUS software packages. To avoid the requirement to couple the root growth modelwith a crop growth model, the modeling approach of Jones et al.
(1991) was slightly modified. In the newly developed root growthmodule, the development of the potential root system is still evaluated as in the standard HYDRUS models. The potential rooting
depth, reached when the development of the root system would be
independent of environmental conditions, is still evaluated using
a time-dependent growth function. However, an actual rooting
depth is evaluated from the potential rooting depth by taking into
account various environmental stress factors, in a similar way tothe approach of Jones et al. (1991). Similarly, an actual spatial pattern of the root length density distribution is evaluated based on
the specified shape functions. Such functions represent the potential root length density distribution under optimal developmental
conditions when no restrictions on root growth occur. Due to suboptimal environmental conditions, the potential root distributionis altered to get an actual root distribution.

Stress Factors
Stress factors Sx (dimensionless) influencing the vertical penetration of roots according to Jones et al. (1991) are listed in Table 1.
Descriptions of the static stress factors accounting for the influenceof Ca deficiency, Al toxicity, and excessive coarse fragments were
provided by Jones et al. (1991) and Hartmann and Šimůnek (2015).
Descriptions of the dynamic stress factors, which can vary withtime and are thus of greater importance for a dynamic alterationof the root system development, are given below.
According to Jones et al. (1991), the bulk density, soil texture, andwater content influence the soil strength and can thus represent
a stress for root development. The soil strength stress factor S St
proposed by Jones et al. (1991) is calculated asæp ö
SSt = S BD sin ççç q f ÷÷÷
è2 ø

[1]

Table 1. Stress factors influencing root growth.
Parameter

Stress factor (S)ST

temperature stressSCa

Ca deficiency

SAl

Al toxicity

S St

soil strengthSAeaeration stress

SCF

excessive coarse fragments

p. 3 of 16where the parameter q f denotes the fractional water content andis calculated for each depth:
qf = 0q-q L
q U -q L
qf =1
qf =

q < qLqU ³ q ³ qL

[2]q > qU

where q is the water content, q L is the lower limit of the plantextractable soil water (conventionally corresponding to the wilting
point), and q U is the drained upper limit of the soil water. In the
HYDRUS implementation of this stress factor, corresponding
pressure heads (instead of water contents) representing these limits
have to be specified. These limiting pressure heads are then used to
calculate corresponding water contents using the soil water retention curve.
The parameter S BD in Eq. [1] is a factor accounting for the interacting effects of texture and bulk density of the moist soil. The
factor is calculated as
S BD = 1

BD < BDOBD X - BD
S BD =BD X - BDO
S BD = 0

BDO ³ BD ³ BD X

[3]

BD > BD X

The parameters BD, BDX, and BDO denote the user-defined bulk
density and the bulk densities above which root growth is first
affected and completely inhibited. The parameters BDX and BDO
are calculated using the user-defined percentage (w/w) of sand
(WPS) of the considered soil:
BD X = 1.6 + 0.004 WPS
BDO = 1.1 + 0.005WPS[4]

The aeration stress factor SAe considers the effect of poor aerationin highly saturated soils. The influence of this factor is calculated as
S Ae = S FT + (1-q e )S Ae = 1

1- S FT
1-q eCq e ³ q eC

An overview of several approaches to describe the influence oftemperature on root growth rates using different mathematicalexpressions was provided by Yan and Hunt (1999) and Li et al.
(2008). In our study, two functions were considered to accountfor the crop-specific dependence of root growth on temperatureand to define the temperature stress factor. The approach of Yin
et al. (1995) considers a flexible, bell-shaped, nonlinear function:ctïìï é T -T t ù é T t -T ù (Topt -Tmin )/(Tc -Topt ) ïüï
(
)
(
)
c
min
ï
ï
úê
ú
S T ( t ,T ) = í êêýïï ê Tc -Topt úú êê Topt -Tmin úúïï [6]ûëûïî ëïþ
S T ( t ,T ) = 0
for T ( t ) < Tmin and T ( t ) > Tc

where T is the temperature at the current time level, Tmin is theminimum temperature when root growth starts, Tc and Topt are
the critical and optimum temperatures for root growth, respectively, and ct denotes the shape coefficient of the bell function.
The second function is a sinusoidal stress function, which was
described by Jones et al. (1991):

é p T ( t ) -T ù
min ú
S T ( t ,T ) = sin êê
ú
T
T
2
opt
min úûêë
S T ( t ,T ) = 0
for T ( t ) < Tmin
S T ( t ,T ) = 1[7]

for T ( t ) > Topt

where T is again temperature at the current time level, and Tmin
and Topt are the minimum and optimum temperatures for root
growth, respectively. HYDRUS users can select whether to use airtemperature or temperature in the root zone in the temperaturestress functions Eq. [6] and [7].
The difference in shapes between the two considered stress factor
functions is displayed in Fig. 1. Both approaches use genotypebased parameters such as the minimum, optimum, and criticaltemperatures for root growth. These parameters, also referred[5]q e < q eC

where q e is the water-filled porosity, which is defined as the ratioof the actual water content and the saturated water content, andq eC is the critical water-filled porosity above which the root system
is affected by poor aeration. This value can be either defined by
the user or calculated based on the approach given by Jones et al.
(1991) from soil texture.
The parameter S FT (0–1), which is crop type dependent, defines
the fraction of normal root growth when the pore space is saturated. Its value has to be defined by the user. Setting SFT equal to
1 implies that flooding has no effect on root growth. A value of 0means that there is no growth in saturated soils.

VZJ | Advancing Critical Zone ScienceFig. 1. Temperature stress functions according to Yin et al. (1995) and
Jones et al. (1991). The positions of the cardinal temperatures—theminimum (Tmin), optimum (Topt), and critical (Tc) temperatures forroot growth—influence the skewness of the graphs.

p. 4 of 16

to as cardinal temperatures, are crop specific and can be found in
the literature for several crop types (e.g., Yan and Hunt, 1999).
The equation by Yin et al. (1995) accounts for a reduction inroot growth due to high temperatures. The equation of Jones etal. (1991) considers only a reduction in root growth due to lowtemperatures and assumes that root growth is not affected by temperatures above the optimum temperature.
It is possible to enable or disable the consideration of each individual stress factor when using the newly developed HYDRUS add-on
module. Consequently, not all stress factors have to be consideredsimultaneously during simulations.

Root Growth
Two approaches are implemented to calculate the actual rooting
depth. In the first approach, the ordinary differential equation
describing population growth is used to calculate the time-dependent development of the rooting depth:

æ L ö
d Lr
= rLr çç1- r ÷÷÷çè Lm ÷ødt

[8]

where Lr [L] is the potential rooting depth when no stresses occur,
Lm [L] is the maximum possible rooting depth, t is time [T], andr [T−1] is the growth rate. The actual rooting depth can be calculated by incorporating the stress factors directly into Eq. [8]:

æ L ö
d La
= rS ( t ) La çç1- a ÷÷÷çè Lm ÷ødt

[9]where S(t) (dimensionless) represents environmental stresses and
La [L] is the actual rooting depth subject to environmental stresses.
By rearranging this equation and solving it using an explicit scheme,
La can be computed at each time step:

La ( t i ) = La ( t i-1 ) + r min éë S1 ( t i ), S2 ( t i )ùû La ( t i-1 )é L (t )ù
´ êê1- a i-1 úú ( t i - t i-1 )
Lm û
ë

[10]

where S1 and S2 (dimensionless) are stress factors influencing thevertical penetration of roots according to Jones et al. (1991).
With this approach, previous reductions in root growth are taken
directly into consideration when evaluating future root systemdevelopment. Because the growth of the root system depends onthe growth of the shoot, this approach implies that the root systemis the limiting factor for the entire crop development. Thus, rootand shoot growth are assumed to be synchronized, and previous
reductions in root development will affect shoot development,
which in turn affects future root growth.
Because a reduction in the root growth rate does not necessarily imply an equal reduction in the shoot growth rate, the second

VZJ | Advancing Critical Zone Science

approach accounts heuristically for the impact of various stresseson the reduction of root growth. This approach is based on using
a time-dependent growth function to calculate the potential
increase in the rooting depth at each time step. In this second
approach, two time-dependent growth functions are considered
to describe the potential increase in the rooting depth with time.
The first function is the Verhulst–Pearl logistic growth function,
which is the analytical solution of Eq. [8], which has already been
implemented in the standard HYDRUS-1D code and is now alsoavailable in the updated HYDRUS-2D code (Šimůnek et al., 2016):

æö
Lo
÷÷÷ LmLr ( t i ) = ççççè Lo + ( Lm - Lo )exp(-rt i ) ÷ø

[11]

where Lr [L] and Lo [L] are the potential and initial rooting depths,
respectively, Lm [L] is the maximum possible rooting depth, r [T−1]is the growth rate, and ti [T] is the time step. The second function
to describe time-dependent root growth is a sinusoidal function
proposed by Borg and Grimes (1986) based on analysis of 135 fieldobservations of 48 crop species:

é
ù ïü
æt ö
ïì
Lr ( t i ) = Lm ïí0.5 + 0.5sin êê 3.03ççç i ÷÷÷-1.47 úú ïýïïè t m ÷øêëúû ïþï
î

[12]

where Lr and Lm are again the potential and the maximum possible rooting depths, respectively, and tm [T] is the time when the
maximum rooting depth is achieved.
The difference between rooting depths under unstressed conditions at two sequential points in time denotes a potential increasein the rooting depth:d p ( t i ) = Lr ( t i )- Lr ( t i-1 )[13]

where dp [L] is the potential increase in the rooting depth, Lr [L]is the potential rooting depth at a certain point in time calculatedusing a time-dependent root growth function, and ti [T] is thetime step. Thus, dp denotes an increase in the rooting depth whenzero influence by any environmental factor such as temperature
is considered during the entire root growth period. The potential
growth function thus assumes optimal conditions not only at aparticular time but also prior to this time.
The actual increase in the rooting depth dr [L] is calculated using
the stress factors Si (dimensionless):dr ( t i ) = d p ( t i )min éë S1 ( t i ), S2 ( t i ), S0 ( t i-1 )ùû
S1 = min ( ST , SCa , S Al )0.5

S2 = min ( SSt , S Ae , SCF )
S0 =

[14]

La ( t i-1 )Lr ( t i-1 )p. 5 of 16

Integration of actual increases in rooting depths with time leads
to the actual rooting depth La.
The approach of Jones et al. (1991) assumes that the environmental
stresses S St, SAe, and SCF have a lower influence on root growth
due to possible compensation effects and that the most limiting stress factor determines the reduction in root growth. Other
approaches can be found in the literature. For example, Sommaet al. (1998) assumed that the potential increase in the rooting
depth is multiplied by all stress factors accounting for the effects
of soil strength, temperature, and nutrient concentration. Prasad
et al. (2008) also discussed the influence of combined stresses.
They pointed out that applying the minimum stress factor is acommon procedure in many models but also stated that describing
the interactions of different stress factors is still one of the greatest challenges in crop modeling. Because the combined effects ofmultiple stress factors are unknown, the approach proposed byJones et al. (1991) was selected in this study.
The parameter S0 in Eq. [14] introduces an additional active stressfactor into the calculation of the actual increase in the rooting
depth that takes into account the influence of past reductions
in growth. This factor is defined as the ratio of the actual rooting depth to the potential rooting depth at the previous time
step. As a result, past reductions in root growth affect actual root
growth. The same approach is used in HYDRUS-2D for bothvertical and horizontal extents of the rooting system. While in
HYDRUS-1D the stress factors are evaluated at the rooting depth,in the HYDRUS-2D model the stress factors reflect average conditions in the root zone.
We have implemented these two approaches because we believethat they their own advantages and disadvantages and thus canbe complementary to HYDRUS users. For example, the former(theoretical) approach ultimately leads to the actual rooting
depth reaching the maximum rooting depth as long as there
is enough time for that. Note that growth continues as long as[(1 − La(ti−1)/Lm] in Eq. [10] is positive. In this approach, future
growth can thus fully compensate for past stresses. On the other
hand, in the latter (heuristic) approach, earlier reductions in root
growth have a permanent effect on future growth.

Root Length Density Distribution
Five additional functions (Jones et al., 1991; Vrugt et al., 2001b;
Hao et al., 2005) describing the spatial distribution of the root
length density have been implemented into the root growthmodule of HYDRUS-1D (Table 2) in addition to the Hoffman
and van Genuchten (1983) root distribution function, which wasavailable previously. The functions listed in the table are expressedusing a normalized rooting depth, zr. While the functions in Table
2 are dimensionless, they are always internally normalized in the
HYDRUS models so that their integration over depth is equal to
1. In the HYDRUS models, the time-dependent actual rooting

VZJ | Advancing Critical Zone Science

Table 2. Root length density distribution b(zr) models implemented in
HYDRUS-1D.
Model b(zr)†

Reference

1

2.938 − 2.462zr

modified after Zuo et al. (2004)

2

4.522(1 – zr)5.228exp[(9.644zr)2.426

3

2.21 – 3.72zr + 3.46(zr)2 – 1.87(zr)2 Wu et al. (1999)

4(1 – zr)exp[−(Pz/zm)|z* − zr|]

5

1 – zr/3)3

Vrugt et al. (2001a)
Jones et al. (1991)† zr, normalized rooting depth; Pz and z*, empirical parameters.

depth is an important input factor for these functions to calculate the root length density distribution. Of these functions, theroot length density distribution model proposed by Vrugt et al.
(2001b) is the most flexible model. Root length density distributions given by these functions represent the potential root lengthdensity distribution pattern, which would be obtained under optimal conditions when no restrictions of root growth, such as rootsenescence, were considered.
The actual root length density distribution is calculated similarly
to the actual rooting depth. The actual increase in root lengthdensity (RLDga) at each time step is derived from the potential
increase (RLDgp), which is defined as the difference between theroot length density distribution [b(zr)] at the current (ti) and previous (ti – 1) time steps:

RLD gp ( t i ) = b( zr , t i )- b( zr , t i-1 )[15]

An increase in the potential root length density in each node
can be reduced due to root senescence. This approach is usedin HYDRUS-1D as well as in HYDRUS-2D. Similarly to theoriginal version of the approach of Jones et al. (1991), the adapted
version implemented into HYDRUS considers root senescenceduring the life of the plants. Only a few minor changes had to bemade in the modeling approach of Jones et al. (1991) to make itcompatible with HYDRUS. In HYDRUS, the root senescence(Sen) is calculated as

Sen ( t i ) = rd RLD gp ( t i ){1 + max éë1-q f ( t i ),1- S AE ( t i )ùû } [16]where rd is a scaling factor (0–1) defining the influence of rootsenescence. This approach assumes that root senescence can beincreased due to low water availability or poor aeration. Moreinformation regarding the evaluation of root senescence was provided by Hartmann and Šimůnek (2015).
The actual root length density distribution ba(ti) is then calculated asba ( t i ) = ba ( t i-1 ) + RLD ga ( t i )RLD ga ( t i ) = RLD gp ( t i )- Sen( t i )[17]

p. 6 of 16

The implemented root growth module was used to simulate theimpact of various environmental factors on root growth using two
hypothetical examples.

Experiment Description and Measurements
Jones et al. (1991) stated that their approach is only the first approximation of the effects of environmental factors on root growth thatneeds to be verified against experimental data. Experimental datafrom an aeroponic system (described below) were selected to evaluate the described approach to model the influence of temperatureon the vertical penetration of the root system. The experimental
setup in an aeroponic system with no soil and no water movement
has the advantage of evaluating the effects of temperature on root
growth while limiting the effects of other factors.
The experiment was conducted in a greenhouse with bell pepper
at the Jacob Blaustein Institutes for Desert Research, in Midreshet
Ben-Gurion, Israel (30°51¢5² N, 34°47¢0² E, altitude 480 m). The
objective was to evaluate the effect of three different root zone
temperatures on root and plant growth. Six bell pepper plants were
cultivated from 7 Jan. to 20 Feb. 2015 in aeroponic pots mounted
on top of the aeroponic systems. The aeroponic apparatus comprised circular pots made from plastic material with a diameter of
50 cm and a depth of 14 cm. Within each thermally isolated aeroponic pot, misters (Coolnet, Netafim Israel) were fixed to produce
the desired fine mist sprayed directly onto the plant roots. The
computer-controlled spraying varied from 8-s sprayings at 1-min
intervals, depending on the growth stage (size) of the plants andthe temperature of the greenhouse. Three different water temperatures were applied to each treatment, namely 7, 17, and 27°C. The
air and root zone temperatures were measured daily. The air temperatures in the greenhouse were 25°C during the day and 18°C atnight. The treatments were replicated twice, leading to six tanks intotal. The maximum rooting depth of each plant was observed fourtimes during the 44 d of plant growth. In the aeroponic systems,
all other factors that affect the growth and development of roots
were rendered insignificant during the treatment.

Validation of the Temperature-Dependent
Root Growth Modeling Approach
Evaluated Models

The implemented modeling approach to simulate the temperaturedependent vertical root penetration was evaluated by comparing
modeling results with the measured maximum rooting depths
in the experimental aeroponic system with bell pepper. Because
two approaches were implemented in HYDRUS to describe
both time-dependent potential root growth and the temperaturestress factor, four combinations are thus available to describe thetemperature-dependent root growth. All four models were tested
against measured maximum rooting depths to validate their ability to describe the temperature- and time-dependent vertical root
penetration. The four combinations are summarized in Table 3.
Were the these root growth models (see Table 3) able to properly

VZJ | Advancing Critical Zone Science

Table 3. Overview of four models (combinations of time-dependent
potential root growth and the temperature stress factor equations) usedin this study to evaluate the temperature-dependent root growth.
Model

Potential root growth

Temperature stress factor

A

Verhulst–Pearl logistic growth

Jones et al. (1991)

B

Verhulst–Pearl logistic growth

Yin et al. (1995)

C

Borg and Grimes (1986)

Jones et al. (1991)

D

Borg and Grimes (1986)

Yin et al. (1995)

simulate the influence of temperature on vertical root penetration, a single combination of model parameters for each modelthat could reproduce the measured maximum rooting depths for
all temperature treatments would have to exist.
An overview of the model parameters that have to be specified foreach model is given in Table 4. Table 4 shows that, depending on
the model, four or six parameters have to be specified to model
temperature-dependent root growth. In a complex soil water flow
model such as HYDRUS, four to six additional parameters can significantly increase the calibration effort and parameter uncertainties.
The temperature-dependent modeling approach was tested outsideof the HYDRUS implementation and within a MATLAB environment. The goal of the evaluation was to determine whether the
combination of the time-dependent root growth functions and thetemperature stress functions were able to reproduce the measured
rooting depths under the given boundary conditions.

Sensitivity Analysis
A global sensitivity analysis was conducted using the Sobol¢
method to reveal the key parameters of each model and to determine the contribution of the uncertainty of each parameter to theuncertainty of the model output. The Sobol¢ method is based on
variance decomposition and provides the impact of each parameterand its interactions with other parameters on the model output(Sobol¢, 1993). This type of global sensitivity analysis can beapplied to nonlinear and non-monotonic models and is a widelyused tool for sensitivity analysis studies. Its ability to account forinteractions between model parameters is an important advantageof the Sobol¢ method (Rosolem et al., 2012).
Table 4. Parameters of the four models for calculating time- and temperature-dependent root growth.
ModelParameters†A

Lm, r, Tmin, Topt

B

Lm, r, Tmin, Topt, Tc, ct

C

Lm, Tmin, Topt, tm

D

Lm, Tmin, Topt, tm, Tc, ct† Lm, maximum possible rooting depth; r, growth rate; Tmin, minimum
temperature when root growth starts; Topt, optimum temperature for root
growth; Tc, critical temperature for root growth; ct, shape of the bell function; tm, time when the maximum rooting depth is achieved.

p. 7 of 16

An overview of studies using the Sobol¢ method for sensitivity analysis in hydrological modeling was provided by Song et al. (2015).
The method has already been applied with the HYDRUS software
package by Li et al. (2012), Brunetti et al. (2016) (HYDRUS-1D),and Wang et al. (2016) (HYDRUS-2D).
Sobol¢ (1993) proposed that the total variance of the model
output can be decomposed into component variances of individual parameters and their interactions. The first-order sensitivity
index quantifies the main effect of the ith parameter, Xi . Thissensitivity index denotes the part of the total variance due to Xi
without considering the interactions with other parameters. The
total-order sensitivity index additionally includes the proportion of
the variance due to the interactions of Xi with the other parameters.
The values of the indices vary from 0 to 1, where 0 stands for no
influence and 1 for a high influence on the variance.
The first-order (Si) and total-order sensitivity (Sti) indices wereapproximated by a numerical Monte Carlo estimation proposedby Saltelli et al. (2010):
éù
(1 N )å f ( B ) j êê f ( A iB ) j f ( A ) j úú
ë
û
Si =
2
N é
2
ù
(1 N )å j=1 ëê f ( A ) j ûú - f oN
j =1

S ti =éëù2
û

(1 2 N )å Nj=1 êê f ( A ) j - f ( A iB ) j úú
(1 N )å

N é
fj =1 êë

2( A ) j ùú - f o2
û

[18]

[19]

where A and B denote two matrices of data sets generated using
Sobol¢ quasi-random sequences, which consist of a sample ofnumbers between 0 and 1 distributed in a p-dimensional unit
hypercube. In Eq. [18], A Bi denotes a matrix where all columns
are from Matrix A except for the ith column, which is from Matrix
B, and fo is calculated as (1/N) S Nj=1 f (A)j. Each matrix has thedimension N ´ p, whereby p is the number of parameters and N
denotes the number of required simulation runs. Each matrix thus
consists of N sets of parameters.
The parameters were transformed into their defined parameter
space based on the specified minimum and maximum values of theparameters (listed in Table 5). Parameter limits for rooting depths
under favorable environmental conditions were taken from Borg
and Grimes (1986). The limiting values of cardinal temperatures
were specified based on information provided by Saha et al. (2010).
The limits of parameter tm were chosen in a way that common
durations of vegetation periods would lie within these boundaries.
Parameters f (A) and f (B), as well as f (A Bi), represent the model
output in the form of a chosen statistical metric. In this study,
the root mean square error (RMSE) was chosen for the statistical metric:

VZJ | Advancing Critical Zone Science

Table 5. Upper and lower boundaries of the parameter space for theparameters of the four models considered in the sensitivity analysis andoptimization.
Parameter†

Lower boundary

Upper boundary

Value forsensitivity analysis

Lm, cm

70

150

110

r, d−1

0.1

1

0.15

Tmin, °C

019

8Topt, °C

20

32

23

Tc, °C

33

50

35

ct

0.01

100

0.15

tm, d

43

365

150

† Lm, maximum possible rooting depth; r, growth rate; Tmin, minimum
temperature when root growth starts; Topt, optimum temperature for root
growth; Tc, critical temperature for root growth; ct, shape of the bell function; tm, time when the maximum rooting depth is achieved.

2

RMSE =

å( LObs - LSim )n[20]which describes the square root of the average squared differences
between the observed (LObs) and simulated (LSim) rooting depths.
The RMSE measure was selected to evaluate prediction errors
because only one output of the root growth models, i.e., the rooting depth, was analyzed in the sensitivity analysis.
The number of parameter sets (N) in the sensitivity analysis of thefour models was set to 10,000. This number was initially set higher
than in studies of Brunetti et al. (2016) and Zhang et al. (2013)
to avoid a time-consuming convergence analysis of the sensitivity
analysis and to achieve a higher accuracy of the sensitivity analysis,
which increases with an increasing number of model runs. The p
parameter depends on the considered model and is either four (forModels A and C) or six (for Models B and D). To calculate thesensitivity indices for all i = 1, …, p parameters, Matrix A Bi has tobe evaluated p times. The total number of model runs required to
calculate the sensitivity indices for all parameters of each model
were M = N(2 + p). Archer et al. (1997) suggested using bootstrap
confidence intervals (Efron and Tibshirani, 1993) to evaluate asuitable accuracy of the sensitivity estimates. For this reason, eachestimation of the sensitivity indices was repeated 500 times toevaluate the 25th and 75th percentiles of the sensitivity indices.
The small number of repetitions is due the fact that the sensitivity analysis of four models, each including four or six parameters,
requires a high level of computational effort. Evaluated percentiles
were used only as an additional parameter to assess the sensitivity
indices (Archer et al., 1997).
Rather than using experimental data, which may be subject to various errors and effects of various factors, the sensitivity analysis was
performed using a hypothetical data set in which observed rooting
p. 8 of 16

depths LObs were generated by running all four models (Table 3)with a predefined parameterization. The cardinal temperatures
Tmin, Topt , and Tc were set to 8, 23, and 27°C, respectively. A
complete list of parameter values for the model runs is provided
in Table 5. The boundary conditions were similar to those of the
aeroponic experiments. The temperature was set to 22°C during
the first 14 d and to either 7, 17, or 27°C during the remaining time
period of 150 d. Additionally, a fourth scenario was consideredwith a temperature of 37°C after the first 2 wk to make sure thatthe specified cardinal temperatures would lie within the appliedboundary conditions to determine their influence on the model
output. Each model was thus executed four times with four different temperature boundary conditions to generate data for Lsim
in Eq. [20].

Parameter Optimization
The Differential Evolution Adaptive Metropolis (DREAM) algorithm (Vrugt et al., 2009; Vrugt, 2016) was used for optimizing
the model parameters and for model calibration. The DREAM
algorithm is based on Bayesian statistics; it runs multiple different Markov chains to generate a random walk through the search
space. Based on a proposal distribution, the sampler evolves tothe posterior distribution by iteratively finding solutions withstable frequencies stemming from the fixed probability distribution (Laloy and Vrugt, 2012). The Gaussian likelihood function
was used to summarize the distance between model simulations
and corresponding observations. The residuals were assumed to
be independent (uncorrelated) and normally distributed while
the measurement error was neglected. The latest MATLAB
implementation of DREAM (Vrugt, 2016) was used for model
parameter optimization.
Eight Markov chains were run with a set of 5000 generations. The
initial state of each chain was sampled from a Latin hypercube.
The parameter space of each parameter was defined by using thesame boundaries as were used for the sensitivity analysis (see Table5). These parameter limits also define the search domain for the
predominantly physically based parameters (Vrugt, 2016). The
calculation time for the DREAM optimization of a single model
was approximately 10 min, with no parallelization needed on amachine with the following specifications: Intel Core i7–4710HQ
CPU with 2.50 GHz of RAM and 12 GB of storage.

66Results and Discussion

We evaluating the effects of different factors on root growth withtwo examples using the new root growth module in HYDRUS-1Dand HYDRUS-2D. Then we analyzed the collected experimental
data and the evaluation of the temperature-dependent root growthmodeling approach. Therefore, we first collected data on the experimental outcomes. Second, we carried out the sensitivity analysis
to evaluate the sensitivity of the modeling results to various input
parameters and identified which parameters need to be fitted andwhich can be set to values from the literature. Third, we used theDREAM optimization approach to analyze the collected experimental data while considering the results of the sensitivity analysis.

Applications of the Root Growth Module
We used two hypothetical examples that demonstrate theimplemented root growth model and the impact of various
environmental factors on root growth. In the first example, weused HYDRUS-1D and simulated optimal root growth as wellas root growth restricted due to low water availability, temperature, texture, and bulk density. In the second example, we used
HYDRUS-2D to again simulate optimal root growth and thenroot growth affected by a nonuniform distribution of water contents due to asymmetrical irrigation.

HYDRUS-1D Examples
Figure 2 shows examples of the development of simulated rootsystems under the influence of various environmental factors(Scenarios 1–4) compared with the potential development of theFig. 2. Simulated potential (0) and actual (1–4) root length density distributions at different times (different shades of gray): (1) under drought
conditions, (2) under the influence of texture and bulk density, (3a and 3b) under the influence of different root zone temperatures, and (4) in a heterogeneous soil system.

VZJ | Advancing Critical Zone Sciencep. 9 of 16

root system, independent of environmental factors (Scenario 0).
The soil profile was considered to be homogenous, consisting of
10% sand, 50% silt, and 40% clay and having the bulk density of1530 kg m−3. The parameters for the soil hydraulic functions of van
Genuchten (1980) and Mualem (1976) were estimated from the
textural information using the Rosetta module of HYDRUS-1D.
Time-dependent root growth was simulated using the functionof Borg and Grimes (1986). Model 2 from Table 2 was selected todescribe the potential root length density distribution. The maximum potential rooting depth (Lm) was set to 120 cm, and root
growth was considered for 90 d (tm). The upper and lower boundary conditions were set to atmospheric boundary conditions withsurface runoff and free drainage, respectively.
In Scenario 1, the development of the root length densities was
simulated for drought conditions, which may cause aeration stressand root senescence (both options were enabled in this simulation).
Potential evaporation and transpiration rates were set to increase
stepwise at the beginning (up to 28 d) and remain constant at theend of the simulation. Irrigation took place at irregular intervals(only every 3–6 d) and covered only 25% of the potential evapotranspiration. Due to low water availability, the root system wasunderdeveloped compared with the reference simulation when norestrictions on root growth were considered (Scenario 0). The root
length density in this scenario is reduced compared with Scenario
0, mainly close to the soil surface because the drought stress ishighest in soil horizons with the highest root water uptake.
Scenario 2 considered the effects of texture and bulk density (theoption of the soil strength stress factor was additionally enabled)on the development of the root length densities and verticalroot penetration. The results show that due to the increased
soil strength (a high clay content and a high bulk density), thepotential maximum rooting depth of 120 cm was not reached.
The additional stress due to the soil strength negatively influenced
the vertical penetration of roots, in addition to the effects of rootsenescence due to the drought stress in Scenario 1.
Scenario 3 considered the influence of different root zone temperatures on the root length density development. In this case, only the
option to consider the influence of temperature on root growth wasenabled. The development of the root system was simulated using
the root growth function of Borg and Grimes (1986) in combination with the temperature stress function of Jones et al. (1991).
The parameters Tmin, Topt, and tm were set to 5°C, 35°C, and 90 d,
respectively. The temperature throughout a soil depth of 150 cm was
set to a constant value of 10°C. The temperature at the soil surface
was set to 25°C during the first 14 d and to 30 and 7°C in Scenarios
3a and 3b, respectively, until the end of the simulation. The results
show that the root system in Scenario 3b, which was exposed totemperatures close to Tmin, was underdeveloped compared with theroot system in Scenario 3a, which was exposed to temperatures closeto Topt. The root system under conditions that were close to the

VZJ | Advancing Critical Zone Science

optimum temperature showed greater growth rates and thus reacheddeeper depths and developed greater root length densities than theroot system under non-optimum temperature conditions.
Finally, an additional soil horizon, consisting of 24% sand, 64%
silt, and 12% clay and having the bulk density of 1290 kg m−3,
was specified between 20 and 50 cm in Scenario 4. The samestress factors as in Scenario 2 were enabled. The results for this
scenario show that, due to the different properties of the additionalsoil horizon, which affect water availability and O2 supply, thedevelopment of the root length density distribution can be highly
influenced. The results show that the root length density can bespatially altered due to a spatial distribution of stress factors. In
this case, unfavorable growing conditions led to a reduction of root
length density in the middle layer.

HYDRUS-2D Examples
Figure 3 shows the simulated potential and actual two-dimensional
development of the root system in a homogeneous soil consistingof 50% sand, 20% silt, and 30% clay and having a bulk density of1510 kg m−3. The parameters for the soil hydraulic functions of
van Genuchten (1980) and Mualem (1976) were estimated using
the Rosetta module of HYDRUS-2D. The lower boundary condition was set to free drainage. The simulation time was set to 60 d,and root growth was considered until the end of the simulation.
Potential transpiration was set to 0.1 mm d−1 between the firstday and Day 11. From Days 12 to 29, potential transpiration was
set to 1 mm d−1; from Day 30 until the end of the simulation, it
was set to 2 mm d−1. Irrigation was applied at a rate of 2 mm d−1
between Days 30 and 50 only in the top left corner (10 cm), while
the rest of the soil surface remained unirrigated. Potential evaporation was specified only for the unirrigated soil surface, with arate of 0.1 mm d−1 from the beginning until Day 15, 1 mm d−1
between Days 16 and 18, and 2 mm d−1 from Day 19 until theend of the simulation. In addition to water availability and aeration, the texture and bulk density were also considered influenceson root growth. Only one plant in the middle of the domain was
specified with a potential rooting depth in the vertical directionof 100 cm and 30 cm in the horizontal direction. The root lengthdensity function of Vrugt et al. (2001a, 2001b) predicts an axiallysymmetrical shape of the potential root length density distribution.
Asymmetrical irrigation results in an unequal development of
the root system. The root system is more developed on the lefthand side of the transport domain than the right-hand side due togreater water availability. Additionally, the root system simulated
under the influence of growth restrictions showed a smaller rooting depth than the reference simulation when no restrictions onroot growth were considered.

Experimental Outcomes
The observed maximum rooting depths significantly differedamong the three treatments. After 44 d, the average rooting depths
p. 10 of 16

Fig. 3. Potential and actual two-dimensional root system development at different times when irrigation was applied in the top left corner.

for the treatments with root zone temperatures of 7, 17, and 27°C
were 19, 47, and 72 cm, respectively. The biomass above groundshowed the same trend.

Sensitivity Analysis
The results of the global sensitivity analysis in the form of the totaland first-order sensitivity indices of all parameters of the four root
growth models are shown in Fig. 4. The light gray bars represent
the contribution of each individual parameter to the variance ofthe model output (the first-order sensitivity). The dark gray bars
represent the individual parameter contributions, as well as thecontributions of their interactions with the other parameters to
the variance of the model output.
The results of the sensitivity analysis (Fig. 4) reveal that the differences among the influences of different parameters in Models
A and C are more pronounced than between different parameters in Models B and D. The differences between the sensitivity
indices of different parameters of Models B and D are not very
different from each other. Figure 4 reveals that the plant-specificcardinal temperatures Topt, Tmin, and Tc and the plant-specificparameter Lm (the maximum rooting depth reached when no

VZJ | Advancing Critical Zone Science

environmental stresses are present) have only a minimal impact
on the model outcome. Because these parameters are biologically
based and specific for a particular crop cultivar, their values canbe found in the literature and do not have to be calibrated. For
example, Yan and Hunt (1999) provided cardinal temperatures
for sorghum [Sorghum bicolor (L.) Moench], wheat, barley, bean
(Phaseolus vulgaris L.), and maize, and Li et al. (2008) provided
a list of sources for cardinal temperatures for winter wheat. Borg
and Grimes (1986) provided an overview of Lm for several crops.
Specifying parameters that do not exert a large influence on themodel output using values found in the literature reduces thenumber of calibrated parameters and thus the complexity andeffort of the calibration task. By defining a limit of 0.3 for the total
sensitivity index to neglect the parameter influence, the number ofcalibrated parameters of Model A can be reduced from four to two.
Based on the results in Fig. 4, the number of calibrated parametersof Model B can be reduced from six to one. However, because no
information about parameters r and ct could be found in the literature, these parameters need to be considered within the model
calibration. The number of calibrated parameters could thus only
be reduced to two for Model B. Yan and Hunt (1999) suggested
p. 11 of 16

Fig. 4. Total- and first-order indices for all parameters of the four root growth models. The bars represent the mean value of all evaluations, and theboxplots indicate the value of the median and the 25th and 75th percentiles; ct is the shape coefficient of the bell function, Lm is the maximum possible
rooting depth, r is the growth rate, tm is the time when the maximum rooting depth is achieved, Tmin is the minimum temperature for root growth, andTopt is the optimum temperature for root growth.

that reasonable results can be found by setting ct to 1. They alsosuggested that it may be possible to set Tmin to 0°C, which could
be acceptable for all crops except for summer crops.
The results of the global sensitivity analysis for Models C and D(when the sinusoidal growth function of Borg and Grimes [1986] isused) show that the parameter Lm is less important than in Models
A and B. In these models, the parameter defining the time whenplants reach maturity, tm, demonstrates a high influence on themodel output. Based on the results of the sensitivity analysis, thenumber of calibrated parameters can be reduced from four to asingle parameter for Model C and from six to three for Model D,
while the other parameters can be specified using values from theliterature. In contrast to the other models, the results presented inFig. 4 suggest the parameter Tc should be considered as a calibration parameter for Model D. Due to the lack of knowledge about
the value of ct, this parameter should still be considered in themodel calibration.

Parameter Optimization
The optimized parameter values of all four modeling approaches
are listed in Table 6. In addition to these, the table also provides an overview of the goodness of fit obtained by individualmodels determined by the RMSE. This value indicates how well

VZJ | Advancing Critical Zone Sciencethe rooting depths simulated by the models using the optimized
parameter sets fit the observed maximum rooting depths.
Figure 5 visualizes the maximum rooting depths measured for all
temperature treatments in the aeroponic system and simulatedusing the four models and the optimized parameter sets in Table
6. The upper left figure shows the simulation results for all fourmodels compared with the average observed maximum rooting
depths (averaged for six plants) in each tank. These averaged values
were used as data in the comparison within the optimization task.
Table 6. Parameter values for all four models optimized using theDREAM optimization algorithm, including the maximum possible
rooting depth (Lm), the growth rate (r), the minimum temperaturewhen root growth starts (Tmin), the optimum temperature for root
growth (Topt), the critical temperature for root growth (Tc), the shapeof the bell function (ct), and the time when the maximum rooting
depth is achieved (tm). The RMSE shows the goodness of fit of themodel output to the measured data using the optimized parameter sets.
Model Lm
cm

r

Tmin

d−1

Topt

Tc

————— °C —————

ct

tm

RMSE

d

cm

A

70

0.3

1.9

31.8

–

–

–

7.6

B

71.6

0.31

2.9

30.7

37

0.22

–

7.9

C

76.5

–

4.6

32

–

–

50.3

3.1

D

144.6

–

9.3

22.5

45.5

11.8

65

2.1

p. 12 of 16

Fig. 5. Comparison of the rooting depths simulated using the parameter sets from Table 6 (all parameters were optimized) with those measured in the
aeroponic system. Data (simulated and average observed) displayed in the upper left figure (for all three temperatures ) were used to calculate the goodness of fit. The upper right, lower left, and lower right figures display the measured (for all six plants) and simulated rooting depths for temperatures Tof 7, 17, and 27°C, respectively.

The other three subplots show a comparison between the simulated rooting depths and the measured maximum rooting depthsof all six plants for all three temperature treatments.
Figure 5 and the RMSE in Table 6 indicate that Models C and D
fitted the observed data very well. Models A and B did not seem tocapture the root growth dynamics as well as the other two models.
This is underscored by the RMSE values listed in Table 6. The
graphs produced by Models A and B are quite similar, and bothresulted in an RMSE of 7 to 8 cm. In contrast, Models C and
D produced RMSE values below 4 cm. The best correspondence
between the model and the data, with the smallest RMSE of only
2 cm, was achieved with Model D, which uses the growth functionof Borg and Grimes (1986) and the stress reduction function of Yin
et al. (1995). However, despite the slightly worse performance ofModel C, the results show that both temperature stress approaches
can be used to simulate temperature-dependent root growth whenused in combination with the time-dependent root growth function of Borg and Grimes (1986).
Table 7 lists the optimization results when the numbers of calibration parameters were reduced based on the results of the sensitivity
analysis. In this case, the parameters with a total-order sensitivity

VZJ | Advancing Critical Zone Scienceindex lower than 0.3 were set equal to values from the literature.
The parameters describing cardinal temperatures were taken from
Reddy and Kakani (2007), who investigated cardinal temperaturesof bell pepper for pollen germination. Due to the low sensitivityof cardinal temperatures on the model output, these temperatures
could also be used to model the root growth, even though thishas not yet been scientifically documented. The value for Lm wastaken from Borg and Grimes (1986). Based on the results of theTable 7. Optimized parameter sets obtained when the insensitive
parameters were specified with values taken from the literature, including the maximum possible rooting depth (Lm), the growth rate (r), theminimum temperature when root growth starts (Tmin), the optimum
temperature for root growth (Topt), the critical temperature for root
growth (Tc), the shape of the bell function (ct), and the time when the
maximum rooting depth is achieved (tm). Italic values indicate insensitive parameters that were not considered within the optimization
process.
Model Lm
cm

r

Tmin

d−1

Topt

Tc

————— °C —————

ct

tm

RMSE

d

cm

A

70

0.49

12

30

–

–

–

12

B

90

0.6

12

30

40

1.7

–

12.3

C

90

–

12

30

–

–

52

6.7

D

90

–

12

30

49

0.7

63

5.4

p. 13 of 16

sensitivity analysis, the parameter Tmin was still considered as acalibration parameter in Models A and B.
The optimization results shown in Table 7 indicate that despite thereduction in the number of calibrated parameters, the evaluatedroot growth models, especially Models C and D, produced an
acceptable correspondence with the measured data and associatedgoodness of fit. The other two models (A and B) show a larger
increase in the RMSE, which indicates a drop in the goodness offit. Simulated rooting depths for all temperature treatments using
the four root growth models and the optimized parameter setsfrom Table 7 are shown in Fig. 6.
Figure 6 and the RMSE in Table 7 indicate that Model D also
fit the observed data very well, even when three of the six model
parameters were specified using values taken from the literatureand not calibrated. Model C, which only had one calibration
parameter, provided a sufficient reproduction of the observed
maximum rooting depths.

66Summary and Conclusions

We implemented an adapted root growth modeling approachbased on Jones et al. (1991) in the HYDRUS software packages.
The modeling approach describes root growth and root systemdevelopment under the influence of environmental factors such as
soil hydraulic properties, soil strength, water availability, chemical
conditions, and temperature. The implementation of growth andstress functions in the HYDRUS software opens the opportunity
to derive parameters of these functions from laboratory or field
experimental data using inverse modeling. This option was thendemonstrated using experimental data from an aeroponic system
for the temperature stress part of the root growth model.
The impact of temperature on root growth was describedusing two stress response functions, both combined with two,time-dependent root growth functions. The four resulting combinations were tested against rooting depth data measured using
an aeroponic system experiment. Because the experiment was notinitially set up to evaluate the described root growth module, thecollected data were not fully suitable to identify model parameterswith high accuracy. For example, the use of a wider temperature
range would help to better identify cardinal temperatures and the
maximum potential rooting depth (Lm). Nevertheless, the results
indicate that the temperature-dependent root growth approaches
are well capable of reproducing real root growth data. All four of
the tested root growth models reproduced the measured rooting
depth data very well, with RMSE values ranging between 2.1 and8 cm. Especially the results of the root growth model referred toas Model D, which combined the Borg and Grimes (1986) root

Fig. 6. Comparison of the rooting depths simulated using the parameter sets from Table 7 (only sensitive parameters were optimized) with thosemeasured in the aeroponic system. Data (simulated and average observed) displayed in the upper left figure (for all three temperatures) were used to
calculate the goodness of fit. The upper right, lower left, and lower right figures display measured (for all six plants) and simulated rooting depths fortemperatures T of 7, 17, and 27°C, respectively.

VZJ | Advancing Critical Zone Sciencep. 14 of 16

growth function with the Yin et al. (1995) temperature stress function, showed an excellent agreement with the observed rooting
depth data. The results indicate that the implementation of thetime-dependent root growth function of Borg and Grimes (1986)
into the HYDRUS software was beneficial because this function captured the actual root growth better than the originallyimplemented Verhulst–Pearl logistic growth function. However, acomplete evaluation of the two proposed (theoretical and heuristic)
approaches simulating root growth still requires further testingunder dynamic stress conditions that would occur under field conditions, e.g., by using data provided by Aidoo et al. (2018).
A sensitivity analysis revealed that the biologically based parameters of the models, such as the cardinal temperatures and the
maximum rooting depth, tend to have a small impact on the modeloutcome. Because of the low sensitivity of the model outputs tothe uncertainties originating from these parameters, it seems
reasonable to approximate these parameters with values from theliterature and to reduce the number of parameters that need to becalibrated. We have shown that when values from the literature
are used for the model-insensitive parameters, the number of calibrated parameters can be reduced for Model A from four to two,for Model C from four to one, for Model B from six to four, andfor Model D from six to three. The reduced number of calibrated
parameters leads to a reduced complexity of the calibration task
and a reduced uncertainty in optimized parameters, especially inconnection with the HYDRUS programs. The small number of
additional calibration parameters makes the modeling approachof temperature-dependent root growth a convenient add-on to the
HYDRUS models. However, the implementation of the influenceof other stress factors (other than temperature) on the vertical root
penetration, as well as on the root length density distribution, into
the HYDRUS models still needs to be validated against experimental data.

Acknowledgments

This research was partially supported by the Bayer Fellowship Program in the form of the Jeff
Schell Scholarship, the Israel Ministry of Agriculture and Rural Development (Eugene Kandel
Knowledge Centers) as part of the “Root of the Matter: The root zone knowledge center forleveraging modern agriculture,” and the BonaRes project (BOMA03037514) of the Federal
Ministry of Education and Research (BMBF), Germany.

References

Borg, H., and D.W. Grimes. 1986. Depth development of roots with time: An
empirical description. Trans. ASAE 29:194–197. doi:10.13031/2013.30125
Brouwer, R., and C.T. de Wit. 1968. A simulation model of plant growthwith special attention to root growth and its consequences In: W.J.
Whittington, editor, Proceedings of the 15th Easter School of Agricultural Science, Nottingham, UK. Butterworths, London. p. 224–242.
Brunetti, G., J. Šimůnek, and P. Piro. 2016. Sensitivity analysis and global optimization of soil hydraulic properties of a permeable pavement at thefield scale. J. Hydrol. 540:1146–1161. doi:10.1016/j.jhydrol.2016.07.030
Clausnitzer, V., and J.W. Hopmans. 1994. Simultaneous modeling of transient three-dimensional root growth and soil water flow. Plant Soil
164:299–314. doi:10.1007/BF00010082
Couvreur, V., J. Vanderborght, and M. Javaux. 2012. A simple three-dimensional macroscopic root water uptake model based on the hydraulic architecture approach. Hydrol. Earth Syst. Sci. 16:2957–2971.
doi:10.5194/hess-16-2957-2012
Dexter, A.R. 1986. Model experiments on the behaviour of roots at theinterface between a tilled seed-bed and a compacted sub-soil. Plant
Soil 95:123–133. doi:10.1007/BF02378858
Efron, B., and R.J. Tibshirani. 1993. An introduction to the bootstrap. Chapman and Hall, New York. doi:10.1007/978-1-4899-4541-9
Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulation of field water
use and crop yield. John Wiley & Sons, New York.
Fowler, D.B., L.V. Gusta, and N.J. Tyler. 1981. Selection for winter hardiness in wheat: III. Screening methods. Crop Sci. 21:896–901.
doi:10.2135/cropsci1981.0011183X002100060023x
Garnett, T., V. Conn, and B.N. Kaiser. 2009. Root based approaches to improving nitrogen use efficiency in plants. Plant Cell Environ. 32:1272–
1283. doi:10.1111/j.1365-3040.2009.02011.x
Grant, R.F. 1993. Simulation of soil compaction and root growth. Plant Soil
150:1–14. doi:10.1007/BF00779170
Hao, X., R. Zhang, and A. Kravchenko. 2005. Effects of root density distribution models on root water uptake and water flow under irrigation.
Soil Sci. 170:167–174. doi:10.1097/00010694-200503000-00002
Hartmann, A., and J. Šimůnek. 2015. HYDRUS root growth module, Version 1. HYDRUS Softw. Ser. 5. Dep. of Environmental Sciences, Univ. of
California, Riverside.
Hoffman, G.J., and M.Th. van Genuchten. 1983. Soil properties andefficient water use: Water management for salinity control. In:
H.M. Taylor et al., editors, Limitations to efficient water use incrop production. ASA, CSSA, and SSSA, Madison, WI. p. 73–85.
doi:10.2134/1983.limitationstoefficientwateruse.c5
Jones, C., W.L. Bland, J.T. Ritchie, and J.R. Williams. 1991. Simulation ofroot growth. In: J. Hanks and J.T. Ritchie, editors, Modeling plant andsoil systems. Agron. Monogr. 31. ASA, CSSA, and SSSA, Madison, WI. p.
91–123. doi:10.2134/agronmonogr31.c6
Kaspar, T.C., and W.L. Bland. 1992. Soil temperature and root growth. Soil
Sci. 154:290–299. doi:10.1097/00010694-199210000-00005
Laloy, E., and J.A. Vrugt. 2012. High-dimensional posterior exploration ofhydrologic models using multiple-try DREAM(ZS) and high-performancecomputing. Water Resour. Res. 48:W01526. doi:10.1029/2011WR010608
Leitner, D., S. Klepsch, G. Bodner, and A. Schnepf. 2010. A dynamic rootsystem growth model based on L-Systems. Plant Soil 332:177–192.
doi:10.1007/s11104-010-0284-7
Li, L., G.S. McMaster, Q. Yu, and J. Du. 2008. Simulating winter wheat
development response to temperature: Modifying Malo’s exponential sine equation. Comput. Electron. Agric. 63:274–281.
doi:10.1016/j.compag.2008.03.006

Aidoo, M.K., E. Bdolach, A. Fait, N. Lazarovitch, and S. Rachmilevitch.
2016. Tolerance to high soil temperature in foxtail millet (Setaria italicaL.) is related to shoot and root growth and metabolism. Plant Physiol.
Biochem. 106:73–81. doi:10.1016/j.plaphy.2016.04.038

Li, Y., W. Kinzelbach, J. Zhou, G.D. Cheng, and X. Li. 2012. Modelling irrigated maize with a combination of coupled-model simulation anduncertainty analysis, in the northwest of China. Hydrol. Earth Syst. Sci.
16:1465–1480. doi:10.5194/hess-16-1465-2012

Aidoo, M.K., T. Sherman, J.E. Ephrath, A. Fait, S. Rachmilevitch, and N.
Lazarovitch. 2018. Grafting as a method to increase the tolerance
response of bell pepper to extreme temperatures. Vadose Zone J.
17:170006. doi:10.2136/vzj2017.01.0006

Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
doi:10.1029/WR012i003p00513

Archer, G.E.B., A. Sltelli, and I.M. Sobol¢. 1997. Sensitivity measures, ANOVA-like techniques and the use of bootstrap. J. Stat. Comput. Simul.
58:99–120. doi:10.1080/00949659708811825
Asseng, S., C. Richter, and G. Wessolek. 1997. Modelling root growth ofwheat as the linkage between crop and soil. Plant Soil 190:267–277.
doi:10.1023/A:1004228201299
Bingham, I.J., and L. Wu. 2011. Simulation of wheat growth using the 3D
root architecture model SPACSYS: Validation and sensitivity analysis.
Eur. J. Agron. 34:181–189. doi:10.1016/j.eja.2011.01.003

VZJ | Advancing Critical Zone SciencePrasad, P.V.V., S.A. Staggenborg, and Z. Ristic. 2008. Impacts of drought
and/or heat stress on physiological, development, growth, and yieldprocesses of crop plants. In: L.R. Ahuja et al., editors, Response of crops
to limited water: Understanding and modeling water stress effects onplant growth processes. Adv. Agric. Syst. Model. Ser. 1. ASA, CSSA, and
SSSA, Madison, WI. p. 301–355. doi:10.2134/advagricsystmodel1.c11
Reddy, K.R., and V.G. Kakani. 2007. Screening Capsicum species ofdifferent origins for high temperature tolerance by in vitro pollen germination and pollen tube length. Sci. Hortic. 112:130–135.
doi:10.1016/j.scienta.2006.12.014

p. 15 of 16

Rewald, B., E. Raveh, T. Gendler, J.E. Ephrath, and S. Rachmilevitch. 2012.
Phenotypic plasticity and water flux rates of Citrus root orders undersalinity. J. Exp. Bot. 63:2717–2727. doi:10.1093/jxb/err457

periments and a model simulation. J. Plant Growth Regul. 21:441–458.
doi:10.1007/s00344-003-0006-y

Robertson, M.J., S. Fukai, G.L. Hammer, and M.M. Ludlow. 1993. Modellingroot growth of grain sorghum using the CERES approach. Field Crops
Res. 33:113–130. doi:10.1016/0378-4290(93)90097-7

Valenzuela-Estrada, L.R., V. Vera-Caraballo, L.E. Ruth, and D.M. Eissenstat.
2008. Root anatomy, morphology, and longevity among root orders
in Vaccinium corymbosum (Ericaceae). Am. J. Bot. 95:1506–1514.
doi:10.3732/ajb.0800092

Rosolem, R., H.V. Gupta, W.J. Shuttleworth, and X. Zeng. 2012. A fully
multiple-criteria implementation for parameter sensitivity analysis. J.
Geophys. Res. 117:D07103. doi:10.1029/2011JD016355

van Genuchten, M.Th. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–
898. doi:10.2136/sssaj1980.03615995004400050002x

Saha, S.L., M.M. Hossain, M.M. Rahman, C.G. Kuo, and S. Abdullah. 2010.
Effect of high temperature stress on the performance of twelve sweet
pepper genotypes. Bangladesh J. Agric. Res. 35:525–534.

Vrugt, J.A. 2016. Markov chain Monte Carlo simulation using the DREAMsoftware package: Theory, concepts, and MATLAB implementation.
Environ. Model. Softw. 75:273–316. doi:10.1016/j.envsoft.2015.08.013

Saltelli, A., P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and S. Tarantola. 2010. Variance based sensitivity analysis of model output: Designand estimator for the total sensitivity index. Comput. Phys. Commun.
181:259–270. doi:10.1016/j.cpc.2009.09.018

Vrugt, J.A., J.W. Hopmans, and J. Šimůnek. 2001a. Calibration of a two-dimensional root water uptake model. Soil Sci. Soc. Am. J. 65:1027–1037.
doi:10.2136/sssaj2001.6541027x

Seidel, S.J., S. Rachmilevitch, N. Schütze, and N. Lazarovitch. 2016. Modelling the impact of drought and heat stress on common bean with two
different photosynthesis model approaches. Environ. Modell. Softw.
81:111–121. doi:10.1016/j.envsoft.2016.04.001

Vrugt, J.A., C.J.F. ter Braak, C.G.H. Diks, B.A. Robinson, J.M. Hyman, and
D. Higdon. 2009. Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling. Int. J. Nonlinear Sci. Numer. Simul. 10:271–288.
doi:10.1515/IJNSNS.2009.10.3.273

Šimůnek, J., and J.W. Hopmans. 2009. Modeling compensated root water and nutrient uptake. Ecol. Modell. 220:505–521. doi:10.1016/j.ecolmodel.2008.11.004

Vrugt, J.A., M.T. Wijk, J.W. Hopmans, and J. Šimůnek. 2001b. One-, two-,and three-dimensional root water uptake functions for transient modeling. Water Resour. Res. 37:2457–2470. doi:10.1029/2000WR000027

Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2008. Development andapplications of the HYDRUS and STANMOD software packages andrelated codes. Vadose Zone J. 7:587–600. doi:10.2136/vzj2007.0077

Wang, Z., J. Li, and Y. Li. 2016. Assessing the effects of drip irrigation system uniformity and spatial variability in soil on nitrate leaching throughsimulation. Trans. ASABE 59:279–290. doi:10.13031/trans.59.11488

Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2016. Recent developments and applications of the HYDRUS computer software packages.
Vadose Zone J. 15(7). doi:10.2136/vzj2016.04.0033

Williams, J.R., C.A. Jones, J.R. Kiniry, and D.A. Spanel. 1989. The EPIC crop
growth model. Trans. ASAE 32:497–511. doi:10.13031/2013.31032

Smith, D.M., N.G. Inman-Bamber, and P.J. Thornburn. 2005. Growth andfunction of the sugarcane root system. Field Crops Res. 92:169–183.
doi:10.1016/j.fcr.2005.01.017

Wu, L., M.B. McGechan, C.A. Watson, and J.A. Baddeley. 2005. Developing existing plant root system architecture models tomeet future agriculture challenges. Adv. Agron. 85:181–219.
doi:10.1016/S0065-2113(04)85004-1

Sobol¢, I.M. 1993. Sensitivity estimates for nonlinear mathematical models.
Math. Modell. Comput. Exp. 1:407–414.

Wu, L., R. Zhang, and S. Gui. 1999. Modeling soil water movement with water uptake by roots. Plant Soil 215:7–17. doi:10.1023/A:1004702807951

Somma, F., V. Clausnitzer, and J.W. Hopmans. 1998. Modeling of transient three-dimensional soil water and solute transport with root
growth and water and nutrient uptake. Plant Soil 202:281–293.
doi:10.1023/A:1004378602378

Wu, W., and S. Cheng. 2014. Root genetic research, an opportunityand challenge to rice improvement. Field Crops Res. 165:111–124.
doi:10.1016/j.fcr.2014.04.013

Song, X., J. Zhang, C. Zhan, Y. Xuan, M. Ye, and C. Xu. 2015. Global sensitivity analysis in hydrological modeling: Review of concepts, methods, theoretical framework, and applications. J. Hydrol. 523:739–757.
doi:10.1016/j.jhydrol.2015.02.013
Stone, J.A., T.C. Kaspar, and H.M. Taylor. 1983. Predicting soybean rooting depth as a function of soil temperature. Agron. J. 75:1050–1054.
doi:10.2134/agronj1983.00021962007500060042x
Stone, J.A., and H.M. Taylor. 1983. Temperature and the development of
the taproot and lateral roots of four indeterminate soybean cultivars.
Agron. J. 75:613–618. doi:10.2134/agronj1983.00021962007500040010x
Taiz, L., and E. Zeiger. 2002. Plant physiology. 3rd ed. Sinauer Assoc., Sunderland, MA.
Tsutsumi, D., K. Kosugi, and T. Mizuyama. 2003. Effect of hydrotropism
on root system development in soybean (Glycine max): Growth ex-

VZJ | Advancing Critical Zone ScienceYan, W., and L.A. Hunt. 1999. An equation for modelling the temperature response of plants using only the cardinal temperatures. Ann. Bot.
84:607–614. doi:10.1006/anbo.1999.0955
Yin, X., M.J. Kropff, G. McLaren, and R.M. Visperas. 1995. A nonlinear
model for crop development as a function of temperature. Agric. For.
Meteorol. 77:1–16. doi:10.1016/0168-1923(95)02236-Q
Yu, G.-R., J. Zhuang, K. Nakayama, and J. Yan. 2007. Root water uptake
and profile soil water as affected by vertical root distribution. Plant
Ecol. 189:15–30. doi:10.1007/s11258-006-9163-yZhang, C., J. Chu, and G. Fu. 2013. Sobol¢’s sensitivity analysis for a distributed hydrological model of Yichun River Basin, China. J. Hydrol.
480:58–68. doi:10.1016/j.jhydrol.2012.12.005
Zuo, Q., F. Jie, R. Zhang, and L. Meng. 2004. A generalized function ofwheat’s root length density distributions. Vadose Zone J. 3:271–277.
doi:10.2136/vzj2004.2710

p. 16 of 16