Agricultural Water Management 258 (2021) 107208

Contents lists available at ScienceDirect

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

Assessing the nitrate vulnerability of shallow aquifers under Mediterranean
climate conditions
Paolo Nasta a, *, Giuliano Bonanomi b, c, Jirka Šimůnek d, Nunzio Romano a, eaDepartment of Agricultural Sciences, AFBE Division, University of Naples Federico II, Portici, Naples, Italy
Department of Agricultural Sciences, BPSAF Division, University of Naples Federico II, Portici, Naples, Italyc
Task Force on Microbiome Studies, University of Naples Federico II, Naples, Italyd
Department of Environmental Sciences, University of California Riverside, CA, USAe
The Interdepartmental Center for Environmental Research (C.I.R.AM.), University of Naples Federico II, Naples, Italyb

A R T I C L E I N F O

A B S T R A C T

Handling Editor: Dr. Z. Xiying

The EU Nitrates Directive calls for urgent integration of process-oriented indicators of nitrate fate with map
overlay approaches for assessing nitrate vulnerable zones (NVZs). In the region of Campania (southern Italy),groundwater contamination represents a serious concern because of the presence of intensive agricultural
practices and livestock farming. A protocol was proposed to assess the probability distribution of the following
three indicators of nitrate transport across the vadose zone: i) annual cumulative nitrate flux entering the shallow
aquifer, ii) annual cumulative root nitrate uptake, and iii) nitrate transit time across the vadose zone. This method
involves numerical simulations of soil water flow and solute transport using Hydrus-1D for a representative 10m-thick soil profile beneath an irrigated maize plot located in a study area within the Sele plain, Campania. Two
scenarios are built by running a set of one hundred (20-yr-long) simulations at a daily time resolution: i) 195 kg
N ha-1 of nitrate fertilizer (urea) is applied annually to estimate the probability distributions of annual cumu­
lative nitrate leaching and root nitrate uptake (Scenario 1), ii) 195 kg N ha-1 of nitrate fertilizer (urea) is applied
once to estimate the nitrate transit time distribution across the soil profile (Scenario 2). In each Scenario, thesimulations consider two agricultural practices with either one or three annual nitrate fertilizer applications
subject to randomly generated climate forcing using a Monte Carlo approach. Variations in soil and rainfall
properties are described by the Miller-Miller geometric similitude and the Poisson parameterization,
respectively.
In Scenario 1, roots absorb on average 80.3 kg N ha-1 every year (corresponding to 40% of applied nitrate),and the median annual cumulative nitrate flux across the soil profile bottom is 74.9 kg N ha-1 (representing about
38% of applied nitrate) when urea is applied in a single treatment. In contrast, when fertilizer is applied in threetreatments, 112.4 kg N ha-1 (corresponding to almost 60% of applied nitrate) of nitrate is removed by root water
uptake, and the median annual cumulative nitrate leaching is 52.4 kg N ha-1 (corresponding to about 27% of
applied nitrate). In Scenario 2, mean transit time values are 2741 days, 2707 days, and 2650 days when urea isapplied on April 1st, June 1st, and August 1st, respectively. Our model simulations provide useful indicators ofnitrate transport and can be integrated with map overlay procedures for delineating nitrate vulnerable zones.

Keywords:
Groundwater
Nitrate leaching
Root uptake
Transit time distribution (TTD)
Hydrus-1D
HYPRES database
Campania

1. Introduction
The agri-food system is one of the main socio-economic sectors in
Campania, a region of southern Italy where the most productive areas
extend over the internal and coastal alluvial plains with rather intensive
agricultural practices. Campania is also a renowned region worldwide
for the breeding of buffaloes to produce mozzarella cheese and otherexcellent dairy products. Extensive fertilizers and animal manure ap­
plications boost crop production at the cost of polluting the environment(Infascelli et al., 2009). The rapid deterioration of groundwater quality
observed throughout the world over the last decades represents a serious
concern, particularly in Campania, where great attention is devoted tothe transport toward shallow aquifers of potentially toxic elements(PTEs) that are detrimental to human health. For this reason, European* Corresponding author.
E-mail address: paolo.nasta@unina.it (P. Nasta).
https://doi.org/10.1016/j.agwat.2021.107208
Received 29 April 2021; Received in revised form 21 July 2021; Accepted 22 September 2021
Available online 5 October 2021
0378-3774/© 2021 Elsevier B.V. All rights reserved.

P. Nasta et al.

Agricultural Water Management 258 (2021) 107208

Union (EU) policies aim at protecting shallow aquifers by promoting
sustainable land management, whose efficiency is entrusted to a reliablemapping of groundwater vulnerability (Mylevaganam and Ray, 2016;
Wall et al., 2020). Vulnerability can be conceptualized in different ways,
but in this study, we follow the definitions given by National ResearchCouncil (1993) that groundwater vulnerability is viewed as the easewith which a PTE applied on (or near) the land surface can migrate to
the aquifer under a given set of agricultural agronomic practices,
contaminant characteristics, and hydrogeological sensitivity conditions.
In most cases, and especially with the help of Geographic Informa­tion System (GIS) tools, widespread approaches are based on mapping
topographic, hydrogeological, pedological, land-use features controlling
the propagation of PTEs through the unsaturated and saturated zones.
Scores and weights are assigned to the various environmental factors to
obtain classes of groundwater vulnerability to contamination from acertain source following a map overlay procedure (Doerfliger et al.,
1999; Rupert, 2001; McLay et al., 2001; Antonakos and Lambrakis,2007). The parameter weighting and rating approaches (also known as
Point Count System Models; e.g., DRASTIC, SINTACS) can be easily
employed in large-scale applications, but express groundwater vulner­
ability in a rather empirical way and show the effect of a certain degreeof subjectivity (Van Stempvoort et al., 1993; McLay et al., 2001; Pan­
agopoulos et al., 2006; Huan et al., 2012). In contrast, the use of aprocess-oriented model, simulating water flow and solute transport in
the porous media, enables effective indicators (sometimes also referred
to as “functional performance indicators”) of PTE transport to be iden­tified. This represents a step forward and a more objective and efficientmethod to assess groundwater vulnerability, especially for shallow
aquifers. However, the prediction performance of these types of models
is highly dependent on the quantity and quality of information on input
data and system properties. It seems still rather unfeasible to assess thegroundwater vulnerability over relatively large spatial scales by using
three-dimensional (3D) process-oriented hydrological models. However,a viable way to provide details about solute transport is to implement
one-dimensional (1D) models in spatially representative soil profiles of
the area of interest.
The Nitrates Directive 91/676/EEC (reiterated by the new 2006/
118/EC) promoted by the European Council aims at preventing andreversing groundwater degradation induced by nitrates due to agricul­
tural practices. This directive requires the member states to identify
Nitrate Vulnerable Zones (NVZs), allowing, among other things, for a
threshold nitrate concentration value of 50 mg L-1 in the water bodies.
There is a vigorous debate in Campania because of the update of the
NVZs delimitation resulting from the Regional Council resolution n.
762/2017 and the new regional action plan for nitrates vulnerable areasfrom agricultural sources (Busico et al., 2019). A one-off Action Plan of
100 million euros is about to start reducing the livestock load tovulnerable areas that, since March 2020, have to be subjected to a
spreading limit of 170 kg of nitrogen per hectare per year. Recent studies
have reported that over the four-year period of 2012–2015, nitrateconcentrations in groundwater reached values higher than the drinking
water threshold (50 mg L-1 in Italy), with an increasing trend, in com­
parison to the previous four-year period (2008–2011) (Ducci et al.,2019).
As also pointed out by Gómez-Hernández and Wen (1994), the mainindicators controlling the vulnerability of shallow groundwater are ob­tained in the present study by simulating the nitrate fate across a
10-m-thick layered soil profile using Hydrus-1D (Mattern and Van­clooster, 2010; Šimůnek et al., 2016; Sprenger et al., 2016b; Mokari
et al., 2019). The goal is to build plausible scenario-based model simu­
lations to investigate the fate of nitrate, within a stochastic framework,
under random-generated Mediterranean climate forcing conditions. Ourinvestigations refer to the Sele River alluvial plain in Campania, arenowned agri-food district for agricultural produce and buffalo dairy
products. The probability distributions of three fundamental indicators
are assessed to characterize the nitrate transport: i) annual nitrateleaching, ii) annual root nitrate uptake, iii) vadose zone nitrate transit
time. One hundred 20-yr-long simulations were run for each scenario
using Hydrus-1D:1) Scenario 1 assumes 195 kg N ha-1 of nitrate fertilizer (urea) appliedannually to estimate the probability distributions of annual cumu­
lative nitrate leaching and root nitrate uptake;2) Scenario 2, instead, assumes 195 kg N ha-1 of nitrate fertilizer (urea)
applied only once to estimate the nitrate transit time distribution
across the soil profile.
2. Materials and methods
2.1. Site description and available data set
The Sele alluvial plain covers an area of about 500 km2 that spansaround the Sele River near the city of Salerno, located in the region ofCampania, southern Italy. The climate is typically Mediterranean, withcold rainy winters and hot dry summers. 68%, 23%, and 8% of the
agricultural land of the Sele River Plain is used for crops, fruit orchards,
and grazing pastures, respectively (http://www.agricoltura.regione.
campania.it/pubblicazioni/territorio_rurale.html).
The predominant crops are vegetables (8065 ha), including fennel,
lettuce, sweet pepper, tomatoes, rocket salad, broccoli, cauliflower,
green beans, eggplants, cabbage, and several others. Forages coverabout 7481 ha, and maize is the predominant forage, subject to irriga­
tion and fertilization practices. Considering the excessive fragmentationof vegetables, we assumed the maize as the predominant land cover inthe study’s area, and for this reason, we consider maize in model
simulations.
The hydrogeology of the Sele Plain belongs to the PlioceneQuaternary hydrogeological complex with porous aquifers formed by
epiclastic deposits of alluvial and coastal plains (De Vita et al., 2018;
Tufano et al., 2020). Groundwater fluxes are triggered by a hydraulicgradient of 0.6% in the east-west direction from the mountainous region
situated on the eastern border of the study area towards the coastline.
The shallow aquifer permeability ranges between 1.6 × 10-3 m s-1 and
2.0 × 10-3 m s-1, transmissivity spans between 3.0 × 10-2 m2 s-1 and 4.8
× 10-2 m2 s-1, and storage ranges between 0.3 × 10-4 and 1.5 × 10-4
(Cassa per il Mezzogiorno, 1983).
Among multiple nitrate sources, crop management based on nitrogen
fertilizers (mainly urea and ammonium sulfate) and animal manure
spreading represent the primary threat to groundwater quality (Bon­fante et al., 2019). The Campania Regional Council has recently pro­vided a new delimitation of nitrate vulnerable zones (NVZs) influenced
by agricultural practices. Maps of NVZs (green areas in Fig. 1b) can be
downloaded from the geoportal of the Campania region (https://sit2.
regione.campania.it/content/zone-vulnerabili-ai-nitrati in Italian). To
meet this study’s aim, we considered a portion of the Sele alluvial plain
with available information on depths to the groundwater table (222.72
km2) within the NVZ (Fig. 1d). A geological study based on 11 coredrilling lithostratigraphy surveys (http://www.difesa.suolo.regione.
campania.it/content/view/147/) reports the main characteristics ofthe geological setting. The vadose zone is characterized by silty-clay
alluvial deposits. According to the USDA soil classification, we reporta predominance of Mollic HAPLOXERERTS, Typic HAPLOXERERTS, andTypic CALCIXERERTS.
Soil physical and chemical properties, namely percentages of sand
(Sa), silt (Si), and clay (Cl), soil oven-dry bulk density (ρb), and organiccarbon (OC), are available from a recent study that was conducted on the
farmland in Campania within the Campania Trasparente project (Nastaet al., 2020b). Soil samples were collected on the soil surface (topsoil
samples) and at soil depths spanning between 80 cm and 110 cm (sub­soil samples) to determine soil physical and chemical properties. A totalof 3316 and 500 soil samples in the topsoil and subsoil, respectively,
were collected in Campania (Nasta et al., 2020a). In the study area of the2

