Agricultural Water Management 254 (2021) 106971

Contents lists available at ScienceDirect

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

Numerical modeling to optimize nitrogen fertigation with consideration oftransient drought and nitrogen stressThomas Groenveld a, *, Amir Argaman b, Jiří Šimůnek c, Naftali Lazarovitch ba

The Albert Katz International School for Desert Studies, The Jacob Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev, Sede Boqer Campus,
8499000, IsraelbFrench Associates Institute for Agriculture and Biotechnology of Drylands, Jacob Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev, Sede Boqer
Campus, 8499000, Israelc
Department of Environmental Sciences, University of California Riverside, Riverside, CA, USA

A R T I C L E I N F O

A B S T R A C T

Handling Editor: Dr Z Xiying

Optimization of nitrogen (N) fertigation is a formidable challenge involving complex interactions between waterand N uptake and their effects on crop production. Numerical models can be useful in studying the interaction ofmultiple variables like those found in mechanistic simulations of N fertigation strategies. The physical aspects
can often be accurately represented in soil-plant-atmosphere continuum models, while the biological factors lag
due to their oversimplification. When optimizing N fertigation using numerical models, it is essential to consider
the effects of N and water stresses on the plant size and corresponding feedback on potential transpiration and N
uptake. The HYDRUS (2D/3D) model was modified to allow for active uptake and decay of multiple solutes andreduce potential transpiration due to a limitation in N uptake. Subsequently, we calibrated and validated themodel with a dataset that consisted of 3 nitrate (NO-3) concentration and 6 irrigation levels: a total of 18 distinct
treatments used to fertigate cucumber plants grown in lysimeters. The calibration was based on the treatmentthat received the highest N fertigation. The model was validated by testing its ability to accurately reduce po­tential N uptake and transpiration in water and N deficiency cases. Simulations showed that the N stress function
could explain 82% of the reduction in transpiration measured in the experimental setup. The sensitivity analysis,
evaluating the effects of the root shape and distribution parameters by increasing and decreasing their values by
20%, showed that these parameters had little impact on the results. Following its validation, the model was used
to determine the optimal N concentration in irrigation water and the optimal N application amount to obtainmaximal yield with minimal N loss. The optimal irrigation water NO-3-N concentration and seasonal NO-3-N
application were determined to be 75 mg L− 1 and 40 mg m− 2, respectively.

Keywords:
Transpiration
Nitrogen use efficiency
Active uptake
Nitrate leaching
Horticulture

Other parameter definitions are listed in Table 2.
1. Introduction
Global nitrate (NO-3) leaching is a function of fertilizer applications(Wang et al., 2019) and is correlated to the size of the world population(Erisman et al., 2008). The need for nitrogen (N) fertilizer to meet thegrowing global food demand and the environmental concerns associatedwith NO-3 leaching (Thompson et al., 2020) have led to the recognition
that many of the world’s sustainability goals can be achieved by opti­
mizing agricultural N use (Houlton et al., 2019). Among the agricultural
crops, horticultural crops are prone to the highest NO-3 losses due to highfertilizer application rates and their relatively shallow root zone (Padilla
et al., 2018). Fertigation (fertilizer application through irrigation water)

has been demonstrated to maximize plant N uptake while minimizing
NO-3 leaching (Dahan et al., 2014; Groenveld et al., 2019). Matching the
N supply to the N demand in space and time is a complex problem that
depends on many factors, including the plant water and N demand, theroot distribution, the soil hydraulic and solute transport properties, theN concentration of the irrigation water, the timing of the fertigation, andthe total amount of water and N applied. While simple analytical modelsof yield response to N application are helpful (Cerrato and Blackmer,
1990), they often give contradicting results concerning optimal N
application (Groenveld et al., 2019) and only address the total amount of
N that should be applied (Gallardo et al., 2020). With larger data sets,
more complex analytical models can be developed (Głąb et al., 2020),
but these lack the flexibility to be applied to different crops or envi­
ronmental conditions. The complex interactions in space and time of the* Corresponding author.
E-mail addresses: tomgreen@post.bgu.ac.il (T. Groenveld), jiri.simunek@ucr.edu (J. Šimůnek), lazarovi@bgu.ac.il (N. Lazarovitch).
https://doi.org/10.1016/j.agwat.2021.106971
Received 19 January 2021; Received in revised form 11 May 2021; Accepted 13 May 2021
0378-3774/© 2021 Elsevier B.V. All rights reserved.

T. Groenveld et al.

Agricultural Water Management 254 (2021) 106971

uptake per surface area or increasing the total root area. It is possible toestimate each of these pathways’ contributions using detailed models
such as R-SWMS in combination with measured data. However, thesimultaneous consideration of all of these three factors is limited due tomany parameters required (Rengel, 1993) and the difficulty of deter­mining them (Roose and Schnepf, 2008).
The passage of nutrients from the soil solution into the plant can be
described by passive or active mechanisms, which refer to the mass flowof nutrients into roots with water taken up for transpiration or againstenergy gradients, respectively. Barber (1962) estimated passive uptake
by multiplying the volume of transpired water with the nutrient’s mean
concentration in the bulk soil solution. However, this method may
overestimate an actual ion supply by water mass flow (Jungk, 2002).
Epstein and Hagen (1952) proposed using the Michaelis-Menten model
to describe the characteristics of individual membrane carriers involvedin active nutrient uptake. The uptake rate per root area differs according
to the root order (Rewald et al., 2011). The root order changes as roots
grow, and there is potential compensation of reduced nutrient uptake
from nutrient-poor areas by increased uptake from less nutrient stressed
areas of the root zone (Rubio and Lynch, 2007).
The macroscopic numerical model HYDRUS (2D/3D) considers
transient water flow, solute transport, and compensated plant root waterand nutrient uptake (Šimůnek and Hopmans, 2009). These are calcu­
lated as a function of soil hydraulic properties, irrigation salinity andamount, root distribution, potential evapotranspiration, and potential
nutrient uptake, all of which are described in detail in the model tech­nical and user manuals (Šejna et al., 2018; Šimůnek et al., 2018b). The
water flow algorithm in HYDRUS (2D/3D) has been widely tested andused for different applications (Šimůnek et al., 2016), including root
water uptake (Zhu et al., 2009), N transport and uptake (Bristow et al.,2020; Gärdenäs et al., 2005; Karandish and Šimůnek, 2017; Tao et al.,2021; Wang et al., 2014), and the effects of salinity (Phogat et al., 2014;
Ramos et al., 2012; Shouse et al., 2011; Skaggs et al., 2006).
Wang and Baerenklau (2014) linked crop yield to relative transpi­
ration and N uptake in the dataset they generated with HYDRUS-1D
using a function of Pang and Letey (1998). Levers and Schwabe(2017) took this a step further by calculating the instantaneous yield
reduction by comparing potential and actual water and passive N uptake
throughout the growing season instead of only at the end, using amodified version of HYDRUS-1D. Rezaei et al. (2017) used a plant
growth model [LINGRA-N (Wolf, 2012)] to determine the effect of N
stress on the plant’s potential uptake, which was then fed back to
HYDRUS-1D. The effects of limited nutrient uptake on growth (i.e., atranspiration reduction) have not yet been implemented in the standard
version of HYDRUS. Because biomass accumulation is a function of bothtranspiration and N uptake and determines their potential rates, theseprocesses should be linked so that stress in either process affects theother’s potential rate.
Although the active nutrient uptake function was implemented into
the HYDRUS (2D/3D) model more than ten years ago, it has only beenvalidated for NH+
4 uptake in flooded rice paddies (Li et al., 2015; Yanget al., 2017). While it is understood that nitrate (NO-3) uptake by plantroots can be described by Michaelis-Menten kinetics (Tischner andKaiser, 2007), all published HYDRUS simulations validated with NO-3
uptake data do so using the passive uptake mechanism only. This might
be due to the lack of a function in the model that reduces plant growth
due to N stress, representing a formidable research gap in modeling N
transport and uptake.
Our research objective is to validate the combined passive, active,
and compensated NO-3uptake functions in HYDRUS (2D/3D) using
measured data while considering a function that reduces potential
transpiration due to N stress. Subsequently, we demonstrate the use ofthis model as a tool for optimizing N fertigation in terms of N use effi­
ciency (NUE), yield, and NO-3 leaching.

Nomenclature
A30
A60
A90
Ap
DAT
dP
EC
I/T
N
N2
N2 O
NO
nRMSE
Rp
S
Tac
Tp
u

η
ρ

