RESEARCH ARTICLE
10.1029/2019WR025432
Key Points:• A dynamic plant uptake module iscoupled with the HYDRUS model tosimulate the translocation ofchemicals in the soil‐plantcontinuum• Theoretical and experimental
scenarios are used to validate thecoupled model and demonstrate its
ﬂexibility• The model satisfactorily reproduces
the measured carbamazepine
translocation and transformation inthree vegetables

Correspondence to:
G. Brunetti,
giuseppe.brunetti@boku.ac.at

Citation:
Brunetti, G., Kodešová, R., &Šimůnek, J. (2019). Modeling thetranslocation and transformation ofchemicals in the soil‐plant continuum:
A dynamic plant uptake module for the
hydrus model. Water Resources
Research, 55, 8967–8989. https://doi.
org/10.1029/2019WR025432
Received 25 APR 2019
Accepted 9 SEP 2019
Accepted article online 15 SEP 2019
Published online 13 NOV 2019

Modeling the Translocation and Transformationof Chemicals in the Soil‐Plant Continuum:
A Dynamic Plant Uptake Module
for the HYDRUS Model
Giuseppe Brunetti1, Radka Kodešová2, and Jiří Šimůnek3

1

Institute of Hydraulics and Rural Water Management, University of Natural Resources and Life Sciences, Vienna,
Austria, 2Faculty of Agrobiology, Food and Natural Resources, Department of Soil Science and Soil Protection, Czech
University of Life Sciences Prague, Prague, Czech Republic, 3Department of Environmental Sciences, University ofCalifornia, Riverside, CA, USA

Abstract Food contamination is responsible for thousands of deaths worldwide every year. Plants
represent the most common pathway for chemicals into the human and animal food chain. Although
existing dynamic plant uptake models for chemicals are crucial for the development of reliable mitigation
strategies for food pollution, they nevertheless simplify the description of physicochemical processes in soil
and plants, mass transfer processes between soil and plants and in plants, and transformation in plants. To
ﬁll this scientiﬁc gap, we couple a widely used hydrological model (HYDRUS) with a multicompartment
dynamic plant uptake model, which accounts for differentiated multiple metabolization pathways in plant's
tissues. The developed model is validated ﬁrst theoretically and then experimentally against measured datafrom an experiment on the translocation and transformation of carbamazepine in three vegetables. The
analysis is further enriched by performing a global sensitivity analysis on the soil‐plant model to identify
factors driving the compound's accumulation in plants' shoots, as well as to elucidate the role and theimportance of soil hydraulic properties on the plant uptake process. Results of the multilevel numericalanalysis emphasize the model's ﬂexibility and demonstrate its ability to accurately reproduce
physicochemical processes involved in the dynamic plant uptake of chemicals from contaminated soils.
1. Introduction
Food contamination is a direct consequence of environmental pollution and represents a major risk forhuman health. Recently, the Food and Agricultural Organization of the United Nations (Food andAgricultural Organization, 2019) reported that food‐contaminated natural toxins or chemicals are responsible for more than 339,000 people falling ill and 20,000 dying worldwide every year. Plants represent the mostcommon pathway into the human and animal food chain for environmental contaminants. Chemicals from
polluted environments are uptaken by plants, where they are bioaccumulated and eventually metabolized in
active by‐products (e.g., Mattina et al., 2003; Sabourin et al., 2012; Zhang et al., 2009), some of which can be
dangerous for human health (e.g., Kodešová et al., 2019a; Malchi et al., 2014; So et al., 1994; Warner et al.,1992). Thus, widespread pollution forces the scientiﬁc community to better understand physicochemical
processes within the soil‐plant‐atmosphere continuum in order to develop reliable mitigation strategies
for food risk assessment and food safety.

©2019. American Geophysical Union.
All Rights Reserved.

BRUNETTI ET AL.

Numerical models certainly play a crucial role in reaching this goal. Several plant uptake models have beendeveloped and tested in the past, ranging from empirical to mechanistic models. Topp et al. (1986) and Travis
and Arms (1988) proposed two different empirical models correlating the concentration of a chemical compound in the plant's tissues with its physicochemical properties. Despite the advantages of simple analytical
relationships, these models provide a rather coarse description of the plant uptake processes. Alternatively,
Ryan et al. (1988) developed a semimechanistic screening model for assessing the uptake of nonionic chemicals from soil. However, due to its theoretical limitations (e.g., neglecting plant‐air exchange, metabolism,
and growth dilution), this model was only intended to provide a simple approach for screening chemicals
for which plant uptake may be an important pathway for human exposure. Hung and Mackay (1997) proposed a classic mechanistic approach instead, describing the kinetics of uptake in leaves, stem, and roots
8967

Water Resources Research

10.1029/2019WR025432

using a dynamic multicompartment model, which included xylem and phloem ﬂows, gaseous exchange,
growth dilution, and metabolization.
Similarly, Trapp (2007) and Rein et al. (2011) developed a mechanistic multicompartment model encompassing all major physicochemical processes happening during the dynamic plant uptake (DPU) process.
The model, initially developed for neutral compounds, was later extended to ionic chemicals and tested successfully against measured data under different operating conditions (Trapp, 2009). Hurtado et al. (2016)used the DPU model to inversely estimate kinetic parameters from a greenhouse experiment on dynamic
uptake of emerging contaminants in lettuce. In another study, Trapp (2015) successfully calibrated the
DPU model to simulate uptake of chlorinated organic compounds in radish.
Despite acceptable results, existing DPU models generally lack an accurate description of physicochemical
processes in soils, which are among the most important sources of pollution for plants (Eugenio et al.,2018). The variably saturated hydraulic behavior of the root zone and biogeochemical processes in thenear‐surface environment, such as solute sorption and degradation, strongly affect actual transpiration patterns and solute uptake across the root surfaces. The accurate numerical description of these processes isthus crucial for the development of reliable soil‐plant models. Along these lines, Legind et al. (2012) andTrapp and Eggen (2013) coupled the multicompartment model of Trapp (2007) with a model of water andsolute transport in soil based on the “tipping buckets” approach. Such a numerical approach usuallyassumes that, due to gravity, water can only percolate downwards as Darcian ﬂow (e.g., Schroeder et al.,1994) and simpliﬁes upward ﬂow due to capillarity. Additionally, this approach requires several thresholdvalues as input, including ﬁeld capacity and wilting point, that are often poorly deﬁned (Twarakavi et al.,2009). In particular, the physical meaning of ﬁeld capacity is ambiguous (Twarakavi et al., 2009) and should
be adjusted to site‐speciﬁc conditions, considering its high impact on simulated transport processes.
Due to these theoretical limitations (Emerman, 1995; Scanlon et al., 2003), this approach is mostly suited forspeciﬁc working conditions in which water ﬂow is mainly driven by gravity (i.e., irrigated agriculture withpredominantly downward vertical ﬂow) and is not well suited for many other important environments inwhich capillarity plays an important role (e.g., rainfed agriculture and/or upward ﬂow from shallow groundwater). This approach is often used when water ﬂow in soils is a secondary process in the model, whichinstead focuses on providing a detailed description of a different primary process (e.g., crop growth [Booteet al., 2008; Shelia et al., 2018] and DPU of chemicals [Trapp, 2007], etc.). Its application in traditional modeling scenarios (i.e., uptake of chemicals from contaminated soils under transient atmospheric conditions)
may thus result in poor model generalizability stemming from the theoretical model's inadequacy. On theother hand, hydrological models provide sophisticated hydraulic and hydrologic descriptions while simplifying other model components. Therefore, research efforts should focus on the development of integrated
soil‐plant models which include a more coherent description of physicochemical processes in both soil
and plant domains (Vereecken et al., 2016).
Numerical models from vadose zone hydrology, which have been used extensively in the last decades (Jarvis& Larsbo, 2012; Šimůnek et al., 2016), offer great opportunities to head in this direction. In particular,
Richards‐based approaches have been used successfully in many applications and scientiﬁc ﬁelds (e.g.,
Brunetti et al., 2017, 2019; Brunetti, Porti, & Piro, 2018; Brunetti, Šimůnek, & Bautista, 2018; Cheyns
et al., 2010, Hanson et al., 2006, Jellali et al., 2016, Scanlon et al., 2003). Nolan et al. (2005) compared theperformance of the Richards‐based and tipping bucket models in reproducing measured breakthroughcurves of bromide and atrazine from two experimental facilities in Merced, California, and White River
Basin, Indiana, respectively. Results of this study indicated a better accuracy of the Richards‐based models(i.e., HYDRUS‐2D, LEACHP, RZWQM, and VS2DT) than the tipping‐bucket models (i.e., CALF,
GLEAMS, and PRZM). Furthermore, the authors reported a numerical inﬂuence of grid cell thickness on
the effective solute dispersivity in some tipping bucket models. In another study, Herbst et al. (2005) concluded that the Richards' equation‐based models should perform better than the capacity‐based approaches
in terms of predicting soil moisture contents, drainage amounts, and actual evapotranspiration. Similar conclusions were also obtained in the model intercomparison study by Scanlon, Christman, et al. (2003).
Besides simplifying soil processes, it must be emphasized that existing DPU models also simplify metabolization chains in the plant's tissues. The green liver behavior of plants (Sandermann, 1994) is generally
accounted for using a single ﬁrst‐order decay coefﬁcient that summarizes enzymatic biodegradation of
BRUNETTI ET AL.

8968

Water Resources Research

10.1029/2019WR025432

compounds in the plant. Although obtaining comprehensive measurements of metabolites may be very difﬁcult, such an approach strongly reduces the model ﬂexibility due to its inability to describe the generationof metabolites usually encountered in actual conditions (Kodešová et al., 2019b). On the other hand, the useof parent‐metabolites data is becoming more popular in other scientiﬁc ﬁelds, such as pharmacokinetic
modeling (Bertrand et al., 2011; Stroh et al., 2013), that share some theoretical and practical similarities with
the DPU models. Such joint parent‐metabolites pharmacokinetics models are used to simulate transformation and interaction effects and can be used in conjunction with other numerical techniques (e.g., Dumont
et al., 2013) to increase the model predictive capability. This emphasizes the need to further reﬁne thenumerical description of degradation chains in plants.
We aim to ﬁll the aforementioned scientiﬁc gaps by developing a fully coupled soil‐plant model capable ofproviding a mechanistic comprehensive description of transport and reaction processes in the soil‐plantcontinuum. To this end, we couple the widely used hydrological model, HYDRUS (Šimůnek et al.,2016), with a modiﬁed version of the multicompartment DPU model proposed by Trapp (2007), that
now accounts for differentiated multiple metabolization pathways in plant's tissues. HYDRUS is a
Richards equation‐based model that has been chosen due to its reliability, versatility, success in several
model intercomparison studies (e.g., Nolan et al., 2005; Scanlon, Christman, et al., 2003), and widespread
use among hydrologists (Šimůnek et al., 2016). The problem is addressed in the following way. First, theoretical background to numerically combine the soil and plant models is developed. Then, the modiﬁed
DPU model is theoretically validated against the original implementation of Trapp (2007) and internallycoupled with the HYDRUS code to minimize the computational cost. Finally, the developed soil‐plantmodel is validated against measured data from an experiment on the translocation and transformationof carbamazepine (CBZ) in three types of vegetables. The analysis is further enriched by performing a global sensitivity analysis on the soil‐plant model to identify factors driving the compound's accumulation inplants' shoots, as well as to clarify the role and importance of soil hydraulic properties on the plant
uptake process.

2. Materials and Methods
2.1. An Integrated Soil‐Plant Modeling Environment
To describe the fate of neutral compounds in the soil‐plant system (Figure 1), we couple a widely usedhydrological model (HYDRUS) with a multicompartment DPU model. While HYDRUS (Šimůnek
et al., 2016) is a ﬁnite‐element model for simulating the movement of water, heat, and multiple solutes
in variably saturated porous media, the DPU model accounts for the movement of water and solutes
through the plant while also accounting for differentiated multiple metabolization pathways in plant's
tissues. Additionally, the HYDRUS model has a relatively sophisticated root water and nutrient uptake
model (Šimůnek & Hopmans, 2009), which provides a link between HYDRUS and DPU. Figure 1 summarizes the physicochemical processes in the soil‐plant continuum and their theoretical conceptualization in the coupled soil‐plant model.
2.2. DPU of Neutral Compounds
In this section, we provide an overall description of the model structure and a theoretical background forphysical and biogeochemical processes implemented in the DPU model.
2.2.1. Model Structure
The numerical model is based on the multicompartmental approach proposed by Trapp (2007) and further
developed by Rein et al. (2011) and Legind et al. (2011). The model simulates the translocation and biodegradation of N neutral compounds in plants, which are conceptualized in four compartments: roots, stem,
leaves, and fruits. The choice to focus only on neutral compounds is intended to be only a ﬁrst step in themodel development. Transport processes consist of advection (with xylem ﬂow) and diffusive losses or gains,while reaction processes mainly represent the enzymatic biodegradation in plant's tissue. More speciﬁcally,
the following processes are considered in the DPU model:• Translocation of compounds from the roots to the stem and from the stem to leaves and fruits with thetranspiration stream. The phloem ﬂux is neglected.
• Volatilization of compounds in the stem, leaves, and fruits.
• Gaseous and particle deposition from air to the stem, leaves, and fruits.
BRUNETTI ET AL.