P. Nasta et al.

Agricultural Water Management 258 (2021) 107208

Fig. 1. a) Digital Elevation Model (DEM) of Campania Region, b) nitrate vulnerable zones, NVZs (green area) in Campania (https://sit2.regione.campania.it/conten
t/zone-vulnerabili-ai-nitrati) and the study area (red area) in the Sele alluvial Plain, c) DEM of the study area with 86 topsoil samples (black circles), and d) the depth
to groundwater in the study area. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Sele alluvial plain, a total of 86 topsoil (representing the uppermost
15 cm of soil profile) and 50 subsoil (in the range from 90 cm to 120 cm
soil depths) undisturbed samples were collected (Fig. 1c) (Table 1).
These soil samples were used to map the values of sand, silt, and clay
content, oven-dry soil bulk density, and soil organic carbon (Fig. S1 in
Supplemental Material). The OC values were not measured for the
subsoil cores.
A weather station is located near the city of Battipaglia (33 N
497,979 m, 4,495,588 m, 64 m a.s.l.), where daily values of rainfall (R)and daily minimum and maximum temperature data are available from1999 to 2018 (20-yr long time series). In the recording period, the meanannual rainfall is 89.0 cm, and the mean annual temperature is 17.1 ◦ C.
Daily values of grass-reference potential evapotranspiration (ETo) were
calculated using the temperature-based method proposed by Hargreaves(Hargreaves and Allen, 2003). This relation only requires the minimum
and maximum temperatures, while the extraterrestrial radiation term isestimated using the study site’s latitude and the day of the year. For theSele River plain, Pelosi et al. (2016) have demonstrated that the Har­
greaves equation’s performance proved to be comparable to that of thewell-known physically-based Penman-Monteith equation (Allen et al.,1998). The mean annual ET0 is 111.5 cm. Assuming a uniform land
cover of maize (Zea mays L.), crop-specific potential evapotranspiration(ETp) under standard conditions and without water limitations is
calculated by multiplying ET0 by the specific crop (maize) coefficient, Kc(i.e., ETp = Kc×ET0). The leaf area index (LAI) is used to partition ETp
into potential evaporation (Ep) and potential transpiration (Tp) using thefollowing equation (Ritchie, 1972):

Table 1
Descriptive statistics for oven-dry bulk density (ρb), percentages of sand (Sa), silt(Si), and clay (Cl), and soil organic carbon (OC) content for 86 topsoil and 50
subsoil samples in the experimental area in the Sele alluvial Plain (from Nastaet al., 2020b).
ρb

Sand
%

Silt
%

Clay
%

OC
%

1.37
0.15
11.1
1.08
1.70
1.23
0.20
16.1
0.12
1.44

28.21
10.56
37.4
3.35
64.10
24.86
14.70
59.1
3.00
78.12

38.43
8.94
23.3
18.30
63.45
35.72
10.39
29.1
13.09
55.45

33.36
9.96
29.9
13.22
60.88
39.43
12.98
32.9
8.71
67.19

1.71
1.04
61.1
0.054
5.59

g/
cm3
Topsoil

Subsoil

Mean
Standard deviation
Coefficient of variation
Minimum
Maximum
MeanStandard deviation
Coefficient of variation
Minimum
Maximum

EP = ET p e− 0.463LAI(1)

Potential transpiration, Tp (corresponding to potential root water
uptake), is obtained by subtracting potential evaporation (Eq.1) from
ETp. The crop coefficient (Kc; dimensionless), the leaf area index (LAI;
dimensionless), and the maximum root depth for maize are timevariable and change daily (more details are reported in the Supple­
mentary material) (Nasta and Gates, 2013; Bonfante et al., 2019; Corbari
et al., 2021).
3

P. Nasta et al.

Agricultural Water Management 258 (2021) 107208