30 mg L− 1 NO-3-N treatment
60 mg L− 1 NO-3-N treatment
90 mg L− 1 NO-3-N treatment
Potential active uptake (mg day− 1)
Days after transplant
Effective pore diameter (m)Electrical conductivity (dS m− 1)
Irrigation/Transpiration (-)
Nitrogen
Dinitrogen
Nitrous oxide
Nitric oxide
Normalized root mean squared error (-)
Potential solute uptake (mg day− 1)
Substrate concentration (mg cm− 2)
Actual transpiration (cm2 day− 1)
Potential transpiration (cm2 day− 1)
Specific discharge (m s− 1)
Dynamic viscosity (N s m− 2)
Density of water (kg m− 3)

variables mentioned above make the optimization of N fertigation ahighly non-linear problem, which requires numerical models toapproximate a solution (Azad et al., 2020). An increase in computer
processing power and availability of detailed datasets offers the possi­
bility of opening up empirical “black boxes” in numerical models so thatprocesses can be described in an increasingly mechanistic fashion and ingreater detail. A more mechanistic model will be relevant to a wider
variety of plants and will be flexible to changes in the soil type or irri­
gation frequency.
Several models have been developed to describe N leaching fromagricultural systems. The DNDC model simulates crop growth anddevelopment, as well as carbon (C) and N biogeochemistry in agricul­
tural ecosystems (He et al., 2019). The COUP model is a versatile plat­form of different modules that can be used to simulate a wide range of N
transformations, including plant interactions with mycorrhizal fungi(He et al., 2020). A limitation of these models is that they are
one-dimensional. One-dimensional models are suitable to model only
vegetations horizontally invariant, such as pasture (Eshel and Beeck­man, 2013). However, even in grass fields, lateral dimensions have to beconsidered under certain conditions, such as with sloping land or landwith urine patches left by grazing cattle (Snow et al., 2017). While
multiple dimensions significantly increase the computational timerequired for simulations, they are often essential to describe thecomplexity of irrigated agricultural systems that are not horizontallyinvariant (Russo et al., 2013). Single-root models simulate the rootsystem as a cylinder taking up water and nutrients from the soil volume
around it. An advantage of such models is that flow and transport
through the soil and across the root membrane can be calculated sepa­
rately, yielding a more mechanistic formulation (Eshel and Beeckman,
2013; McMurtrie and Näsholm, 2018). However, the difficulty of
upscaling the results of such models is a disadvantage (Roose andSchnepf, 2008), particularly when studying processes at the root zone
scale.
To simulate solute uptake by roots, the solute’s movement to the root
surface (Rengel, 1993), the root growth, and the demand-controlledaspect of uptake should be considered (Hopmans and Bristow, 2002).
While the first two factors are receiving increasing attention in models
like R-SWMS (Javaux et al., 2008; Schröder et al., 2013), the third one (i.
e., the demand-control aspect) has seen limited implementation in 3Dmodels. An N-limited plant can increase N uptake by increasing the2

T. Groenveld et al.

Agricultural Water Management 254 (2021) 106971

2. Methods

transparent columns with a 15 cm diameter. The columns were slowly
wetted from below to prevent the trapping of air in the columns. The
time it took for a ponded layer of 15 cm of water to drop 12 cm wasmeasured four times for three different columns filled with perlite at abulk density of 0.06 g cm− 3.
A retention curve was obtained using a PVC column of 11 cm indiameter and 100 cm in length that was cut into 5 cm slices. The slices
were stacked to form a column that was sealed airtight with plastic tape.
The column was packed with perlite at the same density as in thegreenhouse experiment (0.06 g cm− 3) and gradually wetted from below
over a 24 h period to ensure that the air was pushed out of the column.
Subsequently, the entire column was removed from the water bath,placed on an electronic scale, and covered with a lid to eliminate
evaporation. The mass readings were recorded over time withdecreasing frequency as less water drained out of the column. When thedrainage water flux stopped after a few days and the column was
assumed to be at a hydrostatic state, the 5 cm slices were cut off thecolumn, wrapped in a silver foil, and weighed. These slices were thendried at 105 ◦ C and reweighed, resulting in a retention curve consistingof 20 pairs of height and volumetric water content (Raij et al., 2016).
This experiment was repeated four times.
Breakthrough curves were obtained for all 18 plots to determine thelongitudinal and transverse dispersivities (λL and λT) of the growing
media used. Irrigation water was applied in 3 L pulses at three evenly
spaced locations in the plot’s centerline where the plants would later belocated. All plots were flushed two days with saline irrigation water usedin all experiments (EC: 2.28 ± 0.15 dS m− 1). Then, on the third day, they
were irrigated with five 3 L pulses of distilled water (EC: 0.11 ± 0.05 dS
m− 1), after which irrigation continued with saline water (as before).
After every 3 L irrigation, pulse drainage water was collected for all 18
plots, and its electrical conductivity was measured.

2.1. Field experiment and measurements
2.1.1. Experimental design
The data used to calibrate and validate the model was taken from
Groenveld et al. (2019), where a more detailed description of theexperimental setup can be found. The effluent water of three recircu­
lating aquaculture systems with NO-3-N concentrations of 30, 60, and 90
mg L− 1 (named A30, A60, and A90, respectively) were used to fertigate
six lysimeters per each N treatment. Three cucumber plants spaced 40
cm apart were grown in each lysimeter in a single row. Each lysimeterwas 100 cm long, 50 cm wide, and 20 cm deep. The growing medium
used was perlite (Perlite 206, Agrikal, Israel), whose grain size distri­bution is listed in Table 1. Aside from N, all other nutrients were kept at
equal concentrations in all treatments. Irrigation water was cooled to24 ◦ C and applied through a 6 mm tube to the base of each of the threeplants at a combined flow rate of 0.67 mm/min. Drainage water wascollected and weighed daily. Transpiration was calculated by subtract­
ing drainage from irrigation, ignoring a possible change in the root zone
water content (Tripler et al., 2012). Each lysimeter had a unique irri­
gation regime, a multiple of measured transpiration (irrigation/­transpiration: I/T between 1 and 5). Overall, there were 18 distinct
fertigation treatments of cucumber plants grown in perlite media(Fig. 1A). Because plant transpiration determined irrigation while theNO-3 concentration remained the same, N fertigation was indirectly also
a function of plant transpiration.
2.1.2. Measurements and analysis
Irrigation water was sampled daily and drainage water of all 18 ly­
simeters three times a week to measure the electrical conductivity (EC)
(DDS 120W, Bante Instruments, China) and NO-3-N concentration(Reflectoquant Nitrate Test 116971, Merck KGaA, Germany and OrionStar A214, Thermo Fisher Scientific, MA, USA). Ripe cucumbers were
harvested daily, starting 20 days after transplanting (DAT). Cucumbers
were weighed individually for each of the three plants in each plot. All
plants were destructively harvested at 40 DAT. The dry biomass of fruits,
stems, and leaves was measured and then tested for N content (OEACHNS Flash, 2000, Thermo Fisher Scientific, MA, USA).
The water and N mass balance components included a precise recordof the amounts of water applied, their concentrations, each fertigation
event’s timing, and daily records of transpiration and NO-3 leached out of
the lysimeter. The growth rate and yield of the cucumber plants wererecorded daily. The total uptake of N and its distribution in different
plant parts was determined at the end of the experiment. Once theexperiment was terminated, the root zone was cut into three blocks ofequal sizes, perlite was washed off the roots through a 3-mm sieve, andthe roots were then dried and weighed. The root N content was notmeasured, but an estimated value of 5% was used for all treatments insubsequent calculations (Groenveld et al., 2019). The amounts of N
taken up by the plants and drained out of the root zone as NO-3 are re­ported for each of the 18 treatments as a function of the amount of
fertigated NO-3-N in Fig. 1B. The difference between the amount of N
fertigated and the sum of N taken up and drained is assumed to be due toemissions of gaseous forms of N.

2.2. Modeling2.2.1. Model description
HYDRUS (2D/3D) describes unsaturated flow using the Richards
equation, solute transport using the convection-dispersion equation, andsolute transformation using first-order decay reactions. Potential plantwater and nutrient uptakes are specified in HYDRUS (2D/3D) as input
time-variable boundary conditions distributed across the root zone ac­cording to a root density distribution. Over time, the root distribution
can be fixed, entered as an array of values, described using a logistic
growth function (Šimůnek et al., 2018a), or linked to various environ­mental factors (Hartmann et al., 2017). Potential water uptake can be
reduced due to drought and salinity stresses. However, these stress ef­
fects do not accumulate over time in the current version of HYDRUS(2D/3D). Nutrient uptake is considered as passive or active mechanisms.
Active uptake is described by a modified version of theMichaelis-Menten equation:
Aa =Table 1
Grain size distribution of the perlite used in lysimeters, as reported by the
manufacturer (Agrikal, Israel).
< 0.36