8969

Water Resources Research

10.1029/2019WR025432

Figure 1. A schematic of the physicochemical processes in the (a) soil‐plant continuum and their theoretical conceptualization in the (b) coupled soil‐plant model.

• Compounds dilution by plant growth in all compartments.
• Compounds metabolization in all compartments.
The DPU model calculates the compounds' concentration in each compartment. The concentration matrix C
has N rows (for each compound, represented below by a subscript i) and four columns (for each compartment, i.e., 1 = roots, 2 = stem, 3 = leaves, and 4 = fruits, represented below by a subscript j) and contains
concentrations expressed as mass of each compound per fresh weight in different plants' compartments(Ms/Mfw, where subscript s refers to solutes and fw to fresh weight).
Transport and reaction processes depend on physical and chemical properties of the plant and compounds,
respectively. In particular, the characteristics of four plant tissues are summarized in seven vectors.
8
8 gr 9
8 9
9
9
8
8 9
8 max 9
8 9
M1 >
L1 >
W1 >
ρ1 >
>
>
>
M 01 >
K1 >
SA1 >
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
=
<W =
<L =
< M max =<ρ >< M0 =< K gr =< SA =
2
2
2
2
2
2
2
0
max
gr
A
;L ¼
;M ¼
;Κ ¼
; (1)
W¼
¼
;M
;S ¼
;ρ ¼
gr
max
0
A
>
>
>
>
>
>
>
>
M3 >
L3 >
ρ3 >
M3 >
K3 >
S3 >
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> W3 >
>
>
>
>
>
>
>
>
>
;
;
;
;
:
: >
:
:
>
>
>
>
>
>
: 0;
: gr ;
: A;
M max
W4
L4
ρ4
M4
K4
S4
4

where W is the water content of the plant (L3/Mfw), L is the lipid content of the plant (Ml/Mfw, where subscript l refers to the lipid content), M0 is the initial plant mass (Mfw), Mmax is the maximum plant mass(Mfw), Kgr is the growth rate coefﬁcient (T−1), SA is the plant's speciﬁc area (L2Mfw−1), and ρ is the plant's
density (Mfw/L3). On the other hand, compounds' properties are described by four vectors.
8
9
8
8
8
9
9
9
logK OW1 >K AW1 >C A1 >
>
m1 >
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
< logK<K<C >< m2 >
=
=
=
=
OW2
AW2
A2
KOW ¼
; KAW ¼
; CA ¼
;m ¼
;
>
>
>
>
>
:
: >
: >
: >
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
:
:
:
>
;
;
;
>
: N;
logK OWN
K AWN
C AN
m

(2)

where KOW is the n‐octanol‐water partition coefﬁcient (−), KAW is the air‐water partition coefﬁcient (L3/L3),
CA is the compound concentration in the air (Ms/L3), and m is the molar mass (M). Other required inputs
are the air temperature T (K) and the relative humidity ϕ (−).
2.2.2. Mass Balance
The mass balance of N compounds in the four plant's compartments is mathematically described using a system of ordinary differential equations that is written in a compact form as
dC¼ QIN þ QDEP þ QGAS −QTR −QVOL þ Ω;
dt

(3)

where QIN is the inﬂow ﬂux (Ms·Mfw−1·T−1; for each compound and each compartment), QDEP is the deposition of a chemical from particle deposition, QGAS is chemical gaseous uptake from the air, QTR is the translocation ﬂux between compartments, QVOL is the volatilization ﬂux, and Ω represents reactions.
BRUNETTI ET AL.

8970

Water Resources Research

10.1029/2019WR025432

In the numerical model, equation (3) is solved implicitly in time using the backward Euler method (Butcher,2016) to increase the stability of the numerical solution.
CtþΔt −Ct
tþΔt
tþΔt
tþΔt
tþΔt
tþΔt
:
¼ QtþΔtIN þ QDEP þ QGAS −QTR −QVOL þ Ω
Δt(4)

Equation (4) is expressed as a linear equation system, which is solved using the iterative quasi‐Newton
Broyden's method (Broyden, 1969) at each time step.
2.2.3. Plant Growth
Many annual crops exhibit a logistic growth characterized by an exponential growth at the beginning and an
asymptotic behavior towards ripening. Accordingly, the change of the plant mass in each compartment M(Mfw) as a function of time t (T) is expressed as


dM
M
:
¼ Kgr M 1−
dt
Mmax

(5)

The plant mass in each compartment can be calculated by integrating the logistic growth function as follows:Mðt Þ ¼

1þ



Mmax

:
exp−Kgr t

Mmax
M0 −1

(6)

The plant surface area A (L2) is related to the plant mass M through the speciﬁc surface area SA (L2/Mfw).
Aj ðt Þ ¼ SAj M j ðt Þ

(7)

where index j refers to four plant compartments. The model can alternatively use measured plant growthdata and interpolate between measured data points.
2.2.4. Chemical Equilibrium and Partitioning
The typical time scale for chemical diffusion of neutral compounds through a cell membrane is negligible
compared to the time scale of soil‐plant processes. Under such circumstances, the rate of the sorption process
is fast enough to consider instantaneous equilibrium, and the capacity of a chemical to diffuse across membranes can be described by the equilibrium partition coefﬁcient. In particular, the equilibrium partitioning
between a hydrophobic phase (lipids, oils, etc.) and water is described by the n‐octanol‐water partition coefﬁcient KOW (−). Measured values are available for many compounds, and they are usually expressed using alogarithmic scale. Similarly, the coefﬁcient KAW (−) describes the partitioning between liquid and air.
Briggs et al. (1982) deﬁned a plant‐water partition coefﬁcient KPW (L3Mfw−1), which is expressed in its
matrix form as
KPWij ¼ W j þ Lj aK bOWi

i ¼ 1; N and j ¼ 1; 4;(8)where Wj and Lj are the water (L3/Mfw) and lipid contents of the plant compartments (Ml/Mfw), respectively,
a is a coefﬁcient (L3/Ml) usually assumed 0.00122 m/kgl (Trapp, 2007), and b is an exponent equal to 0.77
and 0.95 for roots and plant shoots, respectively. Indices i and j refer to matrix columns (to N compounds)
and rows (four compartments), respectively. It must be emphasized that the accuracy of equation (8) in predicting the sorption mechanism in the plant decreases with the polarity of the investigated chemical compound (Li et al., 2005). Therefore, for highly polar compounds, it is recommended to experimentally
characterize the sorption mechanism.
2.2.5. Translocation
Most plants have two transport systems inside: xylem and phloem. Xylem is a nonliving vascular system thatis hydraulically connected to phloem (Brodersen et al., 2019). Water is drawn upwards by physical forces.
Phloem is composed of living cells, the sieve tubes. The ﬂux in phloem generally occurs from the leaves tothe fruits. In particular, phloem transports sugar and other essential substances. Xylem is responsible forthe water translocation from the roots to other plant compartments induced by plant transpiration.
BRUNETTI ET AL.

8971

Water Resources Research

10.1029/2019WR025432

The effect of the plant water storage is neglected (e.g, Javaux et al., 2013; Manoli et al., 2014). Such assumption, which is valid only for small crops (Hartzell et al., 2017; Vogel et al., 2017), implies that the xylem ﬂow
QXYL (L3/T) is equal to the transpiration stream through the roots and stem, QTP, while it is partitioned in
fruits and leaves using corresponding surface areas (Trapp, 2007).
QXYL ¼

A4 QTP
A4 QTP:QTP QTP QTP −A4 þ A3 A4 þ A3

(9)

The solute translocation ﬂux QTR (Ms·Mfw−1·T−1) is expressed as
QTRij ¼

QXYLij Cij:
M j K PWij

(10)

Leaves and fruits are the end of the xylem system; and thus, no further translocation is assumed in thesecompartments. Such an assumption implies that to facilitate the model implementation, we can similarly
deﬁne for each compartment the solute inﬂow QIN (Ms·Mfw−1·T−1), which contains the compounds. In particular, solute inﬂow into the root compartment is represented by actual solute uptake calculated by
HYDRUS for a soil surface area of 1 m2 divided by the roots mass (Ms·Mfw−1·T−1), while the translocation
ﬂux from the roots is the inﬂow for the stem. Finally, the xylem ﬂow from the stem is divided into two
streams that enter leaves and fruits.
2.2.6. Exchange Plant‐Air: Plant Permeability
The plant‐air exchange is regulated by the permeability of the plant's components (Trapp, 2007). A chemicalin the air has to overcome only the air boundary layer between the leaf surface and the turbulent atmosphere. It can then adsorb on the leaf surface, which is a waxy layer known as the cuticle. On the other hand,a chemical transported by the transpiration stream into the leaf can be exchanged with air through cuticles
or stomata. Thus, three different phases provide resistance (aqueous phase, lipid phase, and gas phase) andtwo pathways are in parallel (stomatal and cuticular pathways).
The more lipophilic a compound, the higher its permeability across the cuticle barrier.
PCi ¼ 100:704logK OWi −11:2 ;(11)

where PC is the cuticle conductance (L/T). The next resistance is provided by the air boundary layer surrounding the plant, Pair (L/T). Resistance of 200 s/m was estimated as typical for a chemical with molar
weight m = 300 g/mol (Thompson, 1983). The conductance of the air boundary layer, Pair, for a chemicalwith a molar mass m is thus
Pairi ¼

pﬃﬃﬃﬃﬃﬃﬃﬃ
K AWi 300
pﬃﬃﬃﬃﬃﬃ ;200 mi(12)

where m is the molar mass of the chemical. However, water does need to cross the aqueous layer in the apoplast. The permeability of this layer Paqua (L/T) issﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Paquai ¼ DO2

32

mi

z;

(13)

where z is the thickness of the aqueous layer assumed 0.0005 m, and DO2 is the diffusion coefﬁcient of oxygenin water (2x10−9 m2/s). Thus, the total cuticle permeability PC,tot (L/T) is
1
1
1 :
þ
PC
Pair þ Paqua

PC;tot ¼ 1(14)

The role of the stomata is gas exchange. While the stomata are open, carbon dioxide is taken up, and water
vapor is lost. Details on the calculation of transpiration are provided in the following sections. The conductance of the stomata, PS (L/T), can thus be calculated from the xylem ﬂow QXYL as
BRUNETTI ET AL.

8972

Water Resources Research

10.1029/2019WR025432

pﬃﬃﬃﬃﬃ
18
PSij ¼
pﬃﬃﬃﬃﬃﬃ K AWi ;m
Aj C H2O;sat −ϕCH2O;sati
QXYLj ρwater(15)

where ϕ is the relative air humidity (−), ρwater is the water density (M/L3), and CH2O,sat is the water concentration inside plant's shoots (M/L3), which can be calculated as
CH2O;sat ¼

pH2O;
461:9T(16)

where pH2O is the saturation vapor pressure (ML−1·T−2) at a given temperature T (K), calculated using theempirical Magnus equation:
7:5ðT−273:15Þ

pH2O;sat ¼ 610:7×10237þðT−273:15Þ :(17)

The exchange of the chemical through cuticle and stomata occurs in parallel, and thus, permeabilities areadded to derive the total permeability for the exchange between the plant and air, PP (L/T):
PP ¼ PS þ PC;tot :(18)

2.2.7. Volatilization
Chemicals can leave or enter the plant in several ways. This exchange process is regulated by the plant‐air
partition vector KPA (L3/L3), which is expressed as
KPAij ¼

K PWij ρj
:
K AWi

(19)

No volatilization is assumed in the roots since they are generally not exposed to air. The volatilization ratedepends on the permeability of plant's compartments, PP. In the developed model, only the cuticular pathway is used for the stem, while for fruits and leaves, exchanges of the chemical occur through the cuticle andstomata in parallel (Rein et al., 2011). The volatilization ﬂux vector QVOL (Ms·Mfw−1·T−1) is expressed as
QVOLij ¼