and θr (cm3 cm-3) and θs (cm3 cm-3) are residual and saturated soil watercontents, respectively. The two parameters m and n are related as fol­lows: m = 1–1/n (van Genuchten, 1980).
Considering the degree of saturation, Se = (θ− θr)/(θs− θr), which
varies from 0 (θ = θr) to 1 (θ = θs), the unsaturated hydraulic conduc­tivity function, K(Se) is given by the following equation (van Genuchten,1980):[)m ]2(
K(Se ) = Ks Sτe 1 − 1 − S1/m
(5)e

By referring to the daily data recorded at the Battipaglia weather
station from 1999 to 2018, the annual mean sum of ETp is 57.6 cm,
representing a reduction of about 50% of ET0 (see Fig. S3 in Supple­
mentary material and Table 2). Annual mean sums of Tp and Ep are
40.4 cm and 17.2 cm, respectively, representing about 70% and 30% ofETp, respectively.
Rainfall interception (IR) was calculated using the following formula(Braden, 1985):()
1
IR = aLAI 1 −(2)
1 + bP/aLAI

where Ks (cm d-1) is the saturated hydraulic conductivity and τ (-) is atortuosity parameter that is assumed to be τ = 0.5 (Mualem, 1976).
Given that the mode of the statistical distribution of depths to thewater table is about 10 m within the study area in the Sele alluvial plain(not shown), we set the thickness of the soil profile to be 10 m. The lower
boundary condition is set to a constant pressure head (ψ = 0 cm),
expressing a constant level of the water table of the shallow aquifer. The
upper boundary condition depends on fluxes occurring at the soil surface(Rnet and Ep).
Plant potential transpiration, Tp, determines potential root water
uptake (a sink term in the Richards equation) in the soil profile. The root
distribution is assumed to be linear in the soil profile by varying from amaximum at the soil surface to a minimum at the maximum rooting
depth, which is considered time-variant (see Fig. 3c). Both Ep and Tp are
reduced by water limitation and stresses to actual evaporation (Ea) andactual transpiration (Ta). The actual root water extraction rate, ξ, ismodeled using the method proposed by Feddes et al. (1978). Maize Tp isreduced between ψ = − 500 cm and the wilting point that corresponds
to ψ wp = − 8000 cm. To mitigate excessive root water uptake stress ofmaize during the dry growing season, Hydrus-1D has an option to trigger
irrigation when a user-specified pressure head is reached in a selected
observation node. Irrigation starts after a user-specified lag period at auser-specified irrigation rate. The duration of irrigation is also specified.
When surficial matric pressure head, ψ surf, is less than a thresholdvalue, ψ min, the prescribed flux boundary condition (a Neumann con­dition) specified at the upper boundary of the soil profile automaticallyswitches to a constant pressure head boundary condition (a Dirichletcondition), resulting in a reduction of potential evaporation (Ep) into
actual evaporation (Ea). A proportion of net rainfall, Rnet, is turned into
runoff when ψ surf exceeds the maximum soil surface matric potential
value, ψ max, that represents the nominal depth of surface water ponding
allowable before runoff generation.
Hydrus-1D simultaneously considers the fate and transport of mul­
tiple solutes subjected to first-order degradation reactions. Hanson et al.
(2006) considered a urea–ammonium–nitrate fertigation system subject
to nitrogen (N) transformation processes. Urea, (NH2)2CO, degrades toammonium, NH4+, that is subsequently transformed into nitrate, NO3-,
by the process of nitrification. NO3- is then subject to denitrification. The
partial differential equations governing the one-dimensional transportof N involved in sequential first-order decay chain reactions duringtransient water flow in a variably-saturated rigid porous medium aretaken as:()∂θC1 ∂
∂C1
∂qC1
′
(6)
=
θDw1
−
− μw,1 θC1
∂z
∂t
∂z
∂z

where a (cm d-1) is assumed equal to 0.025 cm d-1 and b = 1− e-0.463LAI.
Interception is subtracted from measured rainfall to obtain net rainfall(Rnet). Descriptive statistics of annual mean sums of the aforementioned
weather components are reported in Table 2.
Farmers in the Sele River plain mostly apply urea as fertilizer. Urea isincorporated into the soil by rototilling to reduce volatilization. We
explored two agricultural practices when urea is applied with a single
treatment (195 kg N ha-1) on April 1st or a sequence of three treatments(65 kg N ha-1) on April 1st, June 1st, and August 1st (Bonfante et al.,2019). Heterotrophic bacteria hydrolyze urea (NH2)2CO to form
ammonium (NH4+), which is sequentially nitrified by autotrophic bac­teria to nitrite (NO2-) and nitrate (NO3-) that, in turn, is subsequently
denitrified to form dinitrogen (N2 and N2O). Since anaerobic conditions
rarely occur in the deep vadose zone profile, we neglect the denitrifi­
cation process (transformation of NO3- into N2 and N2O) (Hanson et al.,2006; Ramos et al., 2012; Bradshaw et al., 2013; Deb et al., 2016; Busico
et al., 2020). Therefore, the following three solute species are consideredin this study: (NH2)2CO → NH4+ → NO3-.
2.2. Soil water balance and nitrate transport simulation with Hydrus-1DThe water balance in the soil-plant-atmosphere system was evalu­
ated using Hydrus-1D (Šimůnek et al., 2016), which numerically solves
the one-dimensional Richards equation:
(
[
)]
∂θ ∂
∂ψ
(3)
=
K(ψ )
+1
− ξ(z, ψ , Tp )∂z∂t ∂z
where t is time, expressed in units of days (d), ψ is the matric pressure
head (cm), z is soil depth (positive upward) (cm), θ (cm3 cm-3) is the soil
water content, and ξ (z, ψ , Tp) is the sink term (d-1) describing the actual
plant root water extraction rate function depending on z, ψ , and po­
tential transpiration (Tp). Hydrus-1D numerically solves the partialdifferential Eq. (3) by using a finite element scheme for the spatial dis­cretization and a finite difference scheme for the time discretization. The
soil water retention function θ(ψ ) is described by van Genuchten’s
equation (van Genuchten, 1980):
θ(ψ ) = θr +θs − θr[1 + (α|ψ |n ) ]m(4)

where α (cm-1), m (-), and n (-) are water retention shape parameters,
Table 2
Descriptive statistics of annual sums of gross rainfall (R), rainfall interception(IR), net rainfall (Rnet), grass-reference potential evapotranspiration (ET0), maize
potential evapotranspiration (ETp), potential evaporation (Ep), and potential
transpiration (Tp).

MeanStandard deviation
Coefficient ofvariation
Minimum
Maximum

R
cm

IR
Cm

Rnet
cm

ET0
cm

ETp
cm

Ep
cm

Tp
cm

89.2
22.1
24.8

1.5
0.6
37.4

87.5
21.7
24.8

111.5
3.5
3.2

57.6
1.8
3.2

17.2
0.5
3.2

40.4
1.3
3.3

43.3
126.8

0.6
2.7

42.1
124.2

107.0
121.1

55.1
62.6

16.6
18.8

38.5
43.9

(

)

′

′

(

)

∂θC2 ∂ρb S2 ∂av g2 ∂
∂C2
∂
∂g2
∂qC2
+
+
=
θDw2
+
a Dg
−
∂z
∂z v 2 ∂z
∂t
∂t
∂t
∂z
∂z
− μw,2 θC2 + μw,1 θC1 − ra,2
(

(7)

)

∂θC3 ∂
∂C3
∂qC3
′
=
θDw3
−
− μw,3 θC3 + μw,2 θC2 − ra,3
∂z
∂t
∂z
∂z

(8)

where C is the solute concentration in the liquid phase (mg L-1), S is thesolute concentration in the solid phase (mg g-1), g is the solute concen­
tration in the gas phase (mg L-1), θ is the volumetric water content (cm34

P. Nasta et al.

Agricultural Water Management 258 (2021) 107208

cm-3), ρb is the oven-dry bulk density (g cm-3), av is the air content (cm3
cm-3), q is the volumetric flux density (cm d-1), μw is the first-order rate
constant for the solute in the liquid phase (d-1), μ′ w is the first-order rate
constant providing the connection between individual nitrogen species
in the liquid phase (d-1), ra is the root nutrient uptake term (mg L-1 d-1),
Dw is the dispersion coefficient (cm2 d-1) for the liquid phase, and Dg isthe diffusion coefficient (cm2 d-1) for the gas phase. The subscripts 1, 2,
and 3 represent (NH2)2CO, NH4+, and NO3-, respectively.
The dispersion coefficient in the liquid phase, Dw, is given as:
Dw = DL|q|+ Dw τw
θ

the Miller-Miller geometric similarity hypothesis (Miller and Miller,
1956; Clausnitzer et al., 1992). In other words, the sparse point-scale soil
hydraulic functions over the heterogeneous spatial domain coalesce intoa single reference soil hydraulic function (a reference water retentionfunction, WRFref, and a reference hydraulic conductivity function,
HCFref), referring to a homogeneous spatially representative soil domain
through a set of scaling factors (Nasta et al., 2013).
A representative soil profile was set with two functional soil layers (atopsoil layer from the soil surface to the depth of 100 cm overlying a
subsoil layer between soil depths of 100 cm and 1000 cm) according tomeasurement depths described beforehand. The five unknown soil hy­
draulic parameters (α, n, θr, θs, and Ks) describing the soil water reten­tion and hydraulic conductivity functions for each of the two soil layers
have been estimated in the study area by combining the HYPRES
pedotransfer function (PTF) (Wösten et al., 1999) with the Miller-Miller
geometric similarity theory (Clausnitzer et al., 1992). Gridded values(2 km grid size) of sand, silt, and clay, oven-dry soil bulk density, andsoil organic carbon have been used as input of the HYPRES-PTF toretrieve the soil hydraulic properties of both topsoil and subsoil layers,
whereas parameter θr was set at zero. More details are presented in
Supplemental Materials (Fig. S5), and Table 4 reports the soil hydraulicparameters attributed to topsoil and subsoil in the representative soil
profile.
Initial conditions in the soil profile were set as the soil water content
values at field capacity, which is a reasonable condition on January 1st,in the wet season. Some farmers were interviewed to investigate their
irrigation management and the factors influencing the irrigation
schedule, the irrigation water volume, and the duration of each irriga­
tion event. The irrigation schedule is based on farmers’ experience indetecting soil dryness. It is impossible to quantify the exact amount ofirrigation water, and many farmers stop irrigation when infiltration
decreases or water ponds are noticed. The irrigation events start at thesoil pressure head corresponding to field capacity and end at the soilpressure head at the maize initial root water stress (Romano and Santini,2002). Therefore, to mimic farmers’ behavior as close as possible, irri­
gation was triggered with zero lag time when the pressure head at thesoil depth of 0.60 m (the depth of a tensiometer and a corresponding
control node in Hydrus-1D) reaches ψ at field capacity (ψ FC) (Bonfante
et al., 2019), given by the following physically-based analytical equation
proposed by Assouline and Or (2014):()
1 n − 1 (1− 2n)/n(10)|ψFC | =nα(9)where Dw is the molecular diffusion coefficient in free water (cm2 d-1), τw
is a unitless tortuosity factor in the liquid phase, |q| is the absolute valueof the Darcian fluid flux density (cm d-1), and DL is the longitudinal
dispersivity (cm).
2.3. Set-up of Hydrus-1D for the Sele alluvial plain
According to Nasta et al. (2020a), rainfall variability was quantified
with synthetic rainfall time series considering an alternation of two
6-month lasting seasons (see Fig. S4 in Supplementary material). A wet
season includes the months of October, November, December, January,
February, and March (when the sum of monthly average R is higher than
the sum of monthly average ETp) while the remaining months define thedry season (i.e., from April to September). In each of these two seasons,rainfall is modeled by a stochastic Poisson point process, with daily
rainfall depth μ (cm) and occurrence λ (d-1). A total of one hundred
synthetic rainfall time series of 7304 daily values (corresponding to 20
years) were generated, keeping constant parameters of the Poissonprocess and daily rainfall parent distribution in each season. In eachmodeling scenario, the synthetic time series was then used as input for
the Hydrus-1D to evaluate the effects of rainfall seasonality on the waterbalance components. The Poisson parameters calculated from the daily
rainfall values (1999–2018) measured at the Battipaglia weather station
are μ = 0.88 cm and λ = 0.39 d-1 for the wet season and μ = 0.67 cm
and λ = 0.21 d-1 for the dry season.
Table 3 reports basic statistics of ETp components available from1999 to 2018 at the Battipaglia weather station. A total of one hundred
synthetic Ep and Tp time series of 7304 daily values (corresponding to 20
years) were then generated using random numbers drawn from a normal
distribution with a mean and standard deviation in each month of theyear reported in Table 3.
A similar-media scaling method was used in this study to charac­terize the soil hydraulic behavior of the study area. This scaling
approach is conceptualized by replacing the heterogeneous porous
domain with a homogeneous soil domain when assuming the validity of

To calculate the effective soil hydraulic parameters in the represen­tative soil profile, the parametric aggregation approach presented in
Nasta and Romano (2016) was employed by obtaining ψ FC = − 207 cm.
The irrigation rate was set to 1.0 cm d-1, and the irrigation duration to 1
day.
Urea and nitrate are assumed to be present only in the dissolved
phase (the distribution coefficient, Kd =0 cm3 g-1), ammonium wasassumed to adsorb to the solid phase by setting Kd = 3.5 cm3 g-1 (Hanson
et al., 2006). The first-order decay coefficients μ’w for urea (hydrolysis)
and ammonium (nitrification) were set to 0.38 d-1 and 0.2 d-1,

Table 3
Descriptive statistics in terms of monthly mean and standard deviation of daily
values of potential evaporation, Ep, potential transpiration, Tp, calculated from
20-yr long daily temperature data at the Battipaglia weather station.
Month
Jan
Feb
Mar
Apr
May
Jun
Jul
Aug
Sep
Oct
Nov
Dec

Ep (cm d-1)

Tp (cm d-1)

MeanStandard deviation

Mean

Standard deviation

0.022
0.031
0.046
0.066
0.082
0.080
0.048
0.052
0.053
0.037
0.029
0.020

0.004
0.005
0.009
0.011
0.015
0.013
0.007
0.009
0.007
0.007
0.005
0.004

0
0
0
0
0.031
0.250
0.484
0.391
0.140
0.019
0
0

0
0
0
0
0.034
0.091
0.070
0.073
0.048
0.015
0
0

Table 4
Soil hydraulic parameters for the reference WRFref and HCFref of the topsoil andsubsoil in the soil profile.
Parameter
TopsoilθsαSubsoil

n
Ks
θs

α

n
Ks

5

Unit
cm3 cm-3
cm-1–
cm d-1
cm3 cm-3
cm-1–
cm d-10.447
0.0348
1.34
24.87
0.602
0.0231
1.29
13.16

P. Nasta et al.

Agricultural Water Management 258 (2021) 107208

respectively. Li et al. (2015) recommended that the molecular diffusion
coefficients in free water (Dw) for NH4+ and NO3- were 1.52 and1.64 cm2 d-1, respectively. We assumed Dw = 1.52 cm2 d-1 also for urea.
The longitudinal dispersivity was assumed to be 10 cm (1/10th of thesurface layer). The tortuosity factor in the liquid phase, τw, is calculatedusing the relationship proposed by Millington and Quirk (1961). Un­
limited passive uptake of both NH4+ and NO3- was allowed in the root
solute uptake model (Šimůnek and Hopmans, 2009) by specifying the
maximum allowed uptake concentration exceeding NO3- concentrations
recorded in the root zone.
Fertilization is given in the form of urea (solute 1), and a total annual
amount of 195 kg N ha-1 y-1 is applied on the soil surface (Bonfante
et al., 2019) in two different agricultural practices:

specific year at the soil surface to transform (into nitrate, solute 3) andmigrate to the soil profile bottom by entering the shallow aquifer(Fig. 3). The transit time corresponds to the difference between the peakof the nitrate breakthrough curve (BTC) at the soil profile bottom andfertilizer application date. In the case of a single treatment (agricultural
practice#1), one fertilizer application is set with 1.95 mg cm-3 of urea,(NH2)2CO (corresponding to 195 kg N ha-1), applied at the soil surface
on April 1st only in the third year. In the case of three-treatment ap­
plications (agricultural practice#2), instead, three individual sets ofmodel simulations in the third year were organized. The third year
(instead of the first year) was selected to minimize the impact of un­
known initial conditions (the first two years of simulation thus represent
a warm-up period). As mentioned above, daily R, Ep, and Tp data werestochastically generated over 20 years. This exercise was repeated one
hundred times to obtain a transit time distribution (TTD). As above, aMatlab script was used to provide model inputs and extract model re­sults automatically.

1) In a single treatment on April 1st (agricultural practice#1);
2) In a sequence of three treatments (65 kg N ha-1) on April 1st, June
1st, and August 1st (agricultural practice#2).

3. Results

Fertilization in Hydrus-1D depends on the time-dependent water flux
(precipitation) and concentration (called “cTop”) at the soil surface
boundary. We imposed Rnet (= 1.0 cm > Ep) to ensure solute entering
the soil profile. For example, the amount of 195 kg N ha-1 corresponds toa solute mass per unit area of S = 0.00195 g cm-2. Therefore, by speci­
fying the water depth, L = Rnet = 1.0 cm, the nitrate concentration inwater is cTop=S/L = 0.00195 g cm-3 (1950 mg L-1).

3.1. Results for Scenario 1: assessment of root nitrate uptake and nitrateleaching
A Monte Carlo analysis is performed with randomly generated
boundary fluxes to determine the probability distributions of annual
cumulative root nitrate uptake and nitrate leaching when fertilizer isapplied every year. Recall that new synthetic time-series of daily pre­cipitation and crop-specific potential evapotranspiration data werestochastically generated by considering the typical climate seasonality
under Mediterranean conditions.
One hundred simulated time series of daily root nitrate uptake andnitrate leaching under agricultural practice#1 and agricultural prac­
tice#2 are illustrated in Fig. 4. In general, root nitrate uptake depends
on the timing and amount of fertilization, time-variant maize charac­
teristics, and water balance, which is influenced by precipitation, cropspecific potential evapotranspiration, and irrigation. Crop-specific po­tential evapotranspiration depends on the time-variant crop coefficientthat is maximum (Kc = 1.05) for a short period from July 16th (DOY197) to August 2nd (DOY 214) during tassel, silking, and blister-kernel
phases (Fig. S2). July’s highest crop coefficient value induces maximum
crop-specific potential evapotranspiration equal to 0.532 cm d-1(Table 3). During the dry season, soil water deficit causes a reduction ofpotential transpiration into actual transpiration.
Irrigation mitigates soil water deficit during the dry season and helps
reduce root water stress by increasing actual transpiration fluxes andconsequently root nitrate uptake. Indeed, while annual mean potential
transpiration is 40.4 cm (Table 2), annual mean actual transpiration
simulated in Hydrus-1D is 32.5 cm (about a 20% reduction of Tp into

2.4. Scenario 1: assessment of root nitrate uptake and nitrate leaching
In Scenario 1, realistic agricultural management is mimicked when
the total annual amount (195 kg N ha-1) of urea, (NH2)2CO, is applied
every year for 20 years, starting from the first year. One set of one
hundred model simulations was considered when the total amount of
urea is applied with a single annual treatment (agricultural practice#1)
on April 1st, and another set of one hundred model simulations wasconsidered when the total amount of urea is applied in three treatments(agricultural practice#2) on April 1st, June 1st, and August 1st (Fig. 2).
Daily R, Ep, and Tp data were stochastically generated over 20 years toensure the impact of the Mediterranean seasonal climate regime on thewater balance. This exercise obtains probability distributions of annual
cumulative nitrate flux at the soil profile bottom and root nitrate uptake.
A Matlab script was used to automatically change daily R, Ep, and Tp
values in the input text file “ATMOSPH.IN” and extract results on nitratefrom the output text files “solute3.out” and “Obs_Node.out”.
2.5. Scenario 2: assessment of nitrate transit time
In Scenario 2, the transit time is defined as the time required for asolute (in this case urea, solute 1) applied on a single day of a single

Fig. 2. Schematic overview of Scenario 1: single and multiple treatments of fertilizer are applied every year to assess probability distributions of root nitrate uptake
and nitrate leaching.
6P. Nasta et al.

Agricultural Water Management 258 (2021) 107208

Fig. 3. Schematic overview of Scenario 2 for assessing probability distributions of nitrate transit time.

Fig. 4. Simulated root nitrate uptake (green line) and nitrate leaching (blue line) obtained using one hundred simulations considering fertilization with a) one
treatment (195 kg N ha-1) on April 1st every year, and b) three treatments (65 kg N ha-1) on April 1st, June 1st, and August 1st every year. (For interpretation of thereferences to color in this figure legend, the reader is referred to the web version of this article.)

Ta). Annual cumulative irrigation simulated in Hydrus-1D is between
36.9 cm (corresponding to 3685 m3 per hectare) and 40.7 cm (corre­
sponding to 4065 m3 per hectare) with an average amount of 38.2 cm(corresponding to 3819 m3 per hectare). While this amount seems
plausible, some farmers might irrigate using even 50% more water(http://www.agricoltura.regione.campania.it/irrigazione/prci_risultati.html). Bonfante et al. (2019) estimated total irrigation amounts on
average between 2819 m3 ha-1 and 4699 m3 ha-1. Under Mediterranean
semiarid climate conditions in northeast Spain, farmers apply on
average 81.8–86.2 cm of irrigation water for maize, considering annual
mean precipitation of 34.7–44.3 cm and grass-reference evapotranspi­
ration of 123.0–133.1 cm (Malik et al., 2019).
Nonetheless, in this study, the difference in root nitrate uptake be­
tween the two agricultural practices is caused solely by the timing andamount of fertilization since the water balance is the same under both
agricultural practices. In agricultural practice#1 (195 kg N ha-1 applied
on April 1st), the applied nitrate easily bypasses the root zone becauseroots start growing on May 9th (DOY = 129, Fig. S2). In contrast, in
agricultural practice#2 (namely the 3-steps application), nitrate derived
from fertilizer applications on June 1st and August 1st is mostly taken up

by roots. Root nitrate uptake increases linearly as a function of time (Eq.
S3), with a concurrent increase in potential evapotranspiration (Fig. S3).
Moreover, the LAI reaches its maximum value (LAI = 5.5) between July
16th (DOY 197) and August 2nd (DOY 214) (Fig. S2), maximizing po­
tential root water uptake and minimizing potential evaporation from thesoil surface (Fig. 5a). Actual transpiration reaches its maximum value atthe beginning of the mature growth phase, with favorable climate con­ditions from middle August (Fig. 5a). This interpretation highlights theimportance of adequately describing various stages of the maize
growing season, which are well-documented (Nasta and Gates, 2013;
Bonfante et al., 2019). Fig. 5b shows median daily root nitrate uptake
values under agricultural practice#1 (red line) and agricultural prac­
tice#2 (blue line) during the year, with distinct maize growth phases.
The total nitrate removed by roots during crop planting, emergence(CPE), tassel, silking, and blister-kernel phases (BK) is 0.589 mg cm-2
(58.9 kg N ha-1) and 0.648 mg cm-2 (64.8 kg N ha-1) under agricultural
practice#1 and agricultural practice#2, respectively. In contrast, themain difference in root nitrate uptake is observed during the mature
growth phase (MG), with 0.213 mg cm-2 (21.3 kg N ha-1) removed
under agricultural practice#1 and 0.476 mg cm-2 (47.6 kg N ha-1)
7

P. Nasta et al.

Agricultural Water Management 258 (2021) 107208

Therefore, the median is presented instead of the mean as a represen­tative value of each frequency distribution. The zero is the most frequentvalue (mode) in annual cumulative nitrate leaching (Fig. 6b and d)
because the nitrate flux across the soil profile bottom remains zero
during the first two years of the numerical simulations. The onetreatment fertilization (see red histograms in Fig. 6a and b) results ina median value of annual cumulative root nitrate uptake of
0.803 mg cm-2 (corresponding to 80.3 kg N ha-1) and annual cumulative
nitrate leaching of 0.749 mg cm-2 (corresponding to 74.9 kg N ha-1). In
other words, roots absorb 41% of applied nitrate, while median annual
cumulative nitrate leaching represents about 38% of applied nitrate. The
three-treatment fertilization (see yellow histograms in Fig. 6c and d)results in a median value of annual cumulative root nitrate uptake of
1.12 mg cm-2 (corresponding to 112 kg N ha-1) and annual cumulative
nitrate leaching of 0.524 mg cm-2 (corresponding to 52.4 kg N ha-1).
Almost 60% of applied nitrate is removed by root nitrate uptake, while
median annual cumulative nitrate leaching represents about 27% of
applied nitrate. Considering the total amount of fertilizer applied during
20 years is 3900 kg N ha-1, the final amount of residual nitrate in theentire 10-m-thick soil profile after 20 years is on average 796 kg N ha-1
(20.4% of applied nitrate) and 612 kg N ha-1 (15.7% of applied nitrate)
for agricultural practice#1 and #2, respectively. Similarly, cumulative
nitrate leaching after 20 years is on average 1498 kg N ha-1 and
1048 kg N ha-1 for agricultural practices #1 and #2, respectively.
In the Po River plain (Lombardy, northern Italy), under annual
rainfall conditions fairly similar to those recorded in our study area,
Perego et al. (2012) reported for maize fields at six experimental sites anannual irrigation supply ranging between 1800 m3 ha-1 and 3150 m3
ha-1. Annual fertilization amounts were between 309 kg N ha-1 and642 kg N ha-1, plant nitrogen removal between 194 kg N ha-1 (57% of
the fertilization amount) and 305 kg N ha-1 (47%), and nitrate leaching
below the root zone (z = 1.5 m) measured by suction cups ranging be­
tween 14 kg N ha-1 (4%) and 321 kg N ha-1 (50%). Demurtas et al.
(2016) reported annual mean precipitation of 56.8 cm, irrigation of
4600 m3 ha-1, mean annual fertilization of 599 kg N ha-1, root nitrate
uptake of 264 kg N ha-1 (47% of applied nitrate), nitrate leaching of
84 kg N ha-1 (15% of applied nitrate), and 250 kg N ha-1 residual nitrate(39% of applied nitrate) under Mediterranean climate conditions.

Fig. 5. Daily median values of simulated a) actual transpiration, Ta, and b) root
nitrate uptake considering single (red line) and triple (blue line) application
treatments. The black asterisks indicate fertilizer applications on April 1st, June
1st, and August 1st. Vertical gray dotted lines delimit a dormant period (DP),maize planting and emergence (CPE), tassel, silking, and blister-kernel phases(BK), and mature growth phase (MG). (For interpretation of the references tocolor in this figure legend, the reader is referred to the web version ofthis article.)

removed under agricultural practice#2. In the latter case, the third
treatment applied on August 1st (third black asterisk on DOY 213 inFig. 5) boosts root nitrate uptake (blue line in Fig. 5b) during the mature
growth phase. The impact of harvest timing (usually from middle
September to middle October) on root nitrate uptake is negligible underboth agricultural practices.
Fig. 6 shows the frequency distribution of annual cumulative root
solute uptake and nitrate leaching based on one hundred forward sim­
ulations. Two normality tests (Lilliefors and Kolmogorov-Smirnov
executed in Matlab) reject the null hypothesis that both annual cumu­
lative root solute uptake and nitrate leaching (under both agricultural
practices) come from a normal distribution at the 5% significance level.

Fig. 6. Frequency distribution of 20 annual cumulative nitrate fluxes in one hundred forward modeling simulations: a) root nitrate uptake in agricultural practice#1;b) nitrate leaching in agricultural practice#1; c) root nitrate uptake in agricultural practice#2; d) nitrate leaching in agricultural practice#2.