0.36–0.5

0.5–1.0

1.0–1.4

1.4–2.0

> 2.0

%

6

7

49

21

16

1

(1)

where Aa is the actual active solute uptake rate (mg day− 1), Ap is thepotential active solute uptake rate, which is the remainder of total po­
tential uptake after subtracting passive uptake (mg day− 1), S is thesubstrate concentration (mg cm− 2), Smin is the minimal concentration for
active uptake (mg cm− 2), and KM (mg cm− 2) is the Michaelis-Menten
coefficient, which is equal to the substrate concentration correspond­
ing to half of the maximum solute uptake rate (Jungk, 2002). The
maximum solute uptake rate and the KM parameter for theMichaelis-Menten model are usually determined based on solute uptake
per root length or area. Using a dimensionless root distribution makes it
possible to cancel these root area units. Kim and Lieth (2012) scaled the
maximum solute uptake rate in their model to the total plant require­ment by multiplying the Michaelis-Menten equation by the plant dry

2.1.3. Parameter estimation experiments
The saturated hydraulic conductivity (KS) was measured with thefalling head method (Wever et al., 2004) using 20 cm cores of perlite inmm

S − SminAp
KM + S − Smin

3

T. Groenveld et al.

Agricultural Water Management 254 (2021) 106971

Fig. 1. (A) The measured amount of nitrate applied as a function of the irrigation water concentration and I/T rates. (B) The measured total amount of N found in theplant and drainage water as a function of the applied nitrate. Each color represents a nitrate concentration treatment, and the hue of each color represents the I/Ttreatment whose value can be read from the x-axis of (A). (For interpretation of the references to colour in this figure legend, the reader is referred to the web versionof this article.)(Reprinted from Journal of Cleaner Production, Vol 212, Groenveld, T., Kohn, Y.Y., Gross, A., Lazarovitch, N., 2019. Optimization of N use efficiency by means offertigation management in an integrated aquaculture-agriculture system. Pages 401–408, Copyright 2019, with permission from Elsevier).

matter over time. Sands and Smethurst (1995) also linked the maximum
solute uptake rate linearly to the plant water uptake rate.
Active solute uptake can be limited by limited uptake of water(Šimůnek and Hopmans, 2009):[]
Tac (t)Ap (t) = max Rp (t)(2)− Pa (t), 0
Tp (t)
where Ap is potential active uptake (mg day− 1) at a certain point in time(t), Rp is potential solute uptake (mg day− 1), Tac and Tp are actual andpotential transpiration (cm2 day− 1), and Pa is actual passive solute up­take (mg day− 1). The cm2 units for transpiration are used for twodimensional simulations; they can be converted to cm3 by multiplying
by the plot’s length. Compensation of water or solute uptake in a lessstressed area of the root zone can be enabled according to the degree
determined by the user.
2.2.2. Domain geometry and fixed boundary conditions
A two-dimensional flow domain 18 cm deep and 25 cm wide wasdiscretized using a finite element mesh of 350 nodes, with smaller ele­ments at the points where water entered and left the domain (Lazar­
ovitch et al., 2005) (Fig. 2). As the discretization had a large effect on themass balance error, the number of nodes was optimized to minimize thewhole data set’s mass balance error. The left side of the domain forms aline of symmetry in the middle of the plot, exactly under the dripline andthe plant’s base, considered as line features. The distance between plants
in the lysimeter was 40 cm, and the drainage holes in the drainage
trench at the bottom of the lysimeter were spaced at 15 cm. Emitters at
the dripline were spaced at 40 cm, which resulted in a small overlap of
the wetting patterns (Kandelous et al., 2011). A time-variable flux
boundary condition was used to represent the ponded area under theirrigation water emitter. Drainage outflow at the bottom occurred
through a 0.5 cm wide seepage face boundary condition. The sides of thelysimeter were represented by a no-flow boundary condition (Fig. 2).

Fig. 2. The root density distribution at DAT 24, after which it stayed the same.
The left side of the domain forms a symmetry line in the middle of the 1 m longplot. There is a 3 cm wide atmospheric boundary condition at the top left corner
where the plant base is located. The bottom of the drainage channel, at thelowest point, has a 0.5 cm wide seepage face boundary condition. The
remaining boundaries are no flow. The domain discretization of 350 nodes canbe seen as faint blue lines. (For interpretation of the references to colour in thisfigure legend, the reader is referred to the web version of this article.)weighting. The molecular diffusion coefficient is usually negligible insoils (Radcliffe and Šimůnek, 2010) and was therefore set to zero for the
perlite medium used. The longitudinal and transverse dispersivities (λLand λT) were calculated from collected breakthrough curves using theinverse mode of HYDRUS (2D/3D). The irrigation water amount andsalinity were used as the time-variable boundary conditions at the top of
the domain. The drainage water electrical conductivity, measured after
each irrigation pulse, was used as the data for the inverse solution. The
dispersivity was determined as the average value of the 18 plots.
The S-shaped water stress function (van Genuchten, 1987) was used
to represent water stress with parameters p50 of − 200 cm and pw3 of 3.
This relatively high value for the p50 parameter was used to ensure highsensitivity to water stress; the typical value for horticultural crops is

2.2.3. Input parameters
Model-specific parameters are required to describe water flow, solutetransport and transformations, and root water and solute uptake. The
van Genuchten-Mualem model (van Genuchten, 1980) was used todescribe soil hydraulic properties, and hysteresis was not considered.
RETC 6.02 (van Genuchten et al., 1991) was used to determine thesaturated water content (Ɵs), residual water content (Ɵr), and the shapeparameters α and n.
For solute transport, the Crank-Nicholson scheme was used for timeweighting, and Galerkin finite element scheme was used for space
4T. Groenveld et al.

Agricultural Water Management 254 (2021) 106971

lower (Groenveld et al., 2013). The S-shaped salinity stress function (vanGenuchten and Hoffman, 1984) was used to model the transpiration
response to salinity, with a c50 value of 3 dS m− 1 and a ps3 value of 3.
These values are typical for cucumber (Maas and Hoffman, 1977). Theinteraction between salinity and water stresses was considered multi­
plicative. The osmotic coefficient (φ) was set to 1 as the irrigation watersalinity units (dS m− 1) were the same as those of the stress response
function (Šejna et al., 2018). The maximum possible N concentration ofwater taken up by the roots (for passive uptake) was set at 1000 mg L− 1
to allow for unlimited passive uptake. The KM parameter of theMichaelis-Menten equation (Eq. (1)) for NO-3-N was assumed to be
1 mg L− 1, which is at the lower end of the 0.3–6 mg L− 1 range reported
for high-affinity transport systems of NO-3-N in horticultural crops (LeBot and Adamowicz, 2006).

Tp− new = Tp− input

( )
Ra
Rp

(4)

where Tp-new is adjusted potential transpiration (cm2 day− 1), Tp-input ismaximal potential transpiration of unstressed plants (cm2 day− 1)entered as input, and (Ra/Rp) is relative N uptake (-). This is based on theunderstanding that relative transpiration is a linear function of therelative growth rate (De Wit, 1958), which is a linear function of therelative N addition rate (Ingestad, 1982).
2.2.6. Comparison with measured data (data analysis)
The simulated results of water flow and solute transport were con­verted into units of mm and g m− 2, respectively, to facilitate the com­
parison to measured data. As transpiration and N leaching fluxes were
measured on a (near-)daily basis, the Root Mean Square Error (RMSE)
was calculated for these 2 data sets over time to assess the deviation
between simulated and measured data for each treatment:√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
∑n (O − M)2
RMSE =(5)in2.2.4. Time variable boundary conditions
The time-variable boundary conditions consisted of the measured
amounts of irrigation water with specified NO-3-N concentrations uniquefor every 18 treatments. We used the highest measured transpiration andN uptake (plot 1, treatment A90, I/T = 3.2) as potential transpiration
and N uptake for all treatments. As N uptake was only measured at theend of the season, potential N uptake was distributed throughout theseason at a fixed ratio in relation to measured transpiration. As totaltranspiration was measured in daily increments, potential transpiration
was fixed throughout 12 h of daylight and set to zero at nighttime. The
initial conditions were not precisely known since the plots had alreadyreceived fertigation for two weeks before transplanting. To eliminate theeffect of starting simulations with homogenous water and solute distri­butions in the root zone, water contents and solute concentrations 5 DAT
were inserted as initial conditions, as transpiration was still minimal,
while the decay reactions were already taking place.
The root zone density distribution was described according to a
formula developed by Vrugt et al. (2001):()()() − pz |z∗ − z|+ px |x∗ − x|
Xm
Zm
z
x
b(x, z) = 1 −
(3)
1−
e
Zm
Xm