Aj PP;ij ρj Cij:
M j K PAij

(20)

2.2.8. Gaseous Uptake and Particle Deposition
The gaseous uptake ﬂux of chemicals from the air, QGAS (Ms·Mfw−1·T−1), is deﬁned as:

QGASij ¼


Aj PP;ij 1−f p CAi
Mj

;

(21)

where fp is the fraction of a chemical adsorbed on the particles (−). The uptake of chemical from particle
deposition, QDEP (Ms·Mfw−1·T−1), is:
QDEPij ¼

Aj vdep f p CAi
;
Mj

(22)

where vdep (L/T) is the particle deposition velocity assumed 0.001 m/s (Trapp & Matthies, 1998).
2.2.9. Metabolization of Neutral Compounds
Plant compounds' metabolism is divided into three stages: transformation (I), conjugation (II), and compartmentation and storage processes (III), although plants have no excretion systems. Cellular storage sites are
the vacuole and the cell wall (Sandermann, 1992). Plants thus have the ability to metabolize chemicals by
exploiting the enzyme degradation. Furthermore, it is well known that millions of bacteria live on, but also
inside of, plants. Most of these bacteria are heterotrophic and use chemicals as a substrate to grow.
Compounds can be used as an electron acceptor, an electron donor, as an energy source, or as a precursor
for other molecules. This opens opportunities—plants can be “vaccinated” with degrader bacteria to degrade
pollutants (phytoremediation) or pesticides (Barac et al., 2004; Glick, 2010).
BRUNETTI ET AL.

8973

Water Resources Research

10.1029/2019WR025432

Figure 2. Examples of (a) double, (b)sequential, and (c) closed metabolization chains of three compounds in the jth compartment (in this case j = 1, i.e., in the roots). Black circles indicate an incomplete conversion of the reactant in by‐
products.

Nevertheless, traditional modeling approaches are based on the use of a single ﬁrst‐order rate coefﬁcient todescribe compound's biodegradation (e.g., Legind et al., 2011; Rein et al., 2011; Trapp, 2007). Such an
approach completely neglects the transformation chain in the plant's tissues, which leads to the formationof new chemical by‐products. For instance, CBZ can generate both hydroxylated and conjugated metabolites,
which can be further degraded (Goldstein et al., 2018). Furthermore, not all compounds are generally
degraded since plants have an upper limit on the compound compartmentation. These metabolites, which
can be chemically active, are important for risk assessment (e.g., So et al., 1994; Warner et al., 1992) and must
be considered in the analysis. In addition to that, recent studies have shown differentiated biodegradation
rates in different plant's compartments (e.g., Fantke & Juraske, 2013; Kodešová et al., 2019b), suggesting
the need to use a more ﬂexible mathematical description. To overcome such limitations in the present model,
the plant metabolization is modeled using multiple metabolization matrices Γj. These matrices consist ofmultiple ﬁrst‐order rate coefﬁcients, which provide connections between parent and daughter species. The
matrix Γj is a square matrix of a dimension N (i.e., a number of metabolites considered).
8>< −τ 11
Γj ¼
:
>
: 0
þτ N1

9
0
þτ 1N >=
N
with ∑ Γ ij ≤0:
⋅
:
>
i¼1
;
:: −τ NN::(23)

where τ are ﬁrst‐order rate constants for the compound (T−1), while τ′ are similar ﬁrst‐order rate constants
providing connections between individual chain species (T−1). The metabolization rate vector for the jthcompartment, Ωj (Ms·Mfw−1·T−1), is thus deﬁned as the dot product of the metabolization and concentration
matrices.S
8>< −τ 11
Ωj ¼ Γj ⋅Cj ¼
:
>
: 0
þτ N1

98
9
0
þτ 1N >< C1j >
=>
=
:
⋅
:
⋅
:
>
>
:
;>
;
CNj
:: −τ NN::(24)

The proposed approach guarantees high modeling ﬂexibility for handling multiple nonsequential degradation chains and simulating the effect of reaction yields. Figure 2 shows a schematic of three metabolization
chains in the jth compartment with their respective model implementations.
2.3. Implementation in the HYDRUS Model
The compounds' translocation in the plant is mainly driven by actual transpiration and solute uptake ﬂuxes,
which are in turn dependent on the pressure head and solute distribution in the root zone (de Jong van Lier
BRUNETTI ET AL.

8974

Water Resources Research

10.1029/2019WR025432

et al., 2013), respectively. Huber et al. (2014) demonstrated that soil moisture heterogeneity strongly inﬂuences plant hydraulics. Hence, an accurate description of the physicochemical processes in the near‐surface
environment is crucial for developing a reliable soil‐plant modeling tool (Vereecken et al., 2016). To this aim,
the mechanistic model HYDRUS is coupled with the DPU model to further extend its modeling capabilities
to the plant domain thus obtaining a fully integrated modeling environment for contaminant risk assessment. In this section, we ﬁrst provide an overview of the HYDRUS model and then a thorough descriptionof the coupling strategy and practical implementation of the coupled soil‐plant‐atmosphere model.
2.3.1. HYDRUS ModelHYDRUS (Šimůnek et al., 2016) is a ﬁnite‐element model for simulating the movement of water, heat, andmultiple solutes in variably saturated porous media. It numerically solves the Richards equation for multidimensional variably saturated water ﬂow.
∂θ
¼ ∇½K⋅∇ðh−zÞ−SðhÞ;
∂t

(25)

where t is time (T), z is the vertical coordinate (L), θ is the volumetric water content (L3/L3), K is the unsaturated hydraulic conductivity (L/T), h is the pressure head (L), and S is a sink term representing root water
uptake (T−1).
Solute transport can be simulated in the liquid, solid, and gaseous phases, though the latter is not of
immediate interest in soil‐plant modeling. Therefore, the transport of the ith solute in the soil is describedusing the advection‐dispersion‐reaction equation, assuming that solutes can exist only in the solid andliquid phases.
∂θCi ∂ρsiþ
¼ ∇ θDW
i ⋅∇C i −∇ðqC i Þ−r a;i −ϕi ;
∂t
∂t

(26)

where Ci is the concentration of the ith solute in the liquid phase (M/L), si is the concentration of the ith
solute in the solid phase (M/M), ρ is the soil density (M/L), DiW is the dispersion tensor for the ith solutein water (L2/T), q is the water ﬂux (L/T), ra is the root solute uptake term (M·L−3·T−1), and ϕi represents
the reaction sink/source term (M·L−3·T−1). Several equilibrium and nonequilibrium models are available
to simulate solute sorption to the solid phase (Šimůnek et al., 2016), while zero‐order or ﬁrst‐order coefﬁcients are used to simulate reaction processes.
2.3.2. Root Water and Solute Uptake
Beer's equation is ﬁrst used to partition reference evapotranspiration into potential transpiration, Sp, andpotential soil evaporation, Ep (e. g., Ritchie, 1972). The Leaf Area Index (LAI) is needed to partition evaporation and transpiration ﬂuxes. Equations (6) and (7) are used to calculate the mass and area of the leaves,respectively, with the latter then converted to LAI. For a detailed explanation of evapotranspiration partitioning in HYDRUS, please refer to Sutanto et al. (2012).
The actual transpiration stream, QTP, is calculated as the integral of actual root water uptake in the root zone
domain, R:
QTP ¼ ∫Sðh; x; y; zÞdxdydz ¼ SP ∫aðh; x; y; zÞbðx; y; zÞdxdydz:
R

(27)

R

where S(h,x,y,z) is actual root water uptake, a(h,x,y,z) is a dimensionless water stress response function, andb(x,y,z) is a root density distribution function. The stress response function a depends on the soil pressure
head h and has a range of values between 0 and 1. It parametrizes the effects of the energy status of soil water
caused by capillary, gravitational, and osmotic forces, as well as the oxygen deﬁcit. Feddes et al. (1978) proposed a water stress response function in which water uptake is assumed to be zero close to soil saturation(P0) and for pressure heads higher (in absolute values) than the wilting point (P3). Water uptake is assumed
to be optimal between two speciﬁc pressure heads (POpt, P2H, or P2L), which depend on the plant type. At
high potential transpiration rates (5 mm/day in the model simulation), stomata start to close at lower pressure heads (P2H; in absolute value) than at low potential transpiration rates (1 mm/day; P2L). Parameters of
the stress response function for a majority of agricultural crops can be found in various databases (e.g.,
Taylor & Ashcroft, 1972; Wesseling et al., 1991).
BRUNETTI ET AL.

8975

Water Resources Research

10.1029/2019WR025432