8

P. Nasta et al.

Agricultural Water Management 258 (2021) 107208

3.2. Results for Scenario 2: assessment of nitrate transit time
A Monte Carlo analysis was performed by generating upper bound­ary fluxes to determine a probabilistic distribution of nitrate transit timevalues (TTD). Scenario #2 involves an ensemble of one hundred nu­
merical simulations with Hydrus-1D subjected to daily rainfall data
generated by the Poisson distribution parameters pertaining to the dry
and wet seasons, and Ep and Tp data with normally distributed param­
eters (mean and standard deviation) in each month. Once again, recall
that the fertilizer is applied only in the third year, either with a single
treatment (agricultural practice #1) or with a three-step treatment(agricultural practice #2).
Fig. 7 refers to agricultural practice #1 (i.e., a single treatment ofurea) and depicts the ensemble of one hundred nitrate concentrations at
soil depths of 300 cm, 600 cm, and 1000 cm (Fig. 7a) and corresponding
normal distributions of transit time values (Fig. 7b). These graphs showthat nitrate concentration gradually propagates downward through thevadose zone, with BTCs characterized by decreasing magnitude andincreasing duration with depth (Fig. 7a). It is hypothesized that below
the root zone (z = 300 cm, blue lines), the nitrate BTC peaks depend on
the combination of daily Rnet and ETp following the nitrate spill. The
BTCs are affected by water infiltration induced by rainfall, mainly
occurring in the wet season. To test this hypothesis, the Pearson’s cor­
relation coefficients between the cumulated rainfall and transit timevalues at the soil depths of 300 cm, 600 cm, and 1000 cm were − 0.80,− 0.85, − 0.85, respectively. This indicates that the TTD values depend
highly on Rnet, accelerating or slowing down nitrate transport throughthe unsaturated zone. The standard deviation of TTDs increases with soil
depth (65.4 days, 145.3 days, and 181.4 days at the soil depths of
300 cm, 600 cm, and 1000 cm, respectively).
When considering agricultural practice #2 (i.e., three applications ofurea), Fig. 8 shows the normal probability distributions of the ensembleof one hundred nitrate concentrations at soil depths of 300 cm, 600 cm,
and 1000 cm. In general, transit times decrease when treatment is
individually applied on April 1st, June 1st, and August 1st. The transit
times relevant to the treatment applied in April are very similar to thosepresented previously and associated with a single treatment of
195 kg N ha-1. For the treatments applied on June 1st and August 1st,
the mean transit times are 2707 days (with a standard deviation of 186days) and 2650 days (with a standard deviation of 185 days),respectively.
It is important to point out that, in the present study, TTDs areaffected only by the upper boundary forcing, while other factors were
kept fixed, namely i) the soil hydraulic parameters in the two functionalsoil layers (influencing flow velocity across the unsaturated zone), andii) the reaction and transport parameters, especially the nitrate