where O and M are observed and simulated values and n the number ofdata pairs. When RMSE values of measurements with different units
were compared (i.e., transpiration (mm) and leached N (g m− 2)), they
were normalized (nRMSE) by dividing the RMSE value by the corre­
sponding maximal measured value of each data set.
To compare cumulative simulated and measured data across alltreatments, a linear reduction of the modeled data as a function of themeasured data was compared with the 1:1 line. The coefficient ofdetermination (R2) refers to the fit of the modelled data to the linear
reduction model.
2.2.7. Calibration
Some parameters could be estimated or calculated based on
measured data, others were taken from the literature, and some were
calibrated using the HYDRUS model. The calibration was carried out
using data from the plot that received the highest fertigation (plot 1,
treatment A90, I/T = 3.2, see Fig. 1). The objective function for thecalibration of optimized parameters was defined as an average of thenRMSE values for transpiration and leached N. The parameters weresequentially optimized by trial and error, preserving the value of pre­viously optimized parameters. This leaves the possibility that a change
in a subsequently optimized parameter could change the optimal valueof a previously optimized parameter.

where b is the spatial distribution of the root density (-) in two di­mensions (x,z), Zm and Xm are the maximal vertical and horizontalrooting distances from the base of the plant, respectively, which werethe limits of the lysimeter (cm), z* and x* are the depth and radius of the
maximal uptake intensity estimated at 5 cm each, and pz and px are
dimensionless shape parameters that describe the non-symmetry in root
water uptake in two dimensions (x,z). The root growth rate (gr(x,z)) was
estimated as 1 cm day− 1 in the x and z directions until the lysimeter’sedge was reached.
The root zone distribution and growth rate in the experiment were
unknown and were thus estimated. An increase in root growth due to N
deficiency observed in the measured data at 40 DAT (Groenveld et al.,2019) was not considered in these simulations.

2.2.8. Validation
After determining the optimal parameters for the data from the plotthat received the highest fertigation, the other 17 plots were run withthe same parameters to validate the function that reduces potential
transpiration as a function of N uptake (Eq. (4)). The linear regression of
the cumulative modeled data was compared with the 1:1 line of themeasured data. The RMSE of transpiration and leached N over time werecompared for each treatment.

2.2.5. Model modification
The HYDRUS (2D/3D) model (version 3.01) was modified to enable
solute transformations and active uptake simultaneously. Three solutes
are considered; the first solute acts like a tracer to represent salinity, thesecond solute represents N in the form of NO-3, which can be taken up by
the plant roots, drain out of the root zone, or decay into a daughter
product, which was the third solute that represents a gaseous denitrifi­
cation product (NO, N2O or N2). Nutrient uptake can be passive or
active, and denitrification is considered the first-order reaction with arate constant kD.
An additional formulation was implemented to consider the reduc­tion in potential transpiration as a function of N stress (considered as aratio of actual and potential cumulative N uptake until that time) thatwas calculated at every time step:

2.2.9. Calculation of active uptake
To calculate the contribution of active N uptake to total N uptake, the
HYDRUS (2D/3D) simulations were run with and without active uptake
enabled and stopped at the end of each daily time step. The difference inN uptake between the two simulations was considered the active uptake
contribution. The simulation was continued in this manner till the end of
the season, with the results of the ’active uptake enabled’ simulation forthe previous day used as the initial conditions for the subsequent day.
Cumulative active uptake was divided by cumulative total uptake andreported as a percentage for each treatment.
These values were compared with a percentage of active uptake
5T. Groenveld et al.

Agricultural Water Management 254 (2021) 106971(Uact) calculated according to a modified version of the Barber (1962)
calculation:
Uact =(Utot − Nav ∙Ttot )
∙100%
Utot

3. Results and discussion
3.1. Parameters estimated from experimental data

(6)

3.1.1. Hydraulic parameters
The water content as a function of the pressure head for perlite isshown with the standard deviation in Fig. 3. The hydraulic parameters
determined using the RETC 6.02 software (van Genuchten et al., 1991)
are reported in Table 2, and these were used to model the unsaturated
hydraulic conductivity as a function of the pressure head (van Gen­uchten, 1980) (Fig. 3). The retention curve and unsaturated hydraulicconductivity of sand and loam soils from the Carsel and Parrish catalog(1988) are plotted in the same graph for comparison. Besides ensuring awell-aerated root zone, the large pore volume of perlite gives it anadditional water holding capacity in the range of pressure heads most
relevant under frequent irrigation. The hydraulic conductivity of perlite
at saturation was 943 cm day− 1 and decreased rapidly with small in­creases in pressure head, similar to sand but less pronounced (Raviv
et al., 2002). Though the Richards equation was developed to describewater flow through unsaturated soil, it has also been used to model
water flow in perlite (Li and Babcock, 2015; Palla et al., 2009; Sandoval
et al., 2017).

where Utot is total nutrient uptake (g), Nav is the spatiotemporal average
nutrient concentration of the soil solution (g L− 1), and Ttot is totaltranspiration (L).
2.2.10. Sensitivity analysis
The sensitivity of all 18 simulations to individual parameters wasdemonstrated by running the simulations with all parameters fixed but
increasing and decreasing the analyzed parameter by 20%. The studied
parameters included z*, x*, pz, and px, which describe the root zone
density distribution (Eq. (3)), the first-order degradation rate constantkD, the Michaelis-Menten coefficient KM from Eq. (1), and the solute andwater uptake compensation parameters ΩC and ΩC. The results wereanalyzed by comparing the slope and intercept of the cumulative tran­
spiration and N leaching data to the 1:1 line and the average nRMSE oftranspiration and leached N.
2.2.11. Optimization
In fertigation of horticultural crops, the factors most easily controlled
by the grower are the amount of irrigation water, its fertilizer concen­
tration, and the irrigation frequency. As in the experimental setup, theirrigation frequency was a function of the amount of applied irrigationwater, distributed evenly throughout the day in pulses of 1 min duration
at a rate of 0.67 mm/min. The other environmental conditions, such as
potential transpiration and N uptake, the media type, and irrigationwater salinity, were the same as in the experiment described above. The
calibrated and validated model was used to find the optimal combina­
tion of N concentration in irrigation water and the total N applied. The
optimal combination should be close to the highest measured transpi­
ration and N uptake with the lowest amount of N lost to leaching or
gaseous emissions. Water and N use efficiencies (WUE and NUE) were
defined as the amounts of water transpired or N taken up divided by theamount applied. NUE is a measure that can be used to compare the N use
between different agricultural systems (van Groenigen et al., 2010).
However, since minimizing the N application maximizes the NUE, anobjective function in optimizing fertilizer application demands
including a condition or constraint that a certain percentage of the
maximal yield has to be achieved.

3.1.2. Solute transport parameters
The measured and simulated breakthrough curve data are shown in
Fig. 4. The average λL and λT estimated were 7.9 and 0.0001 cm,
respectively (Table 2). The ratio of the λT to the depth of the domain(7.9 cm/18 cm = 0.44) is much higher than that recommended for fieldscale experiments in the literature, which is usually only 0.1 (Gelhar,1986). However, Kim et al. (2002) found a ratio of 0.3 for a conservative
tracer in an unconfined sandy aquifer. Similarly, Boudreau et al. (2009)
reported an average ratio of 0.43 for a bromide tracer in pure sphagnum
peat media.
3.2. Modeling results
3.2.1. Calibration
The parameters optimized to fit measured data from the lysimeterwith the highest fertigation treatment are listed as “calibrated” inTable 2. Water compensation (ΩC) was not required and could thus be
kept at 1 (no compensation). Some N uptake compensation (ΩS) was
needed as nRMSE decreased for ΩS values between 1.0 and 0.9 but
remained the same for further decreases of ΩS. The first-order rate
parameter (kD) was tested for a range of values, and a distinct minimumof nRMSE was found at kD = 0.2. The N and water fluxes as a function oftime for the calibrated simulation are shown in Fig. 5. Simulated leachedFig. 3. The water content as a function of the pressure head for perlite, reported with the standard deviation, and the unsaturated hydraulic conductivity calculated
by means of van Genuchten (1980). Soil hydraulic functions for sand and loam using parameters from Carsel and Parrish (1988) are plotted for comparison.
6