The root density function, b(x,y,z), modulates actual root water uptake depending on the root growth anddensity. First, equation (6) is used to calculate the root growth, which is synchronized between HYDRUS
and DPU models. Then, the equation proposed by Hoffman and van Genuchten (1983) is used to determine,at each time step, the root density and actual root water uptake at different soil locations. This numericalapproach is applied to one‐dimensional as well as multidimensional soil‐roots modeling scenarios. No feedback mechanism between the plant and the soil domain is considered. In case of soil heterogeneities, the possibility to have a spatial description of soil processes increases the accuracy of the soil‐plant model predictions(Huber et al., 2014) compared to more simplistic approaches that assume a homogeneous soil proﬁle.
In the case of passive root solute uptake, ra is equal to the product of the sink term S in the water ﬂow equation (25) and the concentration of the sink term Ci (M/L3). In the case of active root solute uptake, ra is calculated using the plant nutrient demand (which is plant and time speciﬁc) and the solute distribution in thesoil (Šimůnek & Hopmans, 2009). All compounds dissolved in water are taken up by plant roots when Ci islower than a maximum allowed solution concentration Cmax (M/L3). Details about the implementation ofroot solute uptake in HYDRUS can be found in Šimůnek and Hopmans (2009).
2.3.3. Coupling Strategy and Implementation
Neutral compounds are generally passively taken up by plants with the transpiration ﬂux. Under suchconditions, the plant uptake process is considered purely advective. As reported by Rein et al. (2011),
the error of approximating solute uptake by only advection is small and not relevant in practice. This ismathematically expressed as a one‐way coupled system, in which the soil and plant domains are sequentially coupled and no feedback mechanism exists between the two domains. Brieﬂy, HYDRUS calculates
actual root water and solute uptake, which are then used as input for the roots compartment where the
DPU starts.
The standard HYDRUS code has been modiﬁed to include the DPU component. In particular, all subroutines related to the DPU model have been collected in a separate FORTRAN ﬁle. This ﬁle contains reading
and writing subroutines, as well as the solver for equation (4), which are called by the main code at each time
step during model execution. It must be emphasized that the following coupling approach can be applied toone‐dimensional soil proﬁles (i.e., HYDRUS‐1D) as well as to multidimensional soil domains (i.e.,
HYDRUS‐2D/3D). The practical implementation can be summarized in a few steps.
1. Initialization: The HYDRUS model is set up in the usual way by deﬁning all simulated processes, domain
properties, and input parameters. At the same time, an American Standard Code for Information
Interchange ﬁle is used to collect plant physical characteristics (equation (1)), compounds' chemical
properties (equation (2)), and metabolization matrices for each compartment (equation (24)). In this step,
multiple HYDRUS and DPU output ﬁles are created. In particular, N (i.e., a number of compounds) ﬁlesare initialized for the DPU model, each of which is used to store time series of simulated compounds' concentrations in each compartment.
2. Coupled model execution: The HYDRUS main code initializes variables, reads the input, and starts thecalculation. At each time step, water ﬂow and heat and solute transport are solved ﬁrst. Simulated actualroot water and solute uptakes are then area‐averaged to unit area and passed to the DPU solver, which inturn calculates the compounds' distribution in plant's tissues. The time step is dictated by the main
HYDRUS solver. Preliminary numerical tests were carried to assess this coupling approach and investigate the stiffness of the sink term in equation (25). Such tests considered compounds with different KOW
and KAW under variable atmospheric conditions. Simulations remained numerically stable in all tests,
thus, suggesting a relatively low stiffness of the DPU model.
3. Output: If the simulation is convergent, preinitialized output ﬁles are populated with results of thenumerical simulation, which continues until the simulation period ends. Finally, mass‐balance errors
are calculated for both soil and plant domains.
2.4. Theoretical Validation
Before the proper experimental validation of the coupled soil‐plant model, the DPU component has beenvalidated against a Microsoft Excel® implementation of the DPU model proposed by Trapp (2007), whichis available on his website (https://homepage.env.dtu.dk/stt/). Theoretical validation is a needed step when
BRUNETTI ET AL.

8976

Water Resources Research
Table 1
Input parameters used in the theoretical validation of the dynamic
plant uptake modelParameter3

W (cm /gfw)
L (gl/gfw)
M (gfw)2
SA (cm /gfw)−1
τ11 (day )Air
3
CA (gs/cm )fp (−)
vdep (cm/day)
Compoundlog KOW (−)
33
KAW (cm /cm )Roots

Shoots

0. 89
0.025
1000
0.0
0.15

0. 85
0.02
1000
50.0
0.15

1e‐13
0.5
8640
2.00
0.00005

10.1029/2019WR025432

testing and developing models since it avoids the inﬂuence of unpredictable sources of uncertainty of experimental data, which could mask inaccuracies in the model implementation.
2.4.1. Synthetic Scenario Description
The synthetic scenario considers the translocation of a contaminant in aperennial shrub conceptualized in two compartments (i.e., roots andshoots). The plant is exposed to contaminated soil water and air for asimulation period of 1,200 days. Two 1‐day injection pulses of 0.01 g ofthe compound are used to load the system at Days 300 and 600. Gaseous
uptake, particle deposition, and compound's volatilization are consideredin the numerical simulation, while plant growth is neglected (i.e., a constant plant mass, M). Constant actual transpiration of 1,000 cm3/day is
assumed. The plant's permeability PP is reduced to a scalar value, whichis assumed to be constant and equal to 8,640 cm/day. The time series ofsimulated solute concentrations in different compartments is used to validate the DPU model. Plant's characteristics and input parameters are summarized in Table 1.

2.5. Experimental Validation
Measured data are taken from a recently published study by Kodešová et al. (2019b), which focused on thetranslocation and biodegradation of selected pharmaceuticals in different plants. In the present study, only
measured data about dynamic uptake of CBZ and its main metabolites (i.e., 10,11‐epoxide carbamazepine[EPX], and oxcarbazepine [OXC] in lamb's lettuce, spinach, and arugula) are used. The modeling study isfurther enriched by a global sensitivity analysis to identify the most important factors inﬂuencing the accumulation of CBZ in spinach shoots. The overall analysis aims to provide a preliminary assessment of thedeveloped model against observed data and to demonstrate its ﬂexibility and use potential. In the following
sections, we provide a general overview of the experiment, modeling settings, and numerical methods used.
For a thorough description of the experiment, please refer to Kodešová et al. (2019b).
2.5.1. Case Study Description
Experiments were carried out under greenhouse conditions, with average air temperature and humidityof 22 °C and 35%, respectively. A single plant was planted in a small pot in ﬁve replicates of each soil
(three soils were used) and treatment, which consisted of an injection of a single compound or a mixtureof three pharmaceuticals. In the present study, only data relative to an injection of a single compound in
Haplic Chernozem soils are considered. The irrigation regime and injected concentrations are reported inTable 2. Each pot was perforated at the bottom and had a height of
5.5 cm and a surface area of 25 cm2 and, additionally, was placed on
Table 2
top of a plastic plate that provided a limited accumulation of water durIrrigation doses for a single plant and concentrations of carbamazepine
ing the experiment. The soil surface was covered with a thin layer of3
Irrigation (cm )quartz‐sand gravel to avoid cross‐contamination of treated soil pots.
Carbamazepine
3
An evaporation pan placed near the pots measured evaporation in theDay
Spinach/Arugula
Lamb's lettuce[g/cm ]
greenhouse, which was generally constant and approximately equal to0
Sprout0.2 cm/day. A pan coefﬁcient of 1.0 was assumed to estimate the refer16
30.0
0.0
6.50e‐07
ence evapotranspiration (Zeng et al., 2009).
18
40.0
0.0
6.80e‐07
21
23
25
26
28
32
33
37
39
40
42
43

30.0
30.0
30.0
33.3
30.0
30.0
33.3
30.0
30.0
0.0
Harvest

BRUNETTI ET AL.

0.0
30.0
40.0
33.3
30.0
30.0
33.3
30.0
30.0
30.0
0
Harvest6.80e‐07
6.40e‐07
6.40e‐07
6.40e‐07
6.60e‐07
6.60e‐07
6.60e‐07
6.60e‐07
6.60e‐07
6.60e‐07

Plants were exposed to natural light. Seeds of arugula, lamb's lettuce, andspinach germinated in a horticultural substrate and were then moved in
the soil pot after sprouting. Each plant was initially irrigated with fresh
water and then exposed to the injection of contaminated water. Sets of ﬁve
plants of each soil were also irrigated with solute‐free water to test the possible impact of the compounds on plant growth. The measured CBZ inﬂow
concentration varied between 6.4e‐07 and 6.8e‐07 g/cm3 (Table 2). After
exposure, plants were removed from the soil, rinsed, and divided into roots
and shoots for the next chemical analysis, which was carried out using theliquid chromatography‐tandem mass spectrometry. Before the sample

8977

Water Resources Research
Table 3
Parameters and their bounds used in the experimental validation andglobal sensitivity analysis of the coupled HYDRUS‐dynamic
plant uptake model
Lamb'slettuceaParameter33θr (cm /cm )33
θs (cm /cm )−1α (cm )n (−)
Ks (cm/day)3
ρb (g/cm )λL (cm)KfCBZ
3β1−β −1(cm ·μg ·g )βCBZ (−)−1μLCBZ (day )−1μSCBZ (day )P0 (cm)POpt (cm)P2H (cm)
P2L (cm)
P3 (cm)Roots
3
W (cm /gfw)
L (gl/gfw)
max
(gfw)
M
0
M (gfw)
gr−1
K (day )−1
τ11 (day )−1
τ21′ (day )−1
τ31′ (day )−1
τ22 (day )Shoots3
W (cm /gfw)
L (gl/gfw)
max
(gfw)
M
0
M (gfw)
gr−1
K (day )−1
τ11 (day )−1
τ21′ (day )−1
τ31′ (day )−1
τ22 (day )2
Aharvest (cm )LAIharvest (−)
Compounds
log KOW(CBZ) (−)
log KOW(EPX) (−)
log KOW(OXC) (−)
m(CBZ) (g/mol)m(EPX) (g/mol)m(EPX) (g/mol)aSpinachArugulaSensitivity
analysis

0.08
0.39
0.05
1.22
52
1.09
2
2.97

(0.06, 0.10)(0.30, 0.50)(0.01, 0.1)
(1.1, 2.0)
(20, 100)(1, 1.3)(1, 4)(2, 5)

−3
−20
−200
−800
−8000

0.88
0.00470.0047
−3
−45
−200
−800
−8000

−3
−45
−200
−800
−8000

(0.7, 1.0)
(0.001, 0.01)(0.001, 0.01)(−10, 0)(−35, −11)(−300, −100)
‐
‐

0.89
0.015
126
4
0.15
0.15
0.11
0.026
0.004

0.85
0.01
206
10
0.165
0.165
0.133
0.026
0.004

0.85
0.01
149
4
0.15
0.026
0.023
0.00260.004(0.8, 0.9)
(0.005, 0.025)(100, 300)
(1, 20)(0.08, 0.26)(0.05, 0.3)
‐
‐
‐

0.89
0.015
440
18
0.15
0.45
0.37
0.07
0.004
47.5
1.9

0.85
0.01
819
50
0.165
0.355
0.315
0.035
0.004
65.6
2.6

0.86
0.01
723
60
0.15
0.026
0.023
0.0026
0.004
48.9
2.0

(0.8, 0.9)
(0.005, 0.025)(500, 1,000)(10, 100)(0.08, 0.26)(0.05, 0.3)
−
−
−
−
−

2.25

(2.0, 2.8)

1.26

‐

1.11

‐

236.27252.28252.28‐‐‐θs and θr are the saturated and residual water contents, respectively,
α and n are the empirical van Genuchten‐Mualem shape parameters, Ks
is the saturated hydraulic conductivity, and P0, POpt, P2H, P2L, and P3
are Feddes' parameters.
CBZ = carbamazepine, EPX = 10,11‐epoxide, LAI = Leaf Area Index.

BRUNETTI ET AL.

10.1029/2019WR025432

extraction, the fresh and dry weights of roots and shoots, as well as theone‐sided leaves area, were determined.
2.5.2. Model Setup
The one‐dimensional model HYDRUS‐1D is used to simulate transport
processes in the soil. The 5.5‐cm deep soil proﬁle is discretized in 100 ﬁnite
elements. An atmospheric boundary condition is applied at the soil surface using (a) irrigation and potential evapotranspiration ﬂuxes, (b) a prescribed zero pressure head (i.e., full saturation) during ponding, and (c)equilibrium between the soil surface pressure head and the atmosphericwater vapor pressure when the atmospheric evaporative demand cannot
be met. A seepage face boundary condition is set at the bottom of thepot. This condition assumes that the boundary ﬂux remains zero as longas the pressure head is below a certain threshold value, hseep. However,when the lower end of the soil proﬁle reaches this value, hseep is imposedat the lower boundary and outﬂow is calculated accordingly. In the present study, hseep is set to 1 cm to simulate the accumulation of water inthe plate under the pot. The concentration ﬂux across the top boundary
is simulated using the classic Cauchy‐type boundary condition, while azero concentration gradient is imposed at the bottom. The initial pressure
head is assumed to be constant and equal to −100 cm in the entire model
domain. The soil is considered solute free at the beginning of
the simulation.
Parameters used in the experimental validation are listed in Table 3. The
unimodal van Genuchten‐Mualem model (van Genuchten, 1980) is used
to describe the soil hydraulic properties, which are obtained from the measured particle size distribution and soil bulk density ρb (grams per centimeter) using pedotransfer functions from ROSETTA (Schaap et al.,2001). The piecewise linear Feddes function describes the root water stress(Feddes et al., 1978) and regulates actual transpiration. Feddes' parameters (i.e., P0, POpt, P2H, P2L, and P3 in Table 3) are taken from the literature (Wesseling et al., 1991) and slightly adjusted to better reproduce
measured data. Reference evapotranspiration is partitioned into potential
evaporation and transpiration using the LAI. Growth rates are in line withvalues reported by Glenn et al. (1984) and Gent (2002). Final measured
roots and shoots wet and dry masses, as well as one‐sided leaves area—
which averaged to an area of 1 m2—are used to constrain equation (6).
The advection‐dispersion‐reaction equation (equation (26)) is used tosimulate the transport of CBZ in the soil. The effect of molecular diffusion
is neglected. The solute dispersivity λL is set to 2 cm. In their study,
Kodešová et al. (2015) reported nonlinear equilibrium adsorption of CBZ
on Haplic Chernozem, which can be described using the Freundlich
adsorption isothermsi ¼ K f Cβi ;(28)

where Kf and β are empirical coefﬁcients. Measured data show that the
EPX is the only CBZ metabolite detected in the soil during the experiment.
Thus, the CBZ degradation process in the soil is simulated using a sequential decay chain with ﬁrst‐order rate coefﬁcients, μLCBZ and μSCBZ (T−1),
that account for the transformation of CBZ in EPX in the liquid and solid
phases, respectively. Degradation of CBZ in this soil was studied by
Kodešová et al. (2016). In this study, the degradation half‐life was larger
than 1,000 d, that is, μLCBZ and μSCBZ were lower than 0.0007 day−1.
8978

Water Resources Research

10.1029/2019WR025432

Figure 3. A schematic of the (a) theoretical degradation pathway of carbamazepine in plants, (b) the dynamic plant
uptake implementation, (c) and the metabolization matrix for the jth compartment. The black circle indicates the compound degradation in by‐products, which are not of interest in the analysis (i.e., reaction yield).

However, ﬁnal concentrations of EPX in soils in the case of the pot greenhouse experiment (Kodešová et al.,2016) indicated higher degradation rates. Therefore a new degradation experiment was carried out withfresh soil, which was preincubated with water (20 °C) before application of pharmaceuticals. Resulting
half‐life was 0.0047 day−1. The same sorption parameters are assumed for EPX, whose further
degradation in the soil is neglected (Koba et al., 2016).
Plants are conceptualized in two compartments: roots and shoots. No CBZ was present in the greenhouseduring the experiment and thus gaseous uptake and particle deposition are neglected. Similarly, volatilization is excluded from modeling due to the nonvolatile behavior of CBZ. Octanol‐water partition coefﬁcientsare obtained using the U.S. Environmental Protection Agency Estimation Programs Interface suite program(U.S. EPA, 2013). Water contents of plants' compartments are estimated from the measured fresh and dryweights of roots and shoots and slightly adjusted so that modeled and measured concentrations better
match. The lipid contents are set based on literature values reported for leafy vegetables (Legind et al.,2011; Legind et al., 2012; Trapp & Eggen, 2013; Wen et al., 2016) and slightly increased to compensate theunderestimation of the sorption introduced by equation (8). Contrary to the soil, multiple CBZ metabolites
were observed in plants. In particular, EPX and OXC were the most signiﬁcant by‐products while only small
concentrations of trans‐10,11‐dihydro‐10,11‐dihydroxy carbamazepine (RTC) and 10,11‐dihydro carbamazepine (DHC) were observed. In the present study, we focus only on the transformation of CBZ in its main
active metabolites, EPX and OXC. RTC and DHC are not of immediate interest in the analysis; and thus, they
are not explicitly included in the DPU metabolization chain but instead considered in the reaction yield of
EPX and CBZ. Theoretical degradation pathways of CBZ in plants and its conceptual implementation in
the DPU model are shown in Figure 3. Input concentrations, as well as sorption parameters, are converted
to molar units before executing the model, while output concentrations are reported in nanogramg per gramfor visualization purposes.
2.5.3. Global Sensitivity Analysis
Accurate integrated models can be used to understand how different components of the simulated processinteract and to identify the key driving factors. A classic example is to perform a sensitivity analysis on themodel to investigate the impact of different parameters on some variables of interest. The sensitivity analysis
is frequently used in modeling studies under different settings and aims (e.g., Brunetti et al., 2016; Brunetti
et al., 2017; Brunetti, Šimůnek, & Bautista, 2018; Brunetti, Šimůnek, Turco, & Piro, 2018). In the present
study, we apply a global sensitivity analysis to identify the most important factors among soil, plant, and chemical parameters that drive the accumulation of CBZ in spinach shoots. To this aim, the developed soil‐plantmodel is coupled with a modiﬁed version of the Morris method (Morris, 1991) proposed by Campolongo
et al. (2007).
BRUNETTI ET AL.

8979

Water Resources Research

10.1029/2019WR025432

The Morris method (Morris, 1991) belongs to the class of screening methods. Screening methods aim to provide qualitative sensitivity measuresfor different factors using a relatively small number of model evaluations.
In general, the Morris method is a one‐factor‐at‐a‐time local method since
it computes the elementary effect by changing only one factor at a time.
However, it can be viewed as a global method since it averages several elementary effects computed at different points in the parameter space, providing a qualitative measure of the importance and interactions ofdifferent factors. The main outcomes of the analysis are two sensitivity
indices, σ and μ*, which summarize the interaction and overall importance of the parameter, respectively. High values of these sensitivity
indices indicate a signiﬁcant inﬂuence of the parameter on the variableof interest. To interpret the results by simultaneously taking into accountboth sensitivity measures, Morris suggested their graphical representationin the (μ*−σ) plane.
The same model settings used in the experimental validation for spinach
are used in the sensitivity analysis. Since the main aim is to identify factors
driving the accumulation of CBZ in shoots, its metabolites are neglected,
Figure 4. A comparison between the simulated solute concentrations in and the metabolization matrix is reduced to a single scalar value for eachroots (black) and shoots (gray) of the Trapp (dashed lines) and dynamic
compartment. Furthermore, preliminary numerical simulations haveplant uptake (DPU; solid lines) models, respectively.
shown that during the experiment, the soil was in nearly saturated conditions; thus, only the “anaerobiosis” and “optimal” pressure head values in
the Feddes' model, P0 and POpt, are included in the sensitivity analysis, which reduces the number ofscreened parameters to 27. Table 3 reports these screen parameters with their ranges of variation. A sample
size of 200 is used to generate the initial sample, leading to a total of 5,400 numerical simulations. The analysis is carried out using the Python programming language. For each sample, the code overwrites the input,
executes the coupled HYDRUS‐DPU model, and stores the value of the ﬁnal CBZ concentration in shoots ina preassigned array. Finally, the sensitivity indices are calculated, and the μ*−σ plot is used to target
important factors.

3. Results
3.1. Theoretical Validation
Results of the theoretical validation are reported in Figure 4, which shows a comparison between the simulated solute concentrations in roots (black) and shoots (gray) for the Trapp (dashed lines) and DPU (solidlines) models, respectively. Simulated concentrations overlap, thus conﬁrming the correctness of the DPU
model implementation. The implemented model can accurately reproduce both the solute accumulationin roots after the pulse injection and the following decrease induced by the metabolization. The low compound translocation in shoots due to the relatively high solute lipophilicity (i.e., log KOW = 2.00) is also correctly described by the model. The implicit time discretization method used to solve the mass‐balance
equation (equation (4)) is responsible for the slight smoothing of results observed in the roots compartment(Figure 4) and explains the negligible deviation between the Trapp and DPU models. Nevertheless, it can beconcluded that the DPU model implementation is theoretically correct, numerically stable, and can be used
for further analyses.
3.2. Experimental Validation
The main outcomes of the experimental validation are reported in Figure 5, which shows a comparison
between the simulated and measured concentrations of CBZ and its main metabolites in the soil, roots,
and shoots of lamb's lettuce, spinach, and arugula. Overall, the coupled HYDRUS‐DPU model satisfactorily
reproduces the ﬁnal accumulation of different compounds in the soil and plants' tissues. The highest deviation is observed in the lamb's lettuce scenario, for which the model underestimates the high enzymatic biodegradation observed in shoots (Kodešová et al., 2019b). Such effect, which is more pronounced for the EPX,
can be mainly related to some bias in the assumed ﬁrst‐order metabolization coefﬁcients, as well as to plausible uncertainties in soil and plant input parameters. Nevertheless, simulated compound concentrations arein the measurement uncertainty ranges, suggesting an overall good performance of the model. This is further
BRUNETTI ET AL.

8980

Water Resources Research

10.1029/2019WR025432

Figure 5. A comparison between the simulated (solid lines) and measured (dots) concentrations of carbamazepine (left) and its main metabolites 10,11‐epoxide(middle) and oxcarbazepine (right) in the soil (red), roots (gray) and shoots (black) of lamb' lettuce (top), spinach (middle), and arugula (bottom). The blue bars
are the solute irrigation ﬂuxes (Table 2). The error bar indicates the estimated 30% uncertainty of the measurement method, according to Kodešová et al. (2019b).
The reaction yield indicates the amount of metabolized carbamazepine converted in 10,11‐epoxide and oxcarbazepine.

BRUNETTI ET AL.

8981

Water Resources Research

10.1029/2019WR025432

Figure 6. (a) Simulated actual transpiration (black line) and irrigation volumes (blue bars); (b) simulated roots (gray line) and shoots (black line) masses; (c) the
average root zone pressure head (gray line) and root solute uptake (red line) in the spinach modeling scenario. P0 and POpt are the “anaerobiosis” and “optimal”
pressure head values in the Feddes' model.

conﬁrmed by the higher ﬁtting accuracy achieved in the spinach and arugula modeling scenarios. Here, themodel accurately reproduces both the biodegradation pathways of CBZ in roots and shoots while someoverestimation is evident in the soil. Such behavior is more pronounced for arugula, for which the model
overpredicts both the CBZ and EPX concentrations, thus suggesting that the bias mainly stems frompotential inaccuracies in the description of soil transport processes. On the other hand, the compounds'
bioaccumulation in spinach is well described in all model compartments, though the CBZ concentrationin the soil is again slightly overestimated.
One of the main advantages of the proposed coupled model is its ability to provide a dynamic description oftransport and degradation processes in the soil‐plant continuum, which is crucial to developing a better
understanding of physical processes happening during the experiment. When irrigation begins in the experiment, DPU starts and CBZ is translocated with transpiration streams to the shoots. At the same time, it isbiodegraded by plant's enzymes in its by‐products, which are compartmented in leaves and eventuallyfurther metabolized. The ﬂuctuations in the simulated concentrations are mainly due to the actual transpiration patterns induced by the variably saturated conditions in the soil. This is clearly shown in Figure 6,
which shows the simulated root zone pressure head and actual root solute uptake for the spinach modelingscenario. Model results suggest that irrigation caused nearly saturated conditions in the pot, which induced
anaerobic stress (as deﬁned by Feddes' stress response function) on plant's roots, thereby limiting water andsolute uptake (Figure 6). This is particularly evident in spinach roots and shoots between Days 20 and 30 of
the simulated period when the CBZ concentration declines (Figure 5) due to enzymatic biodegradation andsimultaneous negligible supply of contaminant with transpiration water. After this period, the pressure headdecreases below the anaerobiosis pressure head value and DPU starts again with a subsequent translocationand compartmentation of EPX and OXC in leaves.
Despite the promising results, it must be emphasized that the amount and information content of measured
data available for validating the model was rather limited. For this reason, we avoid numerical parameter
optimization, which could improve the quality of the ﬁt but would likely lead to unphysical parametervalues and high predictive uncertainty. In this perspective, future research will focus on proper Bayesianmodel validation against highly informative measurements from selected experiments, which are already
BRUNETTI ET AL.

8982

Water Resources Research

10.1029/2019WR025432

ongoing at the Czech University of Life Sciences, Prague. Nevertheless,this preliminary experimental validation demonstrates that the developedmodel shows good promise for describing transport processes in the soil‐
plant continuum.
3.2.1. Global Sensitivity Analysis
The results of the experimental validation suggest an important role of thesoil in the compounds' translocation process in plants. To further investigate this aspect, a global sensitivity analysis is performed to screen factors
inﬂuencing the accumulation of CBZ in shoots. Results of the analysis,
carried out using the Morris method, are summarized in the scatter plot
μ*−σ shown in Figure 7.
Three parameter groups can be identiﬁed from the distribution of sensitivity measures. The ﬁrst group, in order of importance, includes the shoots
growth rate Kgrshoots, the van Genuchten‐Mualem shape parameters, αand n, and the ﬁrst‐order degradation rate in the shoots, τshoots11. The
effect of the latter is expected since it directly regulates simulated biodegradation of CBZ in the shoots, thus modulating the sink term in equation (4). The outcome for α and n holds more interest, further
conﬁrming ﬁndings of the experimental validation on the role of soil
Figure 7. Scatter plots of the Morris sensitivity measures for various soil
hydraulic behavior on solute translocation in plants during the experihydraulic, plant, and solute transport parameters.
ment. Parameters α and n describe the start and magnitude of the desaturation process in the soil‐water retention curve, thus governing thesimulated soil retention capacity. In particular, under nearly saturated conditions, small perturbations ofα, generated by the Morris sampling process, propagate in the simulated root zone pressure head, which
now oscillates near the anaerobiosis point deﬁned in the Feddes' model, thus strongly inﬂuencing actual root
water and solute uptake. Such behavior is exacerbated in combination with simultaneous changes of n, as
demonstrated by the high interaction effect, σ, calculated for these parameters. This suggests that an accurate estimation of these parameters would signiﬁcantly increase the model predictive capabilities and reduce
the overall uncertainty. Similarly, the shoots growth rate, Kgrshoots, exhibits a high interaction potential, as it
is mainly involved in the compound's dilution effect in the plant. In this case, few measurements of the plantmass during its growth stage would be sufﬁcient to characterize the plant's growth and thus to better constrain the model.
The second group includes the initial and ﬁnal shoots masses, M0shoots and Mmaxshoots, the saturated watercontent and hydraulic conductivity, θs and Ks, respectively, and the anaerobiosis and optimal pressure headvalues in the Feddes' model, P0 and POpt, respectively. As expected, the Feddes' parameters have a signiﬁcant impact on CBZ bioaccumulation in the shoots since they regulate simulated actual root water uptake,and thus transpiration streams, in nearly saturated conditions. In such circumstances, θs and Ks strongly
inﬂuence the soil hydraulic behavior, which in turn affects the root uptake process. Finally, the dilution
effect mainly explains the moderate sensitivity of the simulated CBZ accumulation to M0shootsand Mmaxshoots.
The third group includes all remaining parameters, which for physical reasons have proven to have a limited
effect on the variable of interest. For instance, the residual water content θr has a negligible impact on thesoil hydraulic behavior in wet conditions, which were observed during the experiment. Similarly, in the analyzed modeling scenario, the octanol‐water partition coefﬁcient log KOW has a limited inﬂuence since its
range of variation does not signiﬁcantly modify lipophilicity of the compound. While it is likely that wider
bounds could lead to a different sensitivity ranking, its scientiﬁc meaning could be highly questionable
due to the adoption of an unphysical value of log KOW for CBZ, whose value is known and reported in theliterature. Interestingly, solute transport and adsorption parameters exhibit low sensitivity measures. Suchbehavior can be partially explained by the limited volume of the pot, as well as by the overall
experimental conditions.
It should be emphasized that the present sensitivity analysis is aimed at providing insight into the mostimportant factors affecting DPU of CBZ during a speciﬁc experiment. Results should thus not be regarded
BRUNETTI ET AL.

8983

Water Resources Research

10.1029/2019WR025432

as general but should instead be restricted to a particular investigated modeling scenario. Different model
settings can certainly lead to different results. However, the present analysis demonstrates how the developed soil‐plant model can be used in conjunction with advanced statistical techniques to better understandphysical processes.

4. Discussion
4.1. Integrated Soil‐Plant Model
The analysis demonstrates how the implementation of the DPU module in the HYDRUS model signiﬁcantly
extends its capabilities by providing an integrated modeling environment for contaminant risk assessment in
the soil‐plant continuum. The adopted coupling strategy is numerically stable, computationally effective,
and mathematically accurate, as demonstrated by the theoretical and experimental validation.
Furthermore, the comprehensive description of the contaminant transport in the soil and plant can be usedin conjunction with advanced statistical tools, such as global sensitivity and Bayesian analyses, to gain
insights into the investigated process and to better explain interactions between the soil and plant domains.
Additionally, the modeling approach—based on multiple metabolization matrices—increases modelingﬂexibility. To continue these advances, further layers of complexity can be added to the model by reﬁning
the description of enzymatic and bacterial biodegradation in the plant's tissues using the Michaelis‐Menten and Monod kinetics, respectively. The former has been observed for the enzymatic degradation of
cyanide in plants (Larsen et al., 2004; Yu et al., 2004) and the exclusion of salt NaCl and NaF from roots(Clausen et al., 2015; Trapp et al., 2008). Future developments could also include an option to use an optimalplant growth function (equation (6)) and reduce the growth due to various stresses such as saturation, osmotic, or chemical stresses. A comparison with more simplistic numerical approaches (e.g., tipping bucket
models) is recommended to assess, under what conditions, this increase in the soil‐plant model complexity
is justiﬁable.
One of the main limitations of the proposed model is the simpliﬁcation of the plant physiology compared toother soil‐plant models (e.g., Fatichi et al., 2019, 2016; Huang et al., 2018, 2017; Huber et al., 2014; Manoli
et al., 2017; Werner & Dubbert, 2016), some of which provide a detailed multidimensional mechanistic
description of the xylem and phloem ﬂow, photosynthesis, carbon allocation, leaf‐level gas exchange, energybalance, etc. However, this increase in the model complexity requires highly informative measurements andposes computational cost issues, especially for model calibration, that can undermine its usability for practical applications at the plot scale. Such problems are mitigated in the model developed in this study, which
can be used in various applications, such as phytoremediation involving uptake and metabolization processes of toxic compounds from contaminated soils (e.g., Ouyang, 2002), pesticides fate assessment (e.g.,Legind et al., 2011), remediation of contaminated groundwater plumes (e.g., Hong et al., 2001), assessmentof salinity in agricultural arid areas (e.g., Trapp et al., 2008), and air pollution removal by urban vegetation(e.g, Yang et al., 2008). To further extend the range of the model applicability, future research should focus
on the modiﬁcation of the model structure to consider DPU of electrolytes and reﬁne the description of theroot uptake mechanism of large molecules (Kumar & Gupta, 2016).
Nevertheless, this study cannot cover all potential modeling scenarios, and thus, additional research isneeded to assess the robustness of the proposed approach under various operating conditions. In particular,
the role of numerical stiffness on the choice of the time step must be clariﬁed, which in turn affects the stability of the overall scheme. In the proposed coupling strategy, the HYDRUS internal solver dictates the time
step for the DPU module. The time step is selected in HYDRUS based on convergence criteria of the numerical solution of the Richards equation and stability criteria of the numerical solution of the convection‐
dispersion equation. Thus, the choice of the time step during the simulation mainly depends on the
dynamics of soil transport processes, which are assumed to be happening at a ﬁner time scale comparedto the plant translocation. Such an assumption is expected to hold in the majority of modeling scenarios.
However, potential numerical issues could arise when simulating DPU of highly volatile compounds fromsoils in nearly steady‐state conditions. Under such circumstances, the time step chosen by HYDRUS could
not be sufﬁciently small so as to catch the dynamics of the volatilization process in the plant. Thus, werecommend future studies to focus on this aspect to target eventual instability problems, as well as on a better description of nonequilibrium sorption and biodegradation in plants.
BRUNETTI ET AL.

8984

Water Resources Research

10.1029/2019WR025432

4.2. The Key Role of Soil Transport Processes
The results of the experimental validation and global sensitivity analysis emphasize the role of soil transport
processes on the compound's bioaccumulation in plants. Such an effect is mainly explained by the inﬂuenceof unsaturated conditions on actual transpiration streams, which govern the contaminants' translocation inplants. Existing DPU models greatly simplify the unsaturated hydraulic behavior of the root zone using questionable coarse theoretical conceptualizations. A classic example is a use of tipping bucket models (Legind
et al., 2012; Trapp & Eggen, 2013), which suffer from multiple theoretical and practical limitations(Emerman, 1995; Scanlon, Christman, et al., 2003). Our implementation, based on the HYDRUS model,
overcomes such limitations, providing a modeler the possibility to mechanistically simulate single porosityand preferential ﬂows, as well as equilibrium and nonequilibrium solute transport processes. The results ofthis study demonstrate that such modeling features are crucial to understanding how the interaction
between variably saturated ﬂow conditions and the plant's respiration affect dynamic root water and solute
uptake (Figure 6).
4.3. Model Complexity and Uncertainty Assessment
The high modeling ﬂexibility of the developed soil‐plant model leads to high model parameterization.
Several parameters are needed to describe physical and chemical processes in the soil and plant.
Furthermore, the use of multiple metabolization matrices for the simulation of differentiated biodegradationin plant's tissues introduces additional layers of model complexity, which are certainly a valuable tool for themodeler when the main aim of the analysis is the model‐based design (e.g., Lange et al., 2015). Under suchcircumstances, parameters are established from the literature, and the model is used to design the systemand simulate its future behavior under actual conditions.
The situation is radically different when the main goal of the analysis is predictive modeling. Under suchcircumstances, parameters must be inferred from measured data through model calibration. However,highly parameterized models often suffer from poor predictive capabilities due to high parameter uncertainty, mainly stemming from the limited information content of measured data. The situation can befurther worsened when simulating highly nonlinear processes, as in the DPU model. For these reasons,the information content of measured data is critical to strengthening the generalizability of the coupledsoil‐plant model, which must be assessed through a Bayesian calibration framework. In the present study,
we decided to avoid a perfunctory model calibration due to the limited amount of available measured data,
which would lead to poor model generalizability. Therefore, we recommend future research to focus on amore rigorous uncertainty assessment using more comprehensive measured data sets to clarify how the highparameterization of the soil‐plant model affects its predictive capabilities in real case studies. A soil‐plantcolumn experiment under variable atmospheric conditions (e.g., a climate chamber) and different soluteinjections that would include measurements of pressure heads and soil water contents at multiple depths,
as well as solute concentrations in plants, would be of great interest. The information content of the experimental data set could be increased by targeted laboratory measurements, such as unsaturated hydraulicproperties (e.g., Peters & Durner, 2008) and solute sorption on plant's tissues (e.g., Li et al., 2005), which
could then be used as prior information for subsequent inverse parameter estimation using transient data.
We advocate that such an approach should be ubiquitous when using numerical simulations of DPU of
selected compounds. Thus far, research in this ﬁeld has focused on the use of conceptual models in combination with local sensitivity analyses, which—despite good intentions—suffered from theoretical(Emerman, 1995; Scanlon, Christman, et al., 2003) and statistical limitations (Saltelli & Annoni, 2010) to
capture the main physical aspects of the investigated process.

5. Conclusions
The main aim of the present study was to develop a comprehensive model capable of providing a mechanistic description of transport and reaction processes in the soil‐plant continuum. The results of the multilevelnumerical analysis demonstrate that the developed model shows good promise for reproducing physicochemical processes involved in DPU of chemicals from contaminated soils. The possibility of having an integrated mechanistic description of both soil and plant behaviors ﬁlls a signiﬁcant scientiﬁc gap in this ﬁeld,
often characterized by fragmented research experiences. For example, the developed model can be used inenvironmental applications to simulate uptake of chemicals from contaminated sites and their
BRUNETTI ET AL.

8985

Water Resources Research

10.1029/2019WR025432

translocation and transformation in the soil‐plant continuum in a comprehensive contaminant risk assessment. We certainly do not expect this study to be exhaustive in terms of exploring all plausible modeling scenarios nor conclusive in the highly dynamic scientiﬁc ﬁeld of soil‐plant modeling. The numerical descriptionof plant physiology can certainly be improved, and the model stability, usability, and generalizability should
be further explored. The latter is of especially great importance since the high parameterization of the developed model can result in poor predictive capabilities. As we pointed out in the section 4, further model testing and experimental validation are needed to assess this and other aspects. Nevertheless, the high modelingﬂexibility of the proposed numerical approach is undoubtedly a very useful tool for developing reliable mitigation strategies for environmental pollution problems.
Acknowledgments
Authors acknowledge the ﬁnancial
support of the Czech Science
Foundation Project 17‐08937S,
behavior of pharmaceuticals in the soil‐
water‐plant system. The work was also
supported from European Structuraland Investment Funds, Operational
Programme Research, Development
and Education and Ministry of
Education, Youth and Sports projects
supporting the development ofinternational mobility of research staff
at CULS Prague (CZ.02.2.69/0.0/0.0/
16_027/0008366), and NutRisk(CZ.02.1.01/0.0/0.0/16_019/0000845).
The data that support the ﬁndings ofthis study are available at https://zenodo.org/record/3275850#.XSSrHY_
ONPY.

BRUNETTI ET AL.

References
Barac, T., Taghavi, S., Borremans, B., Provoost, A., Oeyen, L., Colpaert, J. V., et al. (2004). Engineered endophytic bacteria improve phytoremediation of water‐soluble, volatile, organic pollutants. Nature Biotechnology. https://doi.org/10.1038/nbt960
Bertrand, J., Laffont, C. M., Mentré, F., Chenel, M., & Comets, E. (2011). Development of a complex parent‐metabolite joint population
pharmacokinetic model. AAPS Journal. https://doi.org/10.1208/s12248‐011‐9282‐9
Boote, K. J., Sau, F., Hoogenboom, G., & Jones, J. W. (2008). Experience with Water Balance, Evapotranspiration, and Predictions of Water
Stress Effects in the CROPGRO Model. In L. R. Ahuja, V. R. Reddy, S. A. Saseendran, & Q. Yu (Eds.), Response of Crops to Limited Water:
Understanding and Modeling Water Stress Effects on Plant Growth Processes, Adv. Agric. Syst. Model. 1 (pp. 59–103). Madison, WI: ASA,CSSA, SSSA. https://doi.org/10.2134/advagricsystmodel1.c3
Briggs, G. G., Bromilow, R. H., & Evans, A. A. (1982). Relationships between lipophilicity and root uptake and translocation of non‐ionised
chemicals by barley. Pesticide Science, 13(5), 495–504. https://doi.org/10.1002/ps.2780130506
Brodersen, C. R., Roddy, A. B., Wason, J. W., & McElrone, A. J. (2019). Functional status of xylem through time. Annual Review of Plant
Biology, 70(1), 407–433. https://doi.org/10.1146/annurev‐arplant‐050718‐100455
Broyden, C. G. (1969). A new method of solving nonlinear simultaneous equations. Computer Journal, 12(1), 94–99. https://doi.org/
10.1093/comjnl/12.1.94
Brunetti, G., Porti, M., & Piro, P. (2018). Multi‐level numerical and statistical analysis of the hygrothermal behavior of a non‐vegetatedgreen roof in a mediterranean climate. Applied Energy, 221, 204–219. https://doi.org/10.1016/J.APENERGY.2018.03.190
Brunetti, G., Saito, H., Saito, T., & Šimůnek, J. (2017). A computationally efﬁcient pseudo‐3D model for the numerical analysis of boreholeheat exchangers. Applied Energy, 208, 1113–1127. https://doi.org/10.1016/j.apenergy.2017.09.042
Brunetti, G., Šimůnek, J., & Bautista, E. (2018). A hybrid ﬁnite volume‐ﬁnite element model for the numerical analysis of furrow irrigationand fertigation. Computers and Electronics in Agriculture, 150, 312–327. https://doi.org/10.1016/J.COMPAG.2018.05.013
Brunetti, G., Šimůnek, J., Bogena, H., Baatz, R., Huisman, J. A., Dahlke, H., & Vereecken, H. (2019). On the information content of cosmic‐
ray neutron data in the inverse estimation of soil hydraulic properties. Vadose Zone Journal, 18(1), 0. https://doi.org/10.2136/
vzj2018.06.0123
Brunetti, G., Šimůnek, J., & Piro, P. (2016). A comprehensive numerical analysis of the hydraulic behavior of a permeable pavement.
Journal of Hydrology, 540, 1146–1161. https://doi.org/10.1016/j.jhydrol.2016.07.030
Brunetti, G., Šimůnek, J., Turco, M., & Piro, P. (2018). On the use of global sensitivity analysis for the numerical analysis of permeable
pavements. Urban Water Journal, 15(3), 269–275. https://doi.org/10.1080/1573062X.2018.1439975
Butcher, J.C., 2016. Numerical methods for ordinary differential equations, numerical methods for ordinary differential equations. https://doi.org/10.1002/9781119121534
Campolongo, F., Cariboni, J., & Saltelli, A. (2007). An effective screening design for sensitivity analysis of large models. Environmental
Modelling & Software, 22(10), 1509–1518. https://doi.org/10.1016/j.envsoft.2006.10.004
Cheyns, K., Mertens, J., Diels, J., Smolders, E., & Springael, D. (2010). Monod kinetics rather than a ﬁrst‐order degradation model explainsatrazine fate in soil mini‐columns: Implications for pesticide fate modelling. Environmental Pollution, 158(5), 1405–1411. https://doi.
org/10.1016/j.envpol.2009.12.041
Clausen, L. P. W., Karlson, U. G., & Trapp, S. (2015). Phytotoxicity of sodium ﬂuoride and uptake of ﬂuoride in willow trees. International
Journal of Phytoremediation, 17(4), 369–376. https://doi.org/10.1080/15226514.2014.910166
de Jong van Lier, Q., van Mentré, J. C., Durigon, A., dos Santos, M. A., & Metselaar, K. (2013). Modeling water potentials and ﬂows in thesoil‐plantsystem comparing hydraulic resistances and transpiration reduction functions. Vadose Zone Journal,
12(3). https://doi.org/10.2136/vzj2013.02.0039
Dumont, C., Mentré, F., Gaynor, C., Brendel, K., Gesson, C., & Chenel, M. (2013). Optimal sampling times for a drug and its metaboliteusing simcyp® simulations as prior information. Clinical Pharmacokinetics, 52(1), 43–57. https://doi.org/10.1007/s40262‐012‐0022‐9
Emerman, S. H. (1995). The tipping bucket equations as a model for macropore ﬂow. Journal of Hydrology, 171(1‐2), 23–47. https://doi.org/
10.1016/0022‐1694(95)02735‐8
U.S. EPA (2013). Estimation Program Interface (EPI) Suite, v4.1. United States Environ. Washington, DC, USA: Prot. Agency.
Eugenio, N. R., McLaughlin, M., & Pennock, D. (2018). Soil pollution: A hidden reality. Rome: FAO.
Fantke, P., & Juraske, R. (2013). Variability of pesticide dissipation half‐lives in plants. Environmental Science & Technology, 47(8),3548–3562. https://doi.org/10.1021/es303525x
FAO, 2019. The future of food safety‐First FAO/WHO/AU international food safety conference. Addis.
Fatichi, S., Manzoni, S., Or, D., & Paschalis, A. (2019). A Mechanistic model of microbially mediated soil biogeochemical processes: A
reality check. Global Biogeochemical Cycles, 33(6), 620–648. https://doi.org/10.1029/2018GB006077
Fatichi, S., Pappas, C., & Ivanov, V. Y. (2016). Modeling plant‐water interactions: An ecohydrological overview from the cell to the global
scale. Wiley Interdisciplinary Reviews Water, 3(3), 327–368. https://doi.org/10.1002/wat2.1125
Feddes, R. A., Kowalik, P. J., & Zaradny, H. (1978). Simulation of ﬁeld water use and crop yield.. Wageningen: PUDOC (Centre for agricultural publishing and documentation).
Gent, M. P. N. (2002). Growth and composition of salad greens as affected by organic compared to nitrate fertilizer and by environment inhigh tunnels. Journal of Plant Nutrition, 25(5), 981–998. https://doi.org/10.1081/PLN‐120003933

8986

Water Resources Research

10.1029/2019WR025432

Glenn, E. P., Cardran, P., & Thompson, T. L. (1984). Seasonal effects of shading on growth of greenhouse lettuce and spinach. Scientia
Horticulturae (Amsterdam)., 24(3‐4), 231–239. https://doi.org/10.1016/0304‐4238(84)90106‐7
Glick, B. R. (2010). Using soil bacteria to facilitate phytoremediation. Biotechnology Advances, 28(3), 367–374. https://doi.org/10.1016/j.
biotechadv.2010.02.001
Goldstein, M., Malchi, T., Shenker, M., & Chefetz, B. (2018). Pharmacokinetics in plants: Carbamazepine and its interactions with lamotrigine. Environmental Science & Technology, 52(12), 6957–6964. https://doi.org/10.1021/acs.est.8b01682
Hanson, B. R., Šimůnek, J., & Hopmans, J. W. (2006). Evaluation of urea–ammonium–nitrate fertigation with drip irrigation using
numerical modeling. Agricultural Water Management, 86(1‐2), 102–113. https://doi.org/10.1016/J.AGWAT.2006.06.013
Hartzell, S., Bartlett, M. S., & Porporato, A. (2017). The role of plant water storage and hydraulic strategies in relation to soil moistureavailability. Plant and Soil, 419(1‐2), 503–521. https://doi.org/10.1007/s11104‐017‐3341‐7
Herbst, M., Fialkiewicz, W., Chen, T., Pütz, T., Thiéry, D., Mouvet, C., et al. (2005). Intercomparison of ﬂow and transport models applied tovertical drainage in cropped lysimeters. Vadose Zone Journal, 4(2). https://doi.org/10.2136/vzj2004.0070
Hoffman, G. J., van Genuchten, M. T., & Metselaar, K. (1983). Soil properties and efﬁcient water use: Water management for salinity
control.. In H. M. Taylor, W. R. Jordan, & T. R. Sinclair (Eds.) Limitations to Efﬁcient Water Use in Crop Production. American Society ofAgronomy: Madison, WI(pp. 73–85). https://doi.org/10.2134/1983.limitationstoefﬁcientwateruse.c5
Hong, M. S., Farmayan, W. F., Dortch, I. J., Chiang, C. Y., McMillan, S. K., & Schnoor, J. L. (2001). Phytoremediation of MTBE from agroundwater plume. Environmental Science & Technology, 35(6), 1231–1239. https://doi.org/10.1021/es001911b
Huang, C. W., Domec, J. C., Palmroth, S., Pockman, W. T., Litvak, M. E., & Katul, G. G. (2018). Transport in a coordinated soil‐root‐xylem‐
phloem leaf system. Advances in Water Resources, 119, 1–16. https://doi.org/10.1016/j.advwatres.2018.06.002
Huang, C. W., Domec, J. C., Ward, E. J., Duman, T., Manoli, G., Parolari, A. J., & Katul, G. G. (2017). The effect of plant water storage onwater ﬂuxes within the coupled soil–plant system. New Phytologist, 213(3), 1093–1106. https://doi.org/10.1111/nph.14273
Huber, K., Vanderborght, J., Javaux, M., Schröder, N., Dodd, I. C., & Vereecken, H. (2014). Modelling the impact of heterogeneous rootzone
water distribution on the regulation of transpiration by hormone transport and/or hydraulic pressures. Plan and Soil. 384(1‐2), 93–112.
https://doi.org/10.1007/s11104‐014‐2188‐4
Hung, H., & Mackay, D. (1997). A novel and simple model of the uptake of organic chemicals by vegetation from air and soil. Chemosphere,
35(5), 959–977. https://doi.org/10.1016/S0045‐6535(97)00182‐3
Hurtado, C., Trapp, S., & Bayona, J. M. (2016). Inverse modeling of the biodegradation of emerging organic contaminants in the soil‐plantsystem. Chemosphere, 156, 236–244. https://doi.org/10.1016/j.chemosphere.2016.04.134
Jarvis, N., & Larsbo, M. (2012). MACRO (v5.2): Model use, calibration, and validation. Transactions of the ASABE, 55(4), 1413–1423.
https://doi.org/10.13031/2013.42251
Javaux, M., Couvreur, V., Vanderborght, J., & Vereecken, H. (2013). Root water uptake: From three‐dimensional biophysical processes tomacroscopic modeling approaches. Vadose Zone Journal, 12(4). https://doi.org/10.2136/vzj2013.02.0042
Jellali, S., Diamantopoulos, E., Haddad, K., Anane, M., Durner, W., & Mlayah, A. (2016). Lead removal from aqueous solutions by raw
sawdust and magnesium pretreated biochar: Experimental investigations and numerical modelling. Journal of Environmental
Economics and Management, 180, 439–449. https://doi.org/10.1016/j.jenvman.2016.05.055
Koba, O., Golovko, O., Kodešová, R., Klement, A., & Grabic, R. (2016). Transformation of atenolol, metoprolol, and carbamazepine in soils:
The identiﬁcation, quantiﬁcation, and stability of the transformation products and further implications for the environment.
Environmental Pollution, 218, 574–585. https://doi.org/10.1016/j.envpol.2016.07.041
Kodešová, R., Grabic, R., Kočárek, M., Klement, A., Golovko, O., Fér, M., et al. (2015). Pharmaceuticals' sorptions relative to properties ofthirteen different soils. Science of the Total Environment, 511, 435–443. https://doi.org/10.1016/j.scitotenv.2014.12.088
Kodešová, R., Klement, A., Golovko, O., Fér, M., Kočárek, M., Nikodem, A., & Grabic, R. (2019a). Soil inﬂuences on uptake and transfer ofpharmaceuticals from sewage sludge amended soils to spinach. Journal of Environmental Economics and Management, 250, 109407.
https://doi.org/10.1016/J.JENVMAN.2019.109407
Kodešová, R., Klement, A., Golovko, O., Fér, M., Nikodem, A., Kočárek, M., & Grabic, R. (2019b). Root uptake of atenolol, sulfamethoxazole and carbamazepine, and their transformation in three soils and four plants. Environmental Science and Pollution Research, 26(10),9876–9891. https://doi.org/10.1007/s11356‐019‐04333‐9
Kodešová, R., Kočárek, M., Klement, A., Golovko, O., Koba, O., Fér, M., et al. (2016). An analysis of the dissipation of pharmaceuticals
under thirteen different soil conditions. Science of the Total Environment, 544, 369–381. https://doi.org/10.1016/j.scitotenv.2015.11.085
Kumar, K., & Gupta, S. C. (2016). A framework to predict uptake of trace organic compounds by plants. Journal of Environmental Quality,
45(2). https://doi.org/10.2134/jeq2015.06.0261
Lange, C., Schneider, M., Mutz, M., Haustein, M., Halle, M., Seidel, M., et al. (2015). Model‐based design for restoration of a small urbanriver. Journal of Hydro‐environment Research, 9(2), 226–236. https://doi.org/10.1016/j.jher.2015.04.003
Larsen, M., Trapp, S., & Pirandello, A. (2004). Removal of cyanide by woody plants. Chemosphere, 54(3), 325–333. https://doi.org/10.1016/
S0045‐6535(03)00662‐3
Legind, C. N., Kennedy, C. M., Rein, A., Snyder, N., & Trapp, S. (2011). Dynamic plant uptake model applied for drip irrigation of an
insecticide to pepper fruit plants. Pest Management Science, 67(5), 521–527. https://doi.org/10.1002/ps.2087
Legind, C. N., Rein, A., Serre, J., Brochier, V., Haudin, C. S., Cambier, P., et al. (2012). Simultaneous simulations of uptake in plants andleaching to groundwater of cadmium and lead for arable land amended with compost or farmyard manure. PLOS One, 7(10), e47002.
https://doi.org/10.1371/journal.pone.0047002
Li, H., Sheng, G., Chiou, C. T., & Xu, O. (2005). Relation of organic contaminant equilibrium sorption and kinetic uptake in plants.
Environmental Science & Technology, 39(13), 4864–4870. https://doi.org/10.1021/es050424z
Malchi, T., Maor, Y., Tadmor, G., Shenker, M., & Chefetz, B. (2014). Irrigation of root vegetables with treated wastewater: Evaluating
uptake of pharmaceuticals and the associated human health risks. Environmental Science & Technology, 48(16), 9325–9333. https://doi.
org/10.1021/es5017894
Manoli, G., Bonetti, S., Domec, J. C., Putti, M., Katul, G., & Marani, M. (2014). Tree root systems competing for soil moisture in a 3D soil‐
plant model. Advances in Water Resources, 66, 32–42. https://doi.org/10.1016/j.advwatres.2014.01.006
Manoli, G., Huang, C. W., Bonetti, S., Domec, J. C., Marani, M., & Katul, G. (2017). Competition for light and water in a coupled soil‐plantsystem. Advances in Water Resources. https://doi.org/10.1016/j.advwatres.2017.08.004
Mattina, M. I., Lannucci‐Berger, W., Musante, C., & White, J. C. (2003). Concurrent plant uptake of heavy metals and persistent organic
pollutants from soil. Environmental Pollution, 124(3), 375–378. https://doi.org/10.1016/s0269‐7491(03)00060‐5
Morris, M. D. (1991). Factorial sampling plans for preliminary computational experiments. Technometrics, 33(2), 161–174. https://doi.org/
10.2307/1269043

BRUNETTI ET AL.

8987

Water Resources Research

10.1029/2019WR025432

Nolan, B. T., Randal Bayless, E., Green, C. T., Garg, S., Voss, F. D., Lampe, D. C., et al. (2005). Evaluation of unsaturated‐zone solute‐transport models for studies of agricultural chemicals. Reston, VA: U.S. Geological Survey.
Ouyang, Y. (2002). Phytoremediation: Modeling plant uptake and contaminant transport in the soil‐plant‐atmosphere continuum. Journalof Hydrology, 266(1‐2), 66–82. https://doi.org/10.1016/S0022‐1694(02)00116‐6
Peters, A., & Durner, W. (2008). Simpliﬁed evaporation method for determining soil hydraulic properties. Journal of Hydrology, 356(1‐2),
147–162. https://doi.org/10.1016/j.jhydrol.2008.04.016
Rein, A., Legind, C. N., & Trapp, S. (2011). New concepts for dynamic plant uptake models. SAR QSAR Environment Research, 22(1‐2),
191–215. https://doi.org/10.1080/1062936X.2010.548829
Ritchie, J. T. (1972). Model for predicting evaporation from a row crop with incomplete cover. Water Resources Research, 8(5), 1204–1213.
https://doi.org/10.1029/WR008i005p01204
Ryan, J. A., Bell, R. M., Davidson, J. M., & O'Connor, G. A. (1988). Plant uptake of non‐ionic organic chemicals from soils. Chemosphere,
17(12), 2299–2323. https://doi.org/10.1016/0045‐6535(88)90142‐7
Sabourin, L., Duenk, P., Bonte‐Gelok, S., Payne, M., Lapen, D. R., & Topp, E. (2012). Uptake of pharmaceuticals, hormones and parabens
into vegetables grown in soil fertilized with municipal biosolids. Science of the Total Environment, 431, 233–236. https://doi.org/10.1016/
j.scitotenv.2012.05.017
Saltelli, A., & Annoni, P. (2010). How to avoid a perfunctory sensitivity analysis. Environmental Modelling & Software, 25(12), 1508–1517.
https://doi.org/10.1016/j.envsoft.2010.04.012
Sandermann, H. (1992). Plant metabolism of xenobiotics. Trends in Biochemical Sciences, 17(2), 82–84. https://doi.org/10.1016/0968‐0004(92)90507‐6
Sandermann, H. (1994). Higher plant metabolism of xenobiotics: The ‘green liver’ concept. Pharmacogenetics, 4(5), 225–241. https://doi.
org/10.1097/00008571‐199410000‐00001
Scanlon, B. R., Christman, M., Reedy, R. C., Porro, I., Simunek, J., & Flerchinger, G. N. (2003). Intercode comparisons for simulating waterbalance of surﬁcial sediments in semiarid regions. Water Resources Research, 38(12), WR001233. https://doi.org/10.1029/2001wr001233
Scanlon, B. R., Keese, K., Reedy, R. C., Simunek, J., & Andraski, B. J. (2003). Variations in ﬂow and transport in thick desert vadose zones in
response to paleoclimatic forcing (0‐90 kyr): Field measurements, modeling, and uncertainties. Water Resources Research, 39(7),WR001604. https://doi.org/10.1029/2002WR001604
Schaap, M. G., Leij, F. J., & Van Genuchten, M. T. (2001). Rosetta: A computer program for estimating soil hydraulic parameters withhierarchical pedotransfer functions. Journal of Hydrology, 251(3‐4), 163–176. https://doi.org/10.1016/S0022‐1694(01)00466‐8
Schroeder, P. R.Aziz, N. M., Lloyd, C. M., & Zappi, P. A. (1994). The Hydrologic Evaluation of Landﬁll Performance (HELP) Model:User's
Guide for Version 3, EPA/600/R‐94/168a. Washington, DC: U.S. Environmental Protection Agency Ofﬁce of Research and Development.
Shelia, V., Šimunek, J., Boote, K., & Hoogenbooom, G. (2018). Coupling DSSAT and HYDRUS‐1D for simulations of soil water dynamics inthe soil‐plant‐atmosphere system. Journal of Hydrology Hydromechanics, 66(2), 232–245. https://doi.org/10.1515/johh‐2017‐0055
Šimůnek, J., & Hopmans, J. W. (2009). Modeling compensated root water and nutrient uptake. Ecological Modelling, 220(4), 505–521.
https://doi.org/10.1016/j.ecolmodel.2008.11.004Šimůnek, J., van Genuchten, M. T., & Šejna, M. (2016). Recent developments and applications of the HYDRUS computer software
packages. Vadose Zone Journal, 15(7). https://doi.org/10.2136/vzj2016.04.0033
So, E. L., Ruggles, K. H., Cascino, G. D., Ahmann, P. A., & Weatherford, K. W. (1994). Seizure exacerbation and status epileptics related tocarbamazepine‐10,11‐epoxide. Annals of Neurology, 35(6), 743–746. https://doi.org/10.1002/ana.410350616
Stroh, M., Hutmacher, M. M., Pang, J., Lutz, R., Magara, H., & Stone, J. (2013). Simultaneous pharmacokinetic model for rolofylline andboth M1‐trans and M1‐cis metabolites. The AAPS Journal, 15(2), 498–504. https://doi.org/10.1208/s12248‐012‐9443‐5
Sutanto, S. J., Wenninger, J., Coenders‐Gerrits, A. M. J., & Uhlenbrook, S. (2012). Partitioning of evaporation into transpiration, soil evaporation and interception: A comparison between isotope measurements and a HYDRUS‐1D model. Hydrology and Earth System
Sciences, 16(8), 2605–2616. https://doi.org/10.5194/hess‐16‐2605‐2012
Taylor, S. A., & Ashcroft, G. L. (1972). Physical edaphology: The physics of irrigated and nonirrigated soils. San Francisco: W.H. Freeman.
Thompson, N. (1983). Diffusion and uptake of chemical vapour volatilising from a sprayed target area. Pesticide Science, 14(1), 33–39.
https://doi.org/10.1002/ps.2780140106
Topp, E., Scheunert, I., Attar, A., & Korte, F. (1986). Factors affecting the uptake of 14C‐labeled organic chemicals by plants from soil.
Ecotoxicology and Environmental Safety, 11(2), 219–228. https://doi.org/10.1016/0147‐6513(86)90066‐7
Trapp, S. (2007). Fruit tree model for uptake of organic compounds from soil and air. SAR and QSAR in Environmental Research, 18(3‐4),
367–387. https://doi.org/10.1080/10629360701303693
Trapp, S. (2009). Bioaccumulation of polar and ionizable compounds in plants. In J. Devillers (Ed.), Ecotoxicology modeling. Emerging
topics in ecotoxicology (Principles, approaches and perspectives), Vol, (Vol. 2, pp. 299–353). Boston, MA: Springer. https://doi.org/
10.1007/978‐1‐4419‐0197‐2_11
Trapp, S. (2015). Calibration of a plant uptake model with plant‐ and site‐speciﬁc data for uptake of chlorinated organic compounds into
radish. Environmental Science & Technology, 49(1), 395–402. https://doi.org/10.1021/es503437p
Trapp, S., & Eggen, T. (2013). Simulation of the plant uptake of organophosphates and other emerging pollutants for greenhouse experiments and ﬁeld conditions. Environmental Science and Pollution Research, 20(6), 4018–4029. https://doi.org/10.1007/s11356‐012‐1337‐7
Trapp, S., Feiﬁcova, D., Rasmussen, N. F., & Bauer‐Gottwein, P. (2008). Plant uptake of NaCl in relation to enzyme kinetics and toxic
effects. Environmental and Experimental Botany, 64(1), 1–7. https://doi.org/10.1016/j.envexpbot.2008.05.001
Trapp, S., & Matthies, M. (1998). Chemodynamics and environmental modeling, chemodynamics and environmental modeling: An introduction. Berlin, Heidelberg: Springer Berlin Heidelberg. https://doi.org/10.1007/978‐3‐642‐80429‐8Travis, C. C., & Arms, A. D. (1988). Bioconcentration of organics in beef, milk, and vegetation. Environmental Science & Technology, 22(3),271–274. https://doi.org/10.1021/es00168a005
Twarakavi, N. K. C., Sakai, M., & Šimůnek, J. (2009). An objective analysis of the dynamic nature of ﬁeld capacity. Water Resources
Research, 45(10), W10410. https://doi.org/10.1029/2009WR007944
van Genuchten, M. T. (1980). A closed‐form equation for predicting the hydraulic conductivity of unsaturated soils1. Soil Science Society ofAmerica Journal, 44(5), 892–898. https://doi.org/10.2136/sssaj1980.03615995004400050002x
Vereecken, H., Schnepf, A., Hopmans, J.W., Javaux, M., Or, D., Roose, T., et al., 2016. Modeling soil processes: Review, key challenges, andnew perspectives. Vadose Zone Journal https://doi.org/10.2136/vzj2015.09.0131, 15, 5
Vogel, T., Votrubova, J., Dohnal, M., & Dusek, J. (2017). A simple representation of plant water storage effects in coupled soil water ﬂow
and transpiration stream modeling. Vadose Zone Journal, 16(5). https://doi.org/10.2136/vzj2016.12.0128

BRUNETTI ET AL.

8988

Water Resources Research

10.1029/2019WR025432

Warner, T., Patsalos, P.N., Prevett, M., Elyas, A.A., Duncan, J.S., 1992. Lamotrigine‐induced carbamazepine toxicity: An interaction withcarbamazepine‐10,11‐epoxide. Epilepsy Research https://doi.org/10.1016/0920‐1211(92)90049‐Y
Wen, B., Wu, Y., Zhang, H., Liu, Y., Hu, X., Huang, H., & Zhang, S. (2016). The roles of protein and lipid in the accumulation and distribution of perﬂuorooctane sulfonate (PFOS) and perﬂuorooctanoate (PFOA) in plants grown in biosolids‐amended soils. Environmental
Pollution, 216, 682–688. https://doi.org/10.1016/j.envpol.2016.06.032
Werner, C., & Dubbert, M. (2016). Resolving rapid dynamics of soil‐plant‐atmosphere interactions. New Phytologist, 210(3), 767–769.
https://doi.org/10.1111/nph.13936
Wesseling, J., Elbers, J., Kabat, P., Van den Broek, B. , 1991. SWATRE: Instructions for input. Intern. Note, Winand Star. Cent.
Yang, J., Yu, Q., & Gong, P. (2008). Quantifying air pollution removal by green roofs in Chicago. Atmospheric Environment, 42(31),7266–7273. https://doi.org/10.1016/j.atmosenv.2008.07.003
Yu, X., Trapp, S., Zhou, P., Wang, C., & Zhou, X. (2004). Metabolism of cyanide by Chinese vegetation. Chemosphere, 56(2), 121–126.
https://doi.org/10.1016/j.chemosphere.2004.02.008
Zeng, C. Z., Bie, Z. L., & Yuan, B. Z. (2009). Determination of optimum irrigation water amount for drip‐irrigated muskmelon (Cucumis
melo L.) in plastic greenhouse. Agricultural Water Management, 96(4), 595–602. https://doi.org/10.1016/j.agwat.2008.09.019
Zhang, H., Chen, J., Ni, Y., Zhang, Q., & Zhao, L. (2009). Uptake by roots and translocation to shoots of polychlorinated dibenzo‐p‐dioxins
and dibenzofurans in typical crop plants. Chemosphere, 76(6), 740–746. https://doi.org/10.1016/j.chemosphere.2009.05.030

BRUNETTI ET AL.

8989