Fig. 8. Probability density function (pdf) of transit times obtained by runningone hundred simulations at z = 300 cm (blue line), z = 600 cm (green line),and z = 1000 cm (cyan line) when the treatment is individually applied on a)
April 1st, b) June 1st, c) August 1st. (For interpretation of the references tocolor in this figure legend, the reader is referred to the web version ofthis article.)longitudinal dispersivity (affecting the solute transport process).
4. Discussion
In the present study, the nitrate fate was evaluated in a representa­
tive 10-m-thick soil profile by using a process-oriented hydrologicalmodel supported by a relatively large number of soil data (Nasta et al.,2020b) as well as parameter values describing the nitrogen transport
and transformation taken from the literature (Veizaga et al., 2016).
However, because of a lack of direct local measurements on nitrate data,
our one-dimensional numerical results should not be strictly considered
for actual agronomic applications. Indeed, our model simulations pro­vide useful indications of nitrate transport in the area and, more
importantly, can be integrated with map overlay procedures for delin­eating nitrate vulnerable zones (Matzeu et al., 2017; Catani et al., 2020;
Vigliotti et al., 2020; Tufano et al., 2020). The method presented in thisstudy identifies process-oriented indicators accounting for the local
climatic and environmental conditions. This is, therefore, more effective
and understandable to landowners and stakeholders. To our knowledge,
this study provides a novel protocol to assess TTDs characterizing nitratetransport within the soil profile to support the mapping of NVZs (Fusco
Fig. 7. a) Simulated concentrations at
z = 300 cm (blue line), z = 600 cm (green line),and z = 1000 cm (cyan line) by using one
hundred simulations when considering one
treatment (195 kg N ha-1) on April 1st in thethird year and b) probability density function(pdf) of transit times obtained by running one
hundred simulations at z = 300 cm (blue line),z = 600 cm (green line), and z = 1000 cm(cyan line). Mean values of transit time distri­butions (TTD) are reported. (For interpretationof the references to color in this figure legend,
the reader is referred to the web version of thisarticle.)9P. Nasta et al.