T. Groenveld et al.

Agricultural Water Management 254 (2021) 106971

Table 2
Parameters describing hydraulic and transport properties of perlite, root distribution, stress response functions, and N uptake.
Parameter
Perlite media
ϴr
ϴs

α

n
l
Ks

ρb
λL
λT
Roots
Root distribution parameters
z*
x*
pz
px
gr (x,z)
Water stress response function parameters
P50
Pw3
PW
ΩC
Salinity stress response function parameters
EC50
ps3
φ
ΩS
Nitrate uptake and transformation
KM
kD
Cmax
Smin

Units

Value

SourceResidual water content
Saturated water content
Shape parameter
Shape parameter
Tortuosity parameter
Saturated hydraulic conductivity
Bulk density
Longitudinal dispersivity
Transverse dispersivity(cm3 cm− 3)(cm3 cm− 3)(cm− 1)(-)(-)(cm day− 1)(g cm− 3)
(cm)
(cm)

0.214
0.74
0.039
2.405
0.5
943
0.06
7.895
0.01

Depth of maximal uptake intensity
Radius of maximal uptake intensity
Vertical non-symmetry
Horizontal non-symmetry
Growth rate

(cm)
(cm)
(-)
(-)
(cm day− 1)

5
5
0.8
0.8
1

Tension for a 50% T reduction
Drought stress exponent
Permanent wilting point
Water uptake compensation

(cm)
(-)
(cm)
(-)

Salinity for a 50% T reduction
Salinity stress exponent
Osmotic coefficient
Solute uptake compensation(dS m− 1)
(-)
(-)
(-)

3
3
1
0.9

From literature
From literature
From literature
Calibrated

Michaelis-Menten coefficient, 50% of uptake
First-order rate constant
Maximal passive uptake
Minimal active uptake(mg cm− 3)(day− 1)
(mg cm− 3)(mg cm− 3)

0.001
0.2
10
0.00001

From literature
Calibrated
Estimated
Estimated

N was similar to measured values (RMSE = 0.166 g day− 1). Simulatedtranspiration was almost perfect (RMSE = 0.003 mm day− 1).

-200
3
-5000
1

Measured data
Measured data
Measured data
Measured data
Measured data
Measured data
Measured data
Measured data
Measured data
Estimated
Estimated
Estimated
Estimated
Estimated
Estimated
From literature
Estimated
Calibrated

The RMSE for transpiration and leached N are shown in Fig. 7C and
7D, respectively. They were low for all treatments and for both vari­ables. For transpiration, the RMSE increased slightly with increasing N
stress. The RMSE for leached N decreased linearly with a reduction inleached N, likely because the amounts of leached N became exceedingly
small.

3.2.2. Validation
Simulated cumulative transpiration and leached N are comparedwith the measured data in Fig. 6A and 6C, respectively. Without afunction to reduce potential transpiration due to a limitation in N uptake(Eq. (4)), simulated transpiration would remain 180 mm in all treat­ments and drop only slightly for the low I/T ratio, when transpiration
would be reduced due to drought and salinity stress (Fig. 1A). The slopeof the linear reduction of simulated transpiration is 0.82 (R2 = 0.976)(Fig. 6A), indicating that the N stress explains most but not all transpi­
ration reduction. The remaining reduction could still be due to the N
stress, which increases plant sensitivity to other stresses, such as nutri­tional imbalance, salinity, or drought. However, this hypothesis is notsupported by research combining salinity and N stress (Shenker et al.,2003). Increasing N stress led to a larger proportion of N being assigned
to the roots (Groenveld et al., 2019), which is biomass that is not
involved in transpiration. While it would be possible to include suchmechanisms into the model, this would require a more complex plantmodel to be implemented into the HYDRUS (2D/3D) software.
With a slope of 0.99 (R2 = 0.954), simulated N uptake matchedmeasured uptake well (Fig. 6B). N uptake in all simulations could not behigher than potential N uptake, taken from the highest fertigation
treatment. N fertigation was reduced incrementally so that the NUE
increased in lower treatments until the plant took up almost all applied
N. The amount of leached N was slightly underestimated for the plotsthat received higher rates of fertigation, getting the slope of 0.86 (R2= 0.976) (Fig. 6C). Simulated gaseous N emissions resulting fromdenitrification were lower than N missing from the N mass balance in thehigher N application range, with a slope of 0.46 (Fig. 6D). In the lower
fertigation range, the emissions were higher, and in general, the corre­
lation of the linear fit to the modeled data was lower (R2 = 0.46).

3.2.3. Active and passive uptake pathways
The percentages of active N uptake of total N uptake simulated by
HYDRUS (2D/3D) and calculated according to Barber (1962) (Eq. (6))
are shown in Fig. 8. When N is not a limiting factor (fertigatedN > 40 mg m− 2), passive uptake is overestimated by the Barber calcu­
lation, as was pointed out previously by Jungk (2002). This underesti­
mation of active uptake becomes more significant under deficit N
fertigation (< 40 mg m− 2), where the HYDRUS (2D/3D) simulations

Fig. 4. Relative (to initial salinities) measured irrigation and drainage watersalinities. Drainage water salinities were simulated using the average dis­persivity values.
7T. Groenveld et al.

Agricultural Water Management 254 (2021) 106971

Fig. 5. Simulated (lines) and measured (dots) daily N (A) and water (B) fluxes for the plot receiving the highest fertigation. The calibrated parameters are used. In A,
measured root N uptake over time is not shown as only total uptake at the end of the season was measured.

Fig. 6. Simulated vs. measured total transpiration (A), N uptake (B), N leaching (C), and N emission (D) for 18 fertigation treatments (distinguished by color, adarker hue indicates higher I/T). The fertigation NO-3 concentrations are listed in the legend. (For interpretation of the references to colour in this figure legend, thereader is referred to the web version of this article.)

was 1 mg L− 1, compared to a range of values from 0.3 m to 6 mg L− 1
found in the literature (Le Bot and Adamowicz, 2006). The 20% change
did not take the KM value out of this range, and, likely, this parameter
will only affect the simulation results when a larger change is imple­
mented. Lowering of the kD parameter resulted in reduced gaseous N
emissions and proportionally higher N leaching. This improved the fit of
the N leaching data but reduced the fit of the transpiration data. The
nRMSE was also reduced slightly, indicating that the kD parameter
optimized based on the entire data may not be suitable for all treat­ments. This may be due to changes in root zone aeration with changes in
irrigation regimes. Increasing the solute uptake compensation (reducingΩS) improved the fit of the transpiration data but decreased the fit ofleached N. The water uptake compensation factor (ΩC) could only be
reduced since the max value of 1 was used in all simulations (nocompensation). Doing so did not improve the fit of measured data andincreased the nRMSE (Table 3).

show that active uptake increases with increasing N stress, while acalculation according to Barber shows the opposite. The difference isdue to the consideration of spatiotemporal N distribution by HYDRUS(2D/3D) instead of using the single value for the soil solution N con­centration used by Barber.
3.2.4. Sensitivity analysis
A 20% change in the parameters describing the shape of the rootzone (z*, x*, px, pz) had only a very small effect on the results and did not
improve the fit of measured to modeled data (Table 3). The solute up­take compensation (ΩS) was 0.9 in these simulations, which could
explain the lack of sensitivity of leached N to the root zone shape pa­
rameters. The low sensitivity of the simulated results to the root zone
density distribution indicates that the lack of root data in this experi­mental setup is not critical. A 20% change in the KM parameter had little
impact on the simulated results. The KM value used in the simulations
8

T. Groenveld et al.

Agricultural Water Management 254 (2021) 106971

Fig. 7. RMSE for transpiration (A), and leached N (B) as a function of modeled transpiration and leached N fluxes. Different colors refer to different fertigation
concentrations, and a darker hue indicates higher I/T. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version ofthis article.)

3.0 g m− 2 of N lost as gaseous emissions, 6.1 g m− 2 of N leaching below
the root zone, the NUE of 0.73, and the WUE of 0.34.
The N fluxes, transpiration, NUE, and WUE for an N application
range of 30–50 g m− 2 are shown for the irrigation water N concentrationof 75 mg L− 1 and compared with those of 45 and 105 mg L− 1 (Fig. 8).
When an N application rate is below 40 g m− 2, higher irrigation water N
concentrations lead to more N uptake (Fig. 8A). The highest N concen­
tration did not produce maximum transpiration, as the reduced amountof irrigation water and reduced irrigation frequency limited transpira­tion (Fig. 8B). This reduction in transpiration also reduced active N
uptake according to Eq. (2). Plant N uptake was slightly lower for theirrigation water concentration of 105 mg L− 1 than 75 mg L− 1 at the N
application rate of 40 mg m− 2 (Fig. 8A). Under deficit N application(below 40 mg m− 2), transpiration remained the highest for the
75 mg L− 1 N concentration. The N emissions were higher for the higher
N application concentrations, whereas the amount of leached N was
higher for the lower N concentrations (Fig. 9CD). The N emissions can bein any form of gaseous N; the fraction of these gaseous forms that pose anenvironmental threat is unknown. The leached N is all in the form of
NO-3, and thus, its cost as a loss of N and environmental impacts can beevaluated. The NUE differs in the lower N application range but is thesame for all N concentrations when the N application is higher than
42 g m− 2 (Fig. 8E). There are, however, large differences in the WUE as
the lower N concentrations require more water to be applied to achieve
the same N application (Fig. 8F).

Fig. 8. Comparison of the percentage of N uptake by active pathway as
simulated by HYDRUS (2D/3D) and calculated according to the Barber method.
Different colors refer to different fertigation concentrations, and a darker hue
indicates higher I/T (see Fig. 1A). (For interpretation of the references to colourin this figure legend, the reader is referred to the web version of this article.)
Table 3The sensitivity of different variables to a 20% increase or decrease in selectedparameter values.
Calibrated values:
x*
–
+
z*
–
+
px
–
+
pz
–
+
KM
–
+
kD
–
+
ΩS
–
+
–
ΩC

T slopeT intercept

N slopeN intercept

nRMSE

0.816
0.002
-0.004-0.0020.000
0.000-0.002
0.000
0.000
0.000
0.000
0.000
-0.0020.016-0.036-0.01630.760
0.025
0.225
0.031
0.168
0.058
0.165
-0.0420.045-0.0070.000
0.406-0.057-0.449
1.881
4.179

0.862
-0.0070.008
0.004-0.004-0.0020.003
0.001-0.001
0.000
0.000
0.048
-0.041-0.0250.113-0.002-0.316-0.0390.042
0.015-0.016-0.0150.013
0.007-0.006-0.001-0.001
0.005-0.006-0.2680.051-0.001

0.191
0.000
0.001
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
-0.0020.003-0.0010.004
0.009

4. Conclusions
Reducing potential transpiration due to a limitation in N uptake in a
modified version of the HYDRUS (2D/3D) model could explain 82% of
the reduction in transpiration measured in a cucumber fertigation
experiment with 18 fertigation rates when an active uptake of NO-3 wasconsidered. Without the N stress reduction function (Eq. (4)), simulatedtranspiration would remain the same for all fertigation treatments(which was not the case in the experiments), demonstrating the neces­sity for such a function to model N fate and transport accurately. The
HYDRUS (2D/3D) simulations showed that the percentage of active N
uptake increased with N stress, in contrast to a calculation based on a
spatiotemporal averaged root zone N concentration. Increased spatio­
temporal resolution of the water and solute distribution in the root zone
by HYDRUS (2D/3D) allows for more precise quantification of theamount of solute transported by mass flow or diffusion. The sensitivityof simulated results to the spatial distribution of the roots in the small
domain studied was low, which indicates that the lack of measured root
distribution data is not critical. The model was then used to determine
the optimal N fertigation concentration and rate for cucumbers grown inperlite under given transpiration conditions and N demand over time.
Addressing the combined effect of drought and N stress on potential

The sensitivity is expressed as the change in the slope and intercept of thesimulated results to the measured data for transpiration (T) and N leached (N)(Fig. 6) for the parameters in Table 2. The reported nRMSE is the average nRMSEof the transpiration and leached N data for all 18 plots.

3.2.5. Optimization
The optimal N application under the described conditions, whileachieving 99.7% of the maximal yield, was determined to be 40 g m− 2with an irrigation water N concentration of 75 mg L− 1. This resulted in9T. Groenveld et al.

Agricultural Water Management 254 (2021) 106971

Fig. 9. Simulations of (A) N uptake, (B) transpiration, (C) N emission, (D) leached N, (E) NUE, and (F) WUE for different irrigation water N concentrations as afunction of the total amount of N applied.

transpiration and N uptake allows for more realistic N fertigation
modeling. This is essential when determining the optimal fertigation
strategy.

Funding
This research was partially supported by the Israel Ministry ofAgriculture and Rural Development (Eugene Kandel Knowledge Cen­ters) as part of the program “The root of the matter: The root zone
knowledge center for leveraging modern agriculture.”

Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.

References
Azad, N., Behmanesh, J., Rezaverdinejad, V., Abbasi, F., Navabian, M., 2020. An analysis
of optimal fertigation implications in different soils on reducing environmental
impacts of agricultural nitrate leaching. Sci. Rep. 10, 1–15. https://doi.org/
10.1038/s41598-020-64856-x.
Barber, S.A., 1962. A diffusion and mass flow concept of soil nutrient availability. Soil
Sci. 93, 39–49.
Boudreau, J., Caron, J., Elrick, D., Fortin, J., Gallichand, J., 2009. Solute transport in subirrigated peat-based growing media. Can. J. Soil Sci. 89. https://doi.org/10.4141/
CJSS08023.
Bristow, K.L., Šimůnek, J., Helalia, S.A., Siyal, A.A., 2020. Numerical simulations of theeffects furrow surface conditions and fertilizer locations have on plant nitrogen andwater use in furrow irrigated systems. Agric. Water Manag. 232. https://doi.org/
10.1016/j.agwat.2020.106044.

Acknowledgments
We would like to acknowledge Asaf Bokish from the Arava R&D
station for his contribution and support in thinking about the simulated
processes and Alon Ben-Gal for his comments on the manuscript. Author
contributions: T.G. and N.L. designed the research; T.G. performed theresearch and analyzed the data; A.A. and J.S. coded the model; T.G.
wrote the paper, and all authors contributed to the review of the paper.

10T. Groenveld et al.