Agricultural Water Management 258 (2021) 107208

et al., 2020). If model simulations are assumed valid, the nitrate transit
time can be viewed as a “dynamic” indicator of groundwater vulnera­
bility and an efficient way to evaluate the protective function exerted by
the vadose zone (Matteau et al., 2019; Turkeltaub et al., 2020).
Quantification of nitrate fate and transport in the vadose zone de­pends on several factors that may significantly affect the accuracy of theestimated nitrate transit time (Mattern and Vanclooster, 2010). Water
flow depends on the depth to the groundwater table, soil hydraulicproperties of the vadose zone, and maize characteristics influencing
evapotranspiration fluxes. Nitrogen species fate and transport depend on
their transport and reaction parameters, initial concentrations
throughout the soil profile, and the amount and frequency of fertilizer
applications.
It is very difficult to characterize soil layering and its soil hydraulicproperties throughout large-scale domains. In this study, only informa­tion on soil physical properties within the topmost meter of the vadosezone was used to artificially build two functional soil layers over a 10-mthick soil profile. HYPRES-PTF was employed in each soil layer toretrieve five unknown soil hydraulic parameters (α, n, θr, θs, and Ks)
describing the soil water retention and hydraulic conductivity functions.
Other existing PTFs developed on European datasets could be tested aswell (Weynants et al., 2009; Szabó et al., 2021; Nasta et al., 2021). The
eleven lithological core drilling soil profiles evidence relatively uniform
soil physical properties in the vadose zone. Nevertheless, the impact ofsoil vertical variability on the decay of saturated hydraulic conductivity
induced by soil compaction in the deepest soil layers is expected (Tur­keltaub et al., 2021). The scenario-based simulations are not supported
by monitoring key state variables (water contents, pressure head data)or water fluxes throughout the soil profile. The cost of the instrumen­
tation might be high and unfeasible. Farmers do not recommend theinstallation of sensors beneath the maize fields because of tillage oper­
ations. Some studies interestingly calibrate physical models using soft
data approaches such as isotope-based techniques (Stumpp et al., 2012;
Sacchi et al., 2013; Sprenger et al., 2016a; Boumaiza et al., 2021).
Future improvements will integrate key measurements of nitrateabsorbed by roots, nitrate leaching into the shallow aquifer, and nitrateconcentrations measured throughout the soil profile based on stable
isotope data. In a subsequent study, farmers’ collaboration will be
requested to record the exact amount and frequency of fertilization,
fertilizer type, and irrigation depths during the growing season. Our
results produce annual root nitrate uptake comparable to previous
studies in cornfields (Tafteh and Sepaskhah, 2012).
In a subsequent step, the available gridded soil hydraulic properties
will be used in the Hydrus-1D package by assuming the spatial domain
represented by independent columns. Hydrus-1D will be coupled to
MODFLOW (Beegum et al., 2018, 2019) to describe the nitrate transport
influenced by topography and depth to the water table in the 3D domain(Prudic et al., 2004; Zeng et al., 2019).
Another future study will investigate the impact of point source
pollution caused by livestock farming (Piccini et al., 2016). The
contribution of localized applications of high amounts of buffalo manureon nitrate mineralization will be assessed. In this regard, a refinement ofavailable models is needed to consider the process of nitrate immobili­
zation into the microbial biomass that commonly occurs for organic
amendments with a C/N ratio above ~30 (Hodge et al., 2000; Bonanomi
et al., 2014).