Agricultural Water Management 254 (2021) 106971

Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24, 755–769.
Cerrato, M.E., Blackmer, A.M., 1990. Comparison of models for describing; corn yield
response to nitrogen fertilizer. Agron. J. 82, 138–143. https://doi.org/10.2134/
agronj1990.00021962008200010030x.
Dahan, O., Babad, A., Lazarovitch, N., Russak, E.E., Kurtzman, D., 2014. Nitrate leachingfrom intensive organic farms to groundwater. Hydrol. Earth Syst. Sci. 18, 333–341.
https://doi.org/10.5194/hess-18-333-2014.
De Wit, C.T., 1958. Transpiration and crop yields. Versl. Landbouw Onderz. Wagening.
Univ. Neth. 64, 59–84.
Epstein, E., Hagen, C.E., 1952. A kinetic study of the absorption of alkali cations byBarley Roots. Plant Physiol. 27, 457–474. https://doi.org/10.1104/pp.27.3.457.
Erisman, J.W., Sutton, M.A., Galloway, J., Klimont, Z., Winiwarter, W., 2008. How a
century of ammonia synthesis changed the world. Nat. Geosci. 1, 636–639. https://doi.org/10.1038/ngeo325.
Eshel, A., Beeckman, T. (Eds.), 2013. Plant Roots- The Hidden Half, Fourth ed. Taylor &Francis, Boca Raton, FL.
Gallardo, M., Padilla, F.M., Peña-Fleitas, M.T., de Souza, R., Rodríguez, A., Thompson, R.
B., 2020. Crop response of greenhouse soil-grown cucumber to total available N in a
Nitrate Vulnerable Zone. Eur. J. Agron. 114, 125993 https://doi.org/10.1016/j.
eja.2019.125993.
Gärdenäs, A.I., Hopmans, J.W., Hanson, B.R., Šimůnek, J., 2005. Two-dimensionalmodeling of nitrate leaching for various fertigation scenarios under micro-irrigation.
Agric. Water Manag. 74, 219–242. https://doi.org/10.1016/j.agwat.2004.11.011.
Gelhar, L.W., 1986. Hydrology from theory to applications. Water Resour. Res. 22,
135–145.
Głąb, T., Szewczyk, W., Gondek, K., Mierzwa-Hersztek, M., Palmowska, J., Nęcka, K.,
2020. Optimization of turfgrass fertigation rate and frequency. Agric. Water Manag.
234. https://doi.org/10.1016/j.agwat.2020.106107.
Groenveld, T., Ben-Gal, A., Yermiyahu, U., Lazarovitch, N., 2013. Weather determined
relative sensitivity of plants to salinity: quantification and simulation. Vadose Zone
J. 12. https://doi.org/10.2136/vzj2012.0180.
Groenveld, T., Kohn, Y.Y., Gross, A., Lazarovitch, N., 2019. Optimization of nitrogen use
efficiency by means of fertigation management in an integrated aquacultureagriculture system. J. Clean. Prod. 212, 401–408. https://doi.org/10.1016/j.
jclepro.2018.12.031.
Hartmann, A., Šimůnek, J., Aidoo, M.K., Seidel, S.J., Lazarovitch, N., 2017.
Implementation and application of a root growth module in HYDRUS. Vadose Zone
J. 1–16. https://doi.org/10.2136/vzj2017.02.0040.
He, H., Jansson, P.-E., Gärdenäs, A., 2020. CoupModel (v6.0): an ecosystem model forcoupled phosphorus, nitrogen and carbon dynamics – evaluated against empirical
data from a climatic and fertility gradient in Sweden. Geosci. Model Dev. Discuss.
1–55. https://doi.org/10.5194/gmd-2020-65.
He, W., Grant, B.B., Smith, W.N., VanderZaag, A.C., Piquette, S., Qian, B., Jing, Q.,
Rennie, T.J., Bélanger, G., Jégo, G., Deen, B., 2019. Assessing alfalfa production
under historical and future climate in eastern Canada: DNDC model development
and application. Environ. Model. Softw. 122. https://doi.org/10.1016/j.
envsoft.2019.104540.
Hopmans, J.W., Bristow, K.L., 2002. Current capabilities and future needs of root waterand nutrient uptake modeling. Adv. Agron. 77, 103–183. https://doi.org/10.1016/
S0065-2113(02)77014-4.
Houlton, B.Z., Almaraz, M., Aneja, V., Austin, A.T., Bai, E., Cassman, K.G., Compton, J.E.,
Davidson, E.A., Erisman, J.W., Galloway, J.N., Gu, B., Yao, G., Martinelli, L.A.,
Scow, K., Schlesinger, W.H., Tomich, T.P., Wang, C., Zhang, X., 2019. A world ofcobenefits: solving the global nitrogen challenge. Earth’s Future 7, 865–872. https://doi.org/10.1029/2019EF001222.
Ingestad, T., 1982. Relative addition rate and external concentration; Driving variables
used in plant nutrition research. Plant Cell Environ. 5, 443–453. https://doi.org/
10.1111/1365-3040.ep11611714.
Javaux, M., Schröder, T., Vanderborght, J., Vereecken, H., 2008. Use of a threedimensional detailed modeling approach for predicting root water uptake. Vadose
Zone J. 7, 1079–1088. https://doi.org/10.2136/vzj2007.0115.
Jungk, A., 2002. Dynamics of Nutrient Movement at the Soil–RootInterface. In:
Waisel, Y., Eshel, A., Kafkafi, U. (Eds.), Plant Roots:The Hidden half. Marcel Dekker,
New York, pp. 587–616. NY, USA,.
Kandelous, M.M., Šimůnek, J., van Genuchten, M.T., Malek, K., 2011. Soil water content
distributions between two emitters of a subsurface drip irrigation system. Soil Sci.
Soc. Am. J. 75, 488–497. https://doi.org/10.2136/sssaj2010.0181.
Karandish, F., Šimůnek, J., 2017. Two-dimensional modeling of nitrogen and water
dynamics for various N-managed water-saving irrigation strategies using HYDRUS.
Agric. Water Manag. 193, 174–190. https://doi.org/10.1016/j.agwat.2017.07.023.
Kim, D.J., Kim, J.S., Yun, S.T., Lee, S.H., 2002. Determination of longitudinal dispersivity
in an unconfined sandy aquifer. Hydrol. Process. 16, 1955–1964. https://doi.org/
10.1002/hyp.395.
Kim, W.S., Lieth, J.H., 2012. Simulation of year-round plant growth and nutrient uptake
in Rosa hybrida over flowering cycles. Hortic. Environ. Biotechnol. 53, 193–203.
https://doi.org/10.1007/s13580-012-0054-y.
Lazarovitch, N., Šimůnek, J., Shani, U., 2005. System-dependent boundary condition forwater flow from subsurface source. Soil Sci. Soc. Am. J. 69, 46. https://doi.org/
10.2136/sssaj2005.0046.
Le Bot, J., Adamowicz, S., 2006. Nitrogen nutrition and use in horticultural crops. J. Crop
Improv 15, 369–420. https://doi.org/10.1300/J411v15n02.
Levers, L.R., Schwabe, K.A., 2017. Biofuel as an integrated farm drainage management
crop: a bioeconomic analysis. Water Resour. Res. 53, 2940–2955. https://doi.org/
10.1002/2016WR019773.