can be enhanced in the future by investing in monitoring network in­
frastructures. Effective soil hydraulic and nitrate transport parameters
can be estimated using ad hoc inverse modeling procedures supported
by available vadose zone data. The 20 years long model simulations
show that agricultural practices strongly influence root nitrate uptake
and nitrate leaching. In all cases, nitrate transit times are about two
years to bypass the root zone and about 7–8 years to reach the shallow
aquifer. Groundwater nitrate contamination is a consequence of nitrateleaching and will be subject to future investigations. The proposed
approach can produce vulnerability maps for nitrate pollution or otherpollutant types for other areas of the Campanian alluvial plain and othersimilar regions worldwide.
Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Acknowledgments
The study reported in this paper was supported in part by Istituto
Zooprofilattico Sperimentale del Mezzogiorno (IZSM) through the“Transparency in Campania” Action Plan, and MiUR-PRIN Project“WATer mixing in the critical ZONe: observations and predictions underenvironmental changes – WATZON” (grant 2017SL7ABC). We are
grateful to the Managing Director of IZSM, Antonio Limone.
Appendix A. Supporting information
Supplementary data associated with this article can be found in the
online version at doi:10.1016/j.agwat.2021.107208.
References
Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Crop evapotranspiration: guidelines
for computing crop water requirements. Irrigation and Drainage Paper No. 56, Foodand Agriculture Organization of the United Nations (FAO), Rome, Italy, 300p.
Antonakos, A.K., Lambrakis, N.J., 2007. Development and testing of three hybrid
methods for the assessment of aquifer vulnerability to nitrates, based on the drastic
model, an example from NE Korinthia, Greece. J. Hydrol. 333, 288–304.
Assouline, S., Or, D., 2014. The concept of field capacity revisited: defining intrinsic
static and dynamic criteria for soil internal drainage dynamics. Water Resour. Res.
50 https://doi.org/10.1002/2014WR015475.
Beegum, S., Šimůnek, J., Szymkiewicz, A., Sudheer, K.P., Nambi, I.M., 2018. Updating
the Coupling Algorithm between HYDRUS and MODFLOW in the HYDRUS Package
for MODFLOW. Vadose Zone J. 17 (1), 8. https://doi.org/10.2136/vzj2018.02.0034(180034).
Beegum, S., Šimůnek, J., Szymkiewicz, A., Sudheer, K.P., Nambi, I.M., 2019.
Implementation of solute transport in the vadose zone into the ‘HYDRUS package for
MODFLOW’. Groundwater 57 (3), 392–408. https://doi.org/10.1111/gwat.12815.
Bonanomi, G., D’Ascoli, R., Scotti, R., Gaglione, S.A., Caceres, M.G., Sultana, S., Rao, M.,
Zoina, A., 2014. Soil quality recovery and crop yield enhancement by combined
application of compost and wood to vegetables grown under plastic tunnels. Agric.
Ecosyst. Environ. 192, 1–7.
Bonfante, A., Monaco, E., Manna, P., De Mascellis, R., Basile, A., Buonanno, M.,
Cantilena, G., Esposito, A., Tedeschi, A., De Michele, C., Belfiore, O., Catapano, I.,
Ludeno, G., Salinas, K., Brook, A., 2019. LCIS DSS—an irrigation supporting system
for water use efficiency improvement in precision agriculture: a maize case study.
Agric. Syst. 176, 102646.
Boumaiza, L., Chesnaux, R., Walter, J., Stumpp, C., 2021. Constraining a flow model withfield measurements to assess water transit time through a vadose zone. Groundwater
59, 417–427. https://doi.org/10.1111/gwat.13056.
Braden, H., 1985. Ein Energiehaushalts- und Verdunstungsmodell for Wasser und
Stoffhaushaltsuntersuchungen landwirtschaftlich genutzer Einzugsgebiete. Mittel
Dtsch. Bodenkd. Geselschaft 42, 294–299.
Bradshaw, J.K., Radcliffe, D.E., Šimůnek, J., Wunsch, A., McCray, J.E., 2013. Nitrogenfate and transport in a conventional onsite wastewater treatment system installed ina clay soil: a nitrogen chain model. Vadose Zone J. 12 (3), 20. https://doi.org/
10.2136/vzj2012.0150.
Busico, G., Mastrocicco, M., Cuoco, E., Sirna, M., Tedesco, D., 2019. Protection fromnatural and anthropogenic sources: a new rating methodology to delineate “Nitrate
Vulnerable Zones”. Environ. Earth Sci. 78, 104. https://doi.org/10.1007/s12665019-8118-2.

5. Concluding remarks
Numerical simulations of hydrological processes and nitrate trans­
port provide useful process-oriented indicators that can be integrated
with the information retrieved from GIS-based approaches to map NVZs
and quantify potential hazards of detrimental anthropogenic environ­mental impacts. The use of process-oriented models is often prevented
by the lack of key measurements and farmers’ information (e.g., about
irrigation and fertilization). Nevertheless, the use of numerical models
10P. Nasta et al.

Agricultural Water Management 258 (2021) 107208

Busico, G., Kazakis, N., Colombani, N., Khosravi, K., Voudouris, K., Mastrocicco, M.,
2020. The importance of incorporating denitrification in the assessment ofgroundwater vulnerability. Appl. Sci. 10 (2328), 1–11.
Cassa per il Mezzogiorno , 1983. Idrogeologia dell’Italia Centro-Meridionale. Quaderno
nr. 4/2, Roma (in Italian).
Catani, V., Zuzolo, D., Esposito, L., Albanese, S., Pagnozzi, M., Fiorillo, F., De Vico, B.,
Cicchella, D., 2020. A new approach for aquifer vulnerability assessment: the casestudy of Campania Plain. Water Resour. Manag. 34, 819–834. https://doi.org/
10.1007/s11269-019-02476-5.
Clausnitzer, V., Hopmans, J.W., Nielsen, D.R., 1992. Simultaneous scaling of soil water
retention and hydraulic conductivity curves. Water Resour. Res. 28, 19–31.
Corbari, C., Ben Charfi, I., Mancini, M., 2021. Optimizing irrigation water use efficiency
for tomato and maize fields across italy combining remote sensing data and the
AquaCrop Model. Hydrology 8 (8), 39. https://doi.org/10.3390/hydrology8010039.
De Vita, P., Allocca, V., Celico, F., Fabbrocino, S., Cesaria, M., Monacelli, G., Musilli, I.,
Piscopo, V., Scalise, A.R., Summa, G., Tranfaglia, G., Celico, P., 2018. Hydrogeology
of continental southern Italy. J. Maps 14 (2), 230–241.
Deb, S.K., Sharma, P., Shukla, M.K., Ashigh, J., Šimůnek, J., 2016. Numerical evaluationof nitrate distributions in the onion root zone under conventional furrow fertigation.
J. Hydrol. Eng. 21 (2), 05015026 https://doi.org/10.1061/(ASCE)HE.19435584.0001304.
Demurtas, C.E., Seddaiu, G., Ledda, L., Cappaia, C., Doro, L., Carletti, A., Roggero, P.P.,
2016. Replacing organic with mineral N fertilization does not reduce nitrate leachingin double crop forage systems under Mediterranean conditions. Agric. Ecosyst.
Environ. 219, 83–92.
Doerfliger, N., Jeannin, P.Y., Zwahlen, F., 1999. Water vulnerability assessment in karst
environments: a new method of defining protection areas using a multiattribute
approach and GIS tools (EPIK method). Environ. Geol. 39, 165–176.
Ducci, D., Della Morte, R., Mottola, A., Onorati, G., Pugliano, G., 2019. Nitrate trends ingroundwater of the Campania region (southern Italy). Environ. Sci. Pollut. Res. 26,
2120–2131.
Feddes, R.A., Kowalik, P.J., Zarandny, H., 1978. Simulation of field water use and crop
yield. John Wiley & Sons, New York, NY.
Fusco, F., Allocca, V., Coda, S., Cusano, D., Tufano, R., De Vita, P., 2020. Quantitative
assessment of specific vulnerability to nitrate pollution of shallow alluvial aquifers
by process-based and empirical approaches. Water 12 (1), 269. https://doi.org/
10.3390/w12010269.
Gómez-Hernández, J.J., Wen, X.-H., 1994. Probabilistic assessment of travel times ingroundwater modeling. Stoch. Hydrol. Hydraul. 8, 19–55.
Hanson, B.R., Šimůnek, J., Hopmans, J.W., 2006. Evaluation of urea–ammonium–nitratefertigation with drip irrigation using numerical modeling. Agric. Water Manag. 86,
102–113.
Hargreaves, G.H., Allen, R.G., 2003. History and evaluation of Hargreaves
evapotranspiration equation. J. Irrig. Drain. Eng. 129 (1), 53–63. https://doi.org/
10.1061/(ASCE)0733-9437(2003)129:1(53).
Hodge, A., Robinson, D., Fitter, A., 2000. Are microorganisms more effective than plants
at competing for nitrogen? Trends Plant Sci. 5, 304–308.
Huan, H., Wang, J., Teng, Y., 2012. Assessment and validation of groundwater
vulnerability to nitrate based on a modified DRASTIC model: a case study in Jilin
City of northeast China. Sci. Total Environ. 440, 14–23.
Infascelli, R., Pelorosso, R., Boccia, L., 2009. Spatial assessment of animal manure
spreading and groundwater nitrate pollution. Geospat. Health 4 (1), 27–38.
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.
Malik, W., Isla, R., Dechmi, F., 2019. DSSAT-CERES-maize modelling to improveirrigation and nitrogen management practices under Mediterranean conditions.
Agric. Water Manag. 213, 298–308.
Matteau, J.-P., Gumiere, S.J., Gallichand, J., Létourneau, G., Khiari, L., Gasser, M.-O.,
Michaud, A., 2019. Coupling of a nitrate production model with HYDRUS to predict
nitrate leaching. Agric. Water Manag. 213, 616–626.
Mattern, S., Vanclooster, M., 2010. Estimating travel time of recharge water through a
deep vadose zone using a transfer function model. Environ. Fluid Mech. 10,
121–135.
Matzeu, A., Secci, R., Uras, G., 2017. Methodological approach to assessment ofgroundwater contamination risk in an agricultural area. Agric. Water Manag. 184,
46–58.
McLay, C.D.A., Dragten, R., Sparling, G., Selvarajah, N., 2001. Predicting groundwater
nitrate concentrations in a region of mixed agricultural land use: a comparison of
three approaches. Environ. Pollut. 115, 191–204.
Miller, E.E., Miller, R.D., 1956. Physical theory for capillary flow phenomena. J. Appl.
Phys. 27, 324–332.
Millington, R.J., Quirk, J.M., 1961. Permeability of porous solids. Trans. Faraday Soc. 57,
1200–1207.
Mokari, E., Shukla, M.K., Šimůnek, J., Fernandez, J.L., 2019. Numerical modeling of
nitrate in a flood-irrigated pecan orchard. Soil Sci. Soc. Am. J. 83, 555–564.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resour. Res. 12, 513–522.
Mylevaganam, S., Ray, C., 2016. The assessment of groundwater vulnerability due toleaching of chemicals: the review of attenuation factor. Open J. Soil Sci. 6, 9–20.
Nasta, P., Gates, J.B., 2013. Plot-scale modeling of soil water dynamics and impacts ofdrought conditions beneath rainfed maize in Eastern Nebraska. Agric. Water Manag.
128, 120–130.
Nasta, P., Romano, N., 2016. Use of a flux-based field capacity criterion to identify
effective hydraulic parameters of layered soil profiles subjected to synthetic drainage
experiments. Water Resour. Res. 52 https://doi.org/10.1002/2015WR016979.