Li, Y., Babcock, R.W., 2015. Modeling hydrologic performance of a green roof systemwith HYDRUS-2D. J. Environ. Eng. 141, 04015036 https://doi.org/10.1061/(ASCE)
EE.1943-7870.0000976.
Li, Y., Šimůnek, J., Zhang, Z., Jing, L., Ni, L., 2015. Evaluation of nitrogen balance in a
direct-seeded-rice field experiment using Hydrus-1D. Agric. Water Manag. 148,
213–222. https://doi.org/10.1016/j.agwat.2014.10.010.
Maas, E., Hoffman, G., 1977. Crop salt tolerance: evaluation of existing data. Water Irrig.
Proc. Int. Salin.
McMurtrie, R.E., Näsholm, T., 2018. Quantifying the contribution of mass flow tonitrogen acquisition by an individual plant root. New Phytol. 218, 119–130. https://doi.org/10.1111/nph.14927.
Padilla, F.M., Gallardo, M., Manzano-Agugliaro, F., 2018. Global trends in nitrateleaching research in the 1960–2017 period. Sci. Total Environ. 643, 400–413.
https://doi.org/10.1016/j.scitotenv.2018.06.215.
Palla, A., Gnecco, I., Lanza, L.G., 2009. Unsaturated 2D modelling of subsurface waterflow in the coarse-grained porous matrix of a green roof. J. Hydrol. 379, 193–204.
https://doi.org/10.1016/j.jhydrol.2009.10.008.
Pang, X.P., Letey, J., 1998. Development and evaluation of ENVIRO-GRO, an integrated
water, salinity, and nitrogen model. Soil Sci. Soc. Am. J. 62, 1418–1427.
Phogat, V., Skewes, M.A., Cox, J.W., Sanderson, G., Alam, J., Šimůnek, J., 2014. Seasonal
simulation of water, salinity and nitrate dynamics under drip irrigated mandarin(Citrus reticulata) and assessing management options for drainage and nitrateleaching. J. Hydrol. 513, 504–516. https://doi.org/10.1016/j.jhydrol.2014.04.008.
Radcliffe, D., Šimůnek, J., 2010. Soil Physics with HYDRUS - Modeling and Applications.
CRC Press.
Raij, I., Šimůnek, J., Ben-Gal, A., Lazarovitch, N., 2016. Water flow and multicomponent
solute transport in drip- irrigated lysimeters. Water Resour. Res. 52, 6557–6574.
https://doi.org/10.1029/2008WR006912.M.
Ramos, T.B., Šimůnek, J., Goncalves, M.C., Martins, J.C., Prazeres, A., Pereira, L.S., 2012.
Two-dimensional modeling of water and nitrogen fate from sweet sorghum irrigatedwith fresh and blended saline waters. Agric. Water Manag. 111, 87–104. https://doi.
org/10.1016/j.agwat.2012.05.007.
Raviv, M., Wallach, R., Silber, A., Bar-Tal., A., 2002. Substrates and their analysis. In:
Savvas, D., Passam, H. (Eds.), Hydroponic production of vegetables andornamentals. Em- bryo Publ., Athens, pp. 25–101.
Rengel, Z., 1993. Mechanistic simulation models of nutrient uptake: a review. Plant Soil
152, 161–173. https://doi.org/10.3950/jibiinkoka.94.12_1844.
Rewald, B., Ephrath, J.E., Rachmilevitch, S., 2011. A root is a root is a root? Water
uptake rates of Citrus root orders. Plant Cell Environ. 34, 33–42. https://doi.org/
10.1111/j.1365-3040.2010.02223.x.
Rezaei, M., De Pue, J., Seuntjens, P., Joris, I., Cornelis, W., 2017. Quasi 3D modelling ofvadose zone soil-water flow for optimizing irrigation strategies: challenges,
uncertainties and efficiencies. Environ. Model. Softw. 93, 59–77. https://doi.org/
10.1016/j.envsoft.2017.03.008.
Roose, T., Schnepf, A., 2008. Mathematical models of plant-soil interaction. Philos.
Trans. A Math. Phys. Eng. Sci. 366, 4597–4611. https://doi.org/10.1098/
rsta.2008.0198.
Rubio, G., Lynch, J.P., 2007. Compensation among root classes in Phaseolus vulgaris L.
Plant Soil 290, 307–321. https://doi.org/10.1007/s11104-006-9163-7.
Russo, D., Laufer, A., Shapira, R.H., Kurtzman, D., 2013. Assessment of solute fluxes
beneath an orchard irrigated with treated sewage water: a numerical study. Water
Resour. Res. 49, 657–674. https://doi.org/10.1002/wrcr.20085.
Sandoval, V., Bonilla, C.A., Gironás, J., Vera, S., Victorero, F., Bustamante, W., Rojas, V.,
Leiva, E., Pastén, P., Suárez, F., 2017. Porous media characterization to simulate
water and heat transport through green roof substrates. Vadose Zone J. 16, 0.
https://doi.org/10.2136/vzj2016.10.0101.
Sands, P., Smethurst, P., 1995. Modelling nitrogen uptake in ingestad units using
Michaelis-Menten kinetics. Aust. J. Plant Physiol. 22, 823. https://doi.org/10.1071/
PP9950823.
Schröder, N., Lazarovitch, N., Vanderborght, J., Vereecken, H., Javaux, M., 2013.
Linking transpiration reduction to rhizosphere salinity using a 3D coupled soil-plantmodel. Plant Soil 377, 277–293. https://doi.org/10.1007/s11104-013-1990-8.
Šejna, M., Šimůnek, J., van Genuchten, M., 2018. HYDRUS User Manual, Version 3. ed,
PC Process. 〈https://doi.org/10.1007/springerreference_28001〉.
Shenker, M., Ben-Gal, A., Shani, U., 2003. Sweet corn response to combined nitrogen andsalinity environmental stresses. Plant Soil 256, 139–147. https://doi.org/10.1023/
A:1026274015858.
Shouse, P.J., Ayars, J.E., ŠimŮnek, J., 2011. Simulating root water uptake from a shallow
saline groundwater resource. Agric. Water Manag. 98, 784–790. https://doi.org/
10.1016/j.agwat.2010.08.016.
Šimůnek, J., Hopmans, J.W., 2009. Modeling compensated root water and nutrient
uptake. Ecol. Modell. 220, 505–521. https://doi.org/10.1016/j.
ecolmodel.2008.11.004.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2018b. HYDRUS3D Technical Manual,〈https://www.pc-progress.com/en/Default.aspx?hydrus-3d#k3〉.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J. 15,
vzj2016.04.0033 https://doi.org/10.2136/vzj2016.04.0033.
Šimůnek, J., Šejna, M., Van Genuchten, M.T., 2018a. New features of version 3 of the
HYDRUS (2D/3D) computer software package. J. Hydrol. Hydromech. 66, 133–142.
https://doi.org/10.1515/johh-2017-0050.
Skaggs, T.H., Shouse, P.J., Poss, J.A., 2006. Irrigating forage crops with saline waters: 2.
Modeling root uptake and drainage. Vadose Zone J. 5, 824–837. https://doi.org/
10.2136/vzj2005.0120.
Snow, V.O., Cichota, R., McAuliffe, R.J., Hutchings, N.J., Vejlin, J., 2017. Increasing the
spatial scale of process-based agricultural systems models by representing

11

T. Groenveld et al.

Agricultural Water Management 254 (2021) 106971

heterogeneity: the case of urine patches in grazed pastures. Environ. Model. Softw.
90, 89–106. https://doi.org/10.1016/j.envsoft.2017.01.005.
Tao, Y., Li, N., Wang, S., Chen, H., Guan, X., Ji, M., 2021. Simulation study on
performance of nitrogen loss of an improved subsurface drainage system for onetime drainage using HYDRUS-2D. Agric. Water Manag. 246, 1–16. https://doi.org/
10.1016/j.agwat.2020.106698.
Thompson, R.B., Incrocci, L., van Ruijven, J., Massa, D., 2020. Reducing contaminationof water bodies from European vegetable production systems. Agric. Water Manag.
240. https://doi.org/10.1016/j.agwat.2020.106258.
Tischner, R., Kaiser, W., 2007. Nitrate Assimilation in Plants, Biology of the Nitrogen
Cycle. Elsevier B.V. https://doi.org/10.1016/B978-044452857-5.50019-9
Tripler, E., Shani, U., Ben-gal, A., Mualem, Y., 2012. Apparent steady state conditions inhigh resolution weighing-drainage lysimeters containing date palms grown under
different salinities. Agric. Water Manag. 107, 66–73. https://doi.org/10.1016/j.
agwat.2012.01.010.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898. https://doi.org/
10.2136/sssaj1980.03615995004400050002x.
van Genuchten, M.T., 1987. A numerical model for water and solute movement in andbelow the root zone, Unpublished Research Report.
van Genuchten, M.T., Hoffman, G.J., 1984. Analysis of crop production. In: Shainberg, I.,
Shalhevet, J. (Eds.), Soil Salinity Under Irrigation. Springer, New York, pp. 258–271.
van Genuchten, M.T., Leij, F.J., Yates, S.R., 1991. The RETC code for quantifying thehydraulic functions of unsaturated soils. United States Environ. Reseach Lab, p. 93.
https://doi.org/10.1002/9781118616871.
van Groenigen, J.W., Velthof, G.L., Oenema, O., Van Groenigen, K.J., Van Kessel, C.,
2010. Towards an agronomic assessment of N2O emissions: a case study for arable

crops. Eur. J. Soil Sci. 61, 903–913. https://doi.org/10.1111/j.13652389.2009.01217.x.
Vrugt, J.A., Hopmans, J.W., Šimůnek, J., 2001. Calibration of a two-dimensional root
water uptake model. Soil Sci. Soc. Am. J. 65, 1027–1037. https://doi.org/10.2136/
sssaj2001.6541027x.
Wang, J., Baerenklau, K.A., 2014. Crop response functions integrating water, nitrogen,and salinity. Agric. Water Manag. 139, 17–30. https://doi.org/10.1016/j.
agwat.2014.03.009.
Wang, Y., Ying, H., Yin, Y., Zheng, H., Cui, Z., 2019. Estimating soil nitrate leaching of
nitrogen fertilizer from global meta-analysis. Sci. Total Environ. 657, 96–102.
https://doi.org/10.1016/j.scitotenv.2018.12.029.
Wang, Z., Li, J., Li, Y., 2014. Simulation of nitrate leaching under varying drip system
uniformities and precipitation patterns during the growing season of maize in the
North China Plain. Agric. Water Manag. 142, 19–28. https://doi.org/10.1016/j.
agwat.2014.04.013.
Wever, G., Nowak, J.S., Oliveira, O.M.D.S., Winkel, A. Van, 2004. Determination of
hydraulic conductivity in growing media. Acta Hortic. 648, 135–143.
Wolf, J., 2012. LINGRA-N a grassland model for potential, water limited and N limited
conditions (FORTRAN). Wageningen University, Wageningen,, The Netherlands,
p. 65.
Yang, R., Tong, J., Hu, B.X., Li, J., Wei, W., 2017. Simulating water and nitrogen loss
from an irrigated paddy field under continuously flooded condition with Hydrus-1Dmodel. Environ. Sci. Pollut. Res. 24, 15089–15106. https://doi.org/10.1007/
s11356-017-9142-y.
Zhu, Y., Ren, L., Skaggs, T.H., Lu, H., Yu, Z., Wu, Y., Fang, X., 2009. Simulation ofPopulus euphratica root uptake of groundwater in an arid woodland of the Ejina
Basin, China. Hydrol. Process. 23, 2460–2469. https://doi.org/10.1002/hyp.

12