Nasta, P., Romano, N., Assouline, S., Vrugt, J.A., Hopmans, J.W., 2013. Prediction of
spatially-variable unsaturated hydraulic conductivity using scaled particle-size
distribution functions. Water Resour. Res. 49 https://doi.org/10.1002/wrcr.20255.
Nasta, P., Allocca, C., Deidda, R., Romano, N., 2020a. Assessing the impact of seasonalrainfall anomalies on catchment-scale water balance components. Hydrol. Earth
Syst. Sci. 24, 1–17.
Nasta, P., Palladino, M., Sica, B., Pizzolante, A., Trifuoggi, M., Toscanesi, M., Giarra, A.,
D’Auria, J., Nicodemo, F., Mazzitelli, C., Lazzaro, U., Di Fiore, P., Romano, N.,2020b. Evaluating pedotransfer functions for predicting soil bulk density using
hierarchical mapping information in Campania. Geoderma Reg. 21, e00267 https://doi.org/10.1016/j.geodrs.2020.e00267.
Nasta, P., Szabó, B., Romano, N., 2021. Evaluation of Pedotransfer Functions forpredicting soil hydraulic properties: A voyage from regional to field scales acrossEurope. J. Hydrol. Reg. Stud. 37 https://doi.org/10.1016/j.ejrh.2021.100903.
National Research Council, 1993. Ground Water Vulnerability Assessment:
Contamination Potential under Conditions of Uncertainty. National Academy Press,
Washington, DC, USA.
Panagopoulos, G.P., Antonakos, A.K., Lambrakis, N.J., 2006. Optimization of the
DRASTIC method for groundwater vulnerability assessment via the use of simple
statistical methods and GIS. Hydrogeol. J. 14, 894–911.
Pelosi, A., Medina, H., Villani, P., D’Urso, G., Chirico, G.B., 2016. Probabilistic
forecasting of reference evapotranspiration with a limited area ensemble prediction
system. Agric. Water Manag. 178, 106–118.
Perego, A., Basile, A., Bonfante, A., De Mascellis, R., Terribile, F., Brenna, S., Acutis, M.,
2012. Nitrate leaching under maize cropping systems in Po Valley (Italy). Agric.
Ecosyst. Environ. 147, 57–65. https://doi.org/10.1016/j.agee.2011.06.014.
Piccini, C., Di Bene, C., Farina, R., Pennelli, B., Napoli, R., 2016. Assessing nitrogen use
efficiency and nitrogen loss in a forage-based system using a modeling approach.
Agronomy 6, 23. https://doi.org/10.3390/agronomy6020023.
Prudic, D.E., Konikow, L.F., Bant, E.R. , 2004. A new stream- flow routing (SFR1)
package to simulate stream-aquifer interaction with MODFLOW-2000, U.S.
Geological Survey Open-File Report 2004–1042, 95p.
Ramos, T.B., Šimůnek, J., Gonçalves, 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.
Ritchie, J.T., 1972. Model for predicting evaporation from a row crop with incomplete
cover. Water Resour. Res. 8, 1204–1213.
Romano, N., Santini, A., 2002. Water retention and storage: field. In: Dane, J.H., Topp, G.
C. (Eds.), Methods of Soil Analysis, Part 4: Physical Methods, SSSA Book Series N.5.
Wiley, Madison, WI, USA, pp. 721–738. https://doi.org/10.2136/sssabookser5.4.
c26 (ISBN 0-89118-841-X).
Rupert, M.G., 2001. Calibration of the DRASTIC Ground Water Vulnerability Mapping
Method. Ground Water 39, 625–630.
Sacchi, E., Acutis, M., Bartoli, M., Brenna, S., Delconte, C.A., Laini, A., Pennisi, M., 2013.
Origin and fate of nitrates in groundwater from the central Po plain: insights fromisotopic investigations. Appl. Geochem. 34, 164–180.
Šimůnek, J., Hopmans, J.W., 2009. Modeling compensated root water and nutrient
uptake. Ecol. Model. 220, 505–521.
Šimůnek, J., Šejna, M., van Genuchten, M.T., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J. 15, 1–25.
https://doi.org/10.2136/vzj2016.04.0033.
Sprenger, M., Erhardt, M., Riedel, M., Weiler, M., 2016a. Historical tracking of nitrate incontrasting vineyards using water isotopes and nitrate depth profiles. Agric. Ecosyst.
Environ. 222, 185–192.
Sprenger, M., Seeger, S., Blume, T., Weiler, M., 2016b. Travel times in the vadose zone:
variability in space and time. Water Resour. Res. 52, 5727–5754. https://doi.org/
10.1002/2015WR018077.
Stumpp, C., Stichler, W., Kandolf, M., Šimůnek, J., 2012. Effects of land cover andfertilization method on water flow and solute transport in five lysimeters: a longterm study using stable water isotopes. Vadose Zone J. 11 (1) https://doi.org/
10.2136/vzj2011.0075.
Szabó, B., Weynants, M., Weber, T.K., 2021. Updated European hydraulic pedotransfer
functions with communicated uncertainties in the predicted variables (euptfv2).
Geosci. Model Dev. 14, 151–175.
Tafteh, A., Sepaskhah, A.R., 2012. Application of HYDRUS-1D model for simulatingwater and nitrate leaching from continuous and alternate furrow irrigated rapeseed
and maize fields. Agric. Water Manag. 113, 19–29. https://doi.org/10.1016/j.
agwat.2012.06.011.
Tufano, R., Allocca, V., Coda, S., Cusano, D., Fusco, F., Nicodemo, F., Pizzolante, A., De
Vita, P., 2020. Groundwater vulnerability of principal aquifers of the Campania
region (southern Italy). J. Maps 16, 565–576. https://doi.org/10.1080/
17445647.2020.1787887.
Turkeltaub, T., Ascott, M.J., Gooddy, D.C., Jia, X., Shao, M.-A., Binley, A., 2020.
Prediction of regional-scale groundwater recharge and nitrate storage in the vadosezone: a comparison between a global model and a regional model. Hydrol. Proc. 34,
3347–3357. https://doi.org/10.1002/hyp.13834.
Turkeltaub, T., Jia, X., Zhu, Y., Shao, M.-A., Binley, A., 2021. A comparative study ofconceptual model complexity to describe water flow and nitrate transport in deep
unsaturated Loess. Water Resour. Res. 57 https://doi.org/10.1029/2020WR029250.
van Genuchten, M.T., 1980. A closed form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
Van Stempvoort, D., Evert, L., Wassenaar, L., 1993. Aquifer vulnerability index: a GIS
compatible method for groundwater vulnerability mapping. Can. Water Resour. J.
18, 25–37.

11

P. Nasta et al.

Agricultural Water Management 258 (2021) 107208

Veizaga, E.A., Rodriguez, L., Ocampo, C.J., 2016. Investigating nitrate dynamics in afine-textured soil affected by feedlot effluents. J. Contam. Hydrol. 193, 21–34.
Vigliotti, M., Busico, G., Ruberti, D., 2020. Assessment of the vulnerability to agricultural
nitrate in two highly diversified environmental settings. Environments 7, 80.
https://doi.org/10.3390/environments7100080.
Wall, D.P., Delgado, A., O’Sullivan, L., Creamer, R.E., Trajanov, A., Kuzmanovski, V.,
Bugge Henriksen, C., Debeljak, M., 2020. A decision support model for assessing thewater regulation and purification potential of agricultural soils across Europe. Front.
Sustain. Food Syst. 4, 115. https://doi.org/10.3389/fsufs.2020.00115.

Weynants, M., Vereecken, H., Javaux, M., 2009. Revisiting Vereecken pedotransfer
functions: introducing a closed-form hydraulic model. Vadose Zone J. 8 (1), 86–95.
https://doi.org/10.2136/vzj2008.0062.
Wösten, J.H.M., Lilly, A., Nemes, A., Le Bas, C., 1999. Development and use of a databaseof hydraulic properties of European soils. Geoderma 90, 169–185.
Zeng, J., Yang, J., Zha, Y., Shi, L., 2019. Capturing soil-water and groundwaterinteractions with an iterative feedback coupling scheme: new HYDRUS package for
MODFLOW. Hydrol. Earth Syst. Sci. 23, 637–655. https://doi.org/10.5194/hess-23637-2019.

12