Journal of Hydrology 367 (2009) 138–149

Contents lists available at ScienceDirect

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

A root zone modelling approach to estimating groundwater recharge from
irrigated areas
J. Jiménez-Martínez a,b,*, T.H. Skaggs b, M.Th. van Genuchten c, L. Candela aaDepartment of Geotechnical Engineering and Geosciences, Technical University of Catalonia, UPC, Gran Capitan s.n., 08034 Barcelona, Spain
US Salinity Laboratory, USDA-ARS, Riverside, CA, USAc
Department of Mechanical Engineering, LTTC/COPPE, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro, Brazilb

a r t i c l e

i n f oArticle history:
Received 3 June 2008
Received in revised form 22 December 2008
Accepted 9 January 2009

Associate Editor: Magnus Persson
Keywords:
Root zone
Aquifer recharge
Irrigation
Campo de Cartagena
Water content
Pressure heads u m m a r y
In irrigated semi-arid and arid regions, accurate knowledge of groundwater recharge is important for thesustainable management of scarce water resources. The Campo de Cartagena area of southeast Spain is asemi-arid region where irrigation return ﬂow accounts for a substantial portion of recharge. In this study
we estimated irrigation return ﬂow using a root zone modelling approach in which irrigation, evapotranspiration, and soil moisture dynamics for speciﬁc crops and irrigation regimes were simulated with
the HYDRUS-1D software package. The model was calibrated using ﬁeld data collected in an experimental
plot. Good agreement was achieved between the HYDRUS-1D simulations and ﬁeld measurements madeunder melon and lettuce crops. The simulations indicated that water use by the crops was below potential levels despite regular irrigation. The fraction of applied water (irrigation plus precipitation) going torecharge ranged from 22% for a summer melon crop to 68% for a fall lettuce crop. In total, we estimatethat irrigation of annual fruits and vegetables produces 26 hm3 y1 of groundwater recharge to the top
unconﬁned aquifer. This estimate does not include important irrigated perennial crops in the region, suchas artichoke and citrus. Overall, the results suggest a greater amount of irrigation return ﬂow in theCampo de Cartagena region than was previously estimated.
Ó 2009 Elsevier B.V. All rights reserved.

Introduction
Estimating aquifer recharge is important for determining water
resource availability and assessing aquifer vulnerability to pollutants (Scanlon et al., 2002). Recharge estimation can be difﬁcult,particularly in arid and semi-arid regions where water tables are
typically deep and recharge is predominately focused recharge thatemanates from topological depressions such as streams and lakes.
The recharge rate is limited by the availability of water at the land
surface, which is controlled by temporally and spatially variableclimatic factors such as precipitation and evapotranspiration(Scanlon et al., 2002). In some basins the estimation of recharge
is additionally complicated by irrigation, which may simultaneously remove water from focused recharge sources while creating new sources of diffuse recharge. In irrigated regions, accurate
knowledge of recharge, evaporation, and transpiration is especially
important for the sustainable management of scarce water resources (e.g., Gartuza-Payán et al., 1998).

* Corresponding author. Address: Department of Geotechnical Engineering andGeosciences, Technical University of Catalonia, UPC, Gran Capitan s.n., 08034
Barcelona, Spain. Tel.: +34 93 401733; fax: +34 93 4017251.
E-mail address: joaquin.jimenez@upc.edu (J. Jiménez-Martínez).
0022-1694/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved.
doi:10.1016/j.jhydrol.2009.01.002

Several methods have been used to estimate groundwater recharge with varying degrees of success (reviews include Gee andHillel, 1988; Allison et al., 1994; Scanlon et al., 2002; de Vries
and Simmers, 2002). The methods can be loosely grouped into
three categories depending on whether the focus of the method
is surface water, the vadose zone, or the saturated zone. In eachof these cases, physical and tracer techniques are possible, as are
numerical modelling approaches. The best choice for a particularsituation depends upon the spatial and temporal scales being considered and the intended application of the recharge estimate(Scanlon et al., 2002).
In this work, we studied recharge in the Campo de Cartagena
area of southeast Spain, a semi-arid region where irrigated agriculture is prevalent. Our objective was to test a root zone modelling
approach that can be used to estimate aquifer recharge emanatingfrom irrigated farmland. The modelling approach utilized HYDRUS1D (Šimůnek et al., 2005), a well-known computer model thatsimulates water, heat, and solute movement in variably saturated
porous media. A critical element of water balance and modelling
approaches to recharge estimation is determining actual evapotranspiration rates, which can be below potential rates for long
periods of time in arid and semi-arid regions, even in irrigated systems. Among other general recharge modelling efforts (e.g., Ragab
et al., 1997; Finch, 1998; Zhang et al., 1999; Brunner et al., 2004),

J. Jiménez-Martínez et al. / Journal of Hydrology 367 (2009) 138–149

Kendy et al. (2003) evaluated recharge speciﬁcally for irrigatedcropland using a model in which soil water ﬂow was governedby a tipping-bucket-type mechanism and actual transpiration
was computed based on the soil water status using a method introduced by Campbell and Norman (1998). In our work, root zone
moisture dynamics are simulated with the Richards equation androot water uptake and transpiration are calculated according toFeddes et al. (1978).

Campo de Cartagena
The Campo de Cartagena plain comprises an area of
1440 km2 in the Region de Murcia in southeast Spain (Fig. 1).
The climate is Mediterranean with an average annual rainfallof 300 mm and a mean annual temperature of 18 °C. Estimatesof annual potential evapotranspiration (ETP) range from 800 to
1200 mm y1, depending on the estimation method (Sánchez
et al., 1989). Tourism is an important industry due to the region’s mild climate and many beaches, parks, and golf resorts.
However, the dominant industry and land use is agriculture,both irrigated and rainfed. Irrigated farmland comprises an areaof approximately 299 km2, with 128.1 km2 of annual row crops(principally lettuce and melon), 34.1 km2 of perennial vegetables(principally artichoke), and 136.8 km2 of fruit trees (principally
citrus). These crops are an important source of fruits and vegetables for the European Union. Drip irrigation is used widely intheir production due to the scarcity of water resources and theneed for water conservation. Recently, however, drought conditions have worsened, a deterioration that many attribute to climatechange. The future of irrigated agriculture in Campo de Cartagena isin doubt; the United Nations Food and Agriculture Organization has
identiﬁed current water shortage and desertiﬁcation problems insoutheast Spain as possibly being harbingers of what may become
a global food crisis (New York Times, June 3, 2008).
Geologically, the Campo de Cartagena plain consists of Neogene
and Quaternary materials of sedimentary origin (gravel, sand, silt,
and clay) overlying the Betic Complex which is comprised of External Zones (sedimentary materials) and Internal Zones (metamorphic materials). The top unconﬁned aquifer, of Quaternary age,extends over 1135 km2 with an average thickness of 50 m. Piezometric levels measured in the inner part of the study area arearound 20–30 m below the surface, while near the coast they areusually around 1–2 m. Transmissivity values vary widely depending on location and geological material.
A preliminary study by the Spanish Geological Survey (IGME,
1994) estimated that total recharge to the top unconﬁned aquifer
was about 69 hm3 y1, 46 hm3 y1 due to natural recharge and23 hm3 y1 due to irrigation return ﬂow (where 1 hm3 = 106 m3= 810.71 ac ft). The SGS study used the Thornthwaite method(Thornthwaite, 1948) to estimate natural recharge and, depending
on data availability, a combination of methods to estimate irrigation return ﬂows. For regions where crop and irrigation data wereavailable, irrigation return ﬂows were calculated as the difference
between the applied water and the potential crop water use. For
other areas where only irrigation data were available, irrigation
efﬁciency coefﬁcients for different irrigation methods (e.g., drip,
furrow, ﬂooding) were used to determine the fraction of irrigationwater going to recharge. Thus the SGS estimate for return ﬂow was
based on irrigation water applications only and did not consider
water additions due to precipitation. Instead, precipitation-based
recharge from irrigated farmland was implicitly included in the
SGS estimate of natural recharge, which was a single value calculated for the entire region.
In our work, we calculate irrigation return ﬂows based on combined water inputs from irrigation and precipitation. So that we139

may better compare our results with those of the SGS, we makethe following assumptions and extrapolations about the SGS estimates. According to the land use ﬁgures noted above, irrigatedfarmland covers 26% of the top unconﬁned aquifer of Campo deCaragena (299 km2 in an area of 1135 km2). We therefore assume
that 26% of the SGS estimated natural recharge originates from irrigated lands, 0.26  46 hm3 y1  12 hm3 y1. Thus irrigation return
ﬂow accounting for both precipitation and irrigation is
23 hm3 y1 + 12 hm3 y1 = 35 hm3 y1. We can divide approximately this return ﬂow among annual row crops, perennial vegetables, and fruit trees by partitioning it in proportion to theirrespective land areas: 15 hm3 y1 from annual row crops, 4 hm3 y1from perennial vegetables, and 16 hm3 y1 from fruit trees. These
totals are only approximate because it is unlikely that the different
cropping systems have identical water use efﬁciencies. Nevertheless, the approximate values are useful for evaluating our results.
Field site and experiment
A study of root zone soil moisture was conducted on an experimental plot at the Tomas Ferro Agricultural Science Center, a research facility operated by the Technical University of Cartagena.
The plot was managed according to agricultural practices that arecommon in the Campo de Cartagena region, including crop rotations (melon and lettuce) and drip irrigation. Meteorological data
for the site are available from a Servicio de Información Agrariade Murcia (SIAM, 2007) weather station located 235 m from theexperimental plot.
Characterization of the soil properties
An experimental plot measuring 8  3 m2 was established on asilty loam soil (USDA classiﬁcation system). To characterize the soil
physical properties, two soil cores were extracted from the plot toa depth of 2 m. The cores were sectioned into layers and analyzed
in the laboratory to determine soil bulk density (Grossman and
Reinsch, 2002), volumetric water content (Topp and Ferré, 2002),
and percentages of sand, silt, and clay (Gee and Or, 2002). Soilwater retention values at pressure heads greater than 500 cm
were measured using a pressure plate extractor (Dane and Hopmans, 2002).
Field experiment
A drip irrigation system, similar to that used in Campo deCartagena agriculture, was installed on the plot (Fig. 2). The system featured 16 mm inside diameter tubing, 4 L h1 emitters,
and an emitter spacing of 30 cm. In total, 36 emitters were
installed.
The plot was instrumented to monitor soil water dynamics inthe root zone (Fig. 2). Instrumentation consisted of two tensiometers (Soilmoisture Equipment Corp, Goleta, CA, USA) installed vertically at the depths 30, 45, 60, 90 and 120 cm (10 tensiometers
total) (Young and Sisson, 2002), and two 44 mm diameter, 2 m
deep access tubes for soil moisture measurements with a TRIMEFM TDR probe (Imko GmbH, Ettlingen, Germany) (Laurent et al.,2001, 2005).
Estimation of recharge using root zone modelling
Numerical modelling
We simulated water ﬂow and root water uptake using HYDRUS1D (Šimůnek et al., 2005). Assuming (i) that the soil is homogeneous and isotropic, (ii) that the air phase does not affect liquid

140

J. Jiménez-Martínez et al. / Journal of Hydrology 367 (2009) 138–149

ﬂow processes, and (iii) that water ﬂow due to thermal gradients is
negligible, the governing equation for water ﬂow is the 1D Richards equation:

 

@h
@
@h
¼
K
þ1 S
@t @x
@x

aðhÞ ¼

ð1Þ

where h = soil water pressure head (L); h = volumetric water content (L3 L3); t = time (T); x = vertical space coordinate (L);K = unsaturated hydraulic conductivity (L T1); and S = sink term
(L3 L3 T1), deﬁned as the volume of water removed from a unit
volume of soil per unit time due to plant water uptake. The sinkterm is speciﬁed in terms of a potential uptake rate and a stress factor (Feddes et al., 1978):SðhÞ ¼ aðhÞSp

ð2Þ
3

3

1

where Sp is the potential water uptake rate [L L T ] and a(h) is
the dimensionless water stress response function (0 6 a 6 1) thatprescribes the reduction in uptake that occurs due to drought stress.
For a(h), we used the functional form introduced by Feddes et al.
(1978):8 hh4;h4 < h  h3>>
> h3 h4>>>< 1;
h3 < h  h2
hh1
>
>
;
h2 < h  h1>
h2 h1>>>:0;h  h4 or h > h1

ð3Þ

where h1, h2, h3, and h4 are threshold parameters such that uptake isat the potential rate when the pressure head is between h2 and h3,
drops off linearly when h > h2 or h < h3, and becomes zero whenh < h4 or h > h1. Crop-speciﬁc values for these parameters (Table
1) were taken from the database contained in HYDRUS-1D (Šimůnek et al., 2005).
The soil hydraulic properties were modelled using the van
Genuchten–Mualem constitutive relationships (Mualem, 1976;
van Genuchten, 1980):(hðhÞ ¼r
hr þ ½1þjahshjhn 11=n


h<0

hs

h0

11=n g2
KðhÞ ¼ K s Sle f1  ½1  Sn=ðn1ÞeFig. 1. Location and geology of the Campo de Cartagena area in southeast Spain.

ð4Þ
ð5Þ

141

J. Jiménez-Martínez et al. / Journal of Hydrology 367 (2009) 138–149

Fig. 2. Photograph and schematic of the experiment plot and the instrumentation.

Croph1 (cm)h2 (cm)h3 (cm)h4 (cm)was used at the bottom, the latter condition being appropriate
due to fact that the water table was far below the root zone(Šimůnek et al., 2005).

Melon
Lettuce

10
10

25
25

400
500

8000
8000

Potential evapotranspiration and root growth

Table 1
Root water uptake reduction parameters.

where Se is effective saturation:

Se ¼

hðhÞ  hr
hs  hrand where hs = saturated water content; hr = residual water content;
Ks = saturated hydraulic conductivity; a = air entry parameter;n = pore size distribution parameter; and l = pore connectivity
parameter. The parameters a, n, and l are empirical coefﬁcients that
determine the shape of the hydraulic functions. To reduce the number of free parameters, we took l = 0.5, a common assumption which
is based on the work of Mualem (1976).
Running the model required specifying the hydraulic parameters hr, hs, a, n, Ks, and l. We estimated these parameters using
Rosetta (Schaap et al., 2001), a pedotransfer function model thatpredicts hydraulic parameters from soil texture and related data.
Rosetta contains a hierarchy of pedotransfer functions that canbe used depending upon available data. We predicted the hydraulicparameters using data for bulk density and percentages of sand,silt, and clay. The data and estimated parameters are given in Table
2. Reﬁnements to these parameter estimates were made subsequently based on model ﬁtting to a subset of the measured watercontent and pressure head data (details given below).
HYDRUS-1D (Šimůnek et al., 2005) uses the Galerkin ﬁnite element method to solve Eqs. (1)–(5). An atmospheric boundary condition (explained in the next section) was implemented at the soil
surface while a free drainage condition (unit hydraulic gradient)

Implementing the atmospheric boundary condition required
specifying daily irrigation and precipitation rates, as well as thepotential evaporation and transpiration rates. To determine evaporation and transpiration, we calculated a reference evapotranspiration ET0(t) using the Penman–Monteith method (e.g., Kashyap andPanda, 2001). The potential evapotranspiration ETp(t) was thengiven by (Allen et al., 1998):

ET p ðtÞ ¼ K c ðtÞ  ET 0 ðtÞ

ð6Þ

where ET0(t) was discretized in daily time steps and Kc(t) is acrop-speciﬁc coefﬁcient that characterizes plant water uptake andevaporation relative to the reference crop. Fig. 3 illustrates the time
variation of Kc(t) in terms of annual crop growth stage (the initial,
crop development, mid-season, and late season stages). Allenet al. (1998) provide data on the length of the growth stages andthe values of Kc for various crops. We used the Allen et al. (1998)
method and data to specify for each crop the value of Kc during eachgrowth stage (Table 3, Fig. 3).
With ETp given by Eq. (6), potential evaporation Ep(t) can be calculated according to (e.g., Kroes and Van Damm, 2003; Pachepsky
et al., 2004):

Ep ðtÞ ¼ ET p ðtÞ  expbLAIðtÞ
where b (0.4) is the radiation extinction coefﬁcient and LAI(t) is
the leaf area index. However, we lacked LAI(t) data so we instead
calculated

Ep ðtÞ ¼ ET p f ðtÞ

ð7Þ

Table 2
Measured soil textural and bulk density data, along with estimated hydraulic parameters.a
Depth (cm)

0–30
30–60
60–90
90–150a

Textural fractions
Sand (%)

Silt (%)

Clay (%)

18.7
13.8
19.5
10.8

76.0
80.2
77.2
82.0

3.5
6.0
3.3
6.6

Bulk density (g cm3)hr (cm cm3)hs (cm cm3)

Log10 (a) (cm1)

Log10 (n)Log10 (Ks) (cm d1)

1.45
1.52
1.58
1.70

0.04 ± 0.02
0.05 ± 0.03
0.04 ± 0.03
0.04 ± 0.03

0.38 ± 0.06
0.38 ± 0.06
0.35 ± 0.07
0.38 ± 0.08

2.18 ± 0.51
2.16 ± 0.51
2.01 ± 0.66
2.03 ± 0.76

0.21 ± 0.09
0.21 ± 0.09
0.19 ± 0.11
0.18 ± 0.11

1.67 ± 0.49
1.48 ± 0.52
1.48 ± 0.67
1.12 ± 0.80

Conﬁdence intervals are two standard deviations 95%, estimated by Rosetta.

142

J. Jiménez-Martínez et al. / Journal of Hydrology 367 (2009) 138–149

Results and discussion
Model calibration and predictions

Fig. 3. Illustration of crop growth stages and the time variation of the crop
coefﬁcient Kc.

where the function f(t) was speciﬁed based on the following reasoning. When a crop is ﬁrst planted, ground cover is nonexistent, potential evaporation is maximal, transpiration is zero, and thusf(t) = 1. Conversely, when the crop reaches the mid-season growthstage, ground cover is complete, evaporation is effectively zero,and thereafter f(t) = 0. All that remains is specifying the transition
from f(t) = 1 at planting to f(t) = 0 at the beginning of the mid-season growth stage. Because crop growth typically follows an Sshaped pattern (e.g., Overman and Scholtz, 2002), we modelled thistransition using a sigmoid curve.
With the HYDRUS atmospheric boundary condition, water
evaporates from the soil surface at the potential rate Ep (a ﬂux
boundary condition) as long as the pressure head at the surface remains above a threshold value, hcrit. If the soil surface dries out
such that the surface pressure head reaches the threshold value,
the boundary switches to a constant pressure head condition(=hcrit), generally leading to a computed actual evaporation ratethat is well below the potential rate. In our simulations, hcrit was
assumed to be 15 000 cm. The results of our simulations wereinsensitive to this parameter value when speciﬁed in the range
10 000 cm to 15 000 cm because the surface soil remained relatively wet due to regular irrigation and thus remained above thehcrit threshold.
With ETp and Ep given by Eqs. (6) and (7), the potential transpiration Tp(t) was speciﬁed by:

T p ðtÞ ¼ ET p ðtÞ  Ep ðtÞ

ð8Þ

Fig. 4 summarizes the imposed surface boundary condition, showing daily values of precipitation, irrigation, Tp, and Ep.
The modelled rooting depth was assumed to increase with time
according to a logistic growth function (Šimůnek et al., 2005),
achieving a maximum depth at the end of the crop development
stage. Values for the maximum rooting depth for particular crops(Table 3) were derived from Allen et al. (1998).

Running HYDRUS-1D using the Rosetta hydraulic parameter
estimates resulted in simulations that were in poor agreementwith the ﬁeld data. We therefore attempted to calibrate the soil
hydraulic property model Eqs. (4) and (5) using a subset of thedata. The most intensive data collection at the site occurred during
the cultivation of melon from 17 May to 10 September, 2007, comprising day of year (DOY) 137-253. We used this data set and theparameter optimization routines of HYDRUS-1D to calibrate thesoil hydraulic parameters. Several possible parameterizations wereconsidered which differed according to the number of soil layers(from 1 to 4) and the number and type of hydraulic parametersthat were ﬁtted for each layer (different combinations of 1–4
parameters among a, n, Ks, and hs). The initial estimates for theparameters when ﬁtted, or their ﬁxed value when not ﬁtted, were
the Rosetta estimates given in Table 2. If the results of the optimization looked promising, the ﬁtting was repeated using different
initial estimates to ensure that the same optimized parameter values were obtained. The best overall parameterization was determined informally based on diagnostic information provided by
the HYDRUS-1D routines about the model ﬁt and the convergence
behaviour of the inverse algorithm, visual inspection of the modelﬁt to the data (including laboratory water retention data), and theprinciple of parsimony (i.e., if two parameterizations produced aroughly equal ﬁt to the data, we took the simpler of the two, where‘‘simpler” means fewer ﬁtted parameters and/or soil layers). The
best parameterization was found to involve four soil layers with
two adjustable parameters, a and n, for each layer. The ﬁnal ﬁtted
hydraulic parameter values are given in Table 4; the other parameter values are those estimated with Rosetta (Table 2). Fitting morethan two parameters per layer tended to cause the inverse algorithm to fail. Four soil layers produced a better ﬁt to the data than
was possible with fewer layers in the proﬁle. Overall, the numerical solution with four layers and a and n ﬁtted for each layer provided the best correlation between measured and simulated watercontent and pressure head values. Fig. 5 shows the water content
and pressure head data for various depths in the soil proﬁle, alongwith the corresponding ﬁtted model simulation. Table 5 gives
goodness-of-ﬁt measures for the calibrated simulation run anddata, including R2 = 0.90 for the combined water content and pressure head data. Correlation coefﬁcients rx,y computed by HYDRUS1D for the eight ﬁtted parameters were |rx,y| < 0.4, with three
exceptions: ra1 ;a2 ¼ 0:54, r a1 ;a4 ¼ 0:55, and r a3 ;n3 ¼ 0:75 (where
numerical subscripts on a and n indicate the soil layer).
With the ﬁtted parameterization, HYDRUS-1D was next used
to predict the root zone soil moisture dynamics during the growing periods of two lettuce crops at the site, the ﬁrst one grown
between 7 February and 16 May 2007 (DOY 38-136) and the second one between 25 September and 24 December 2007 (DOY

Table 3Growth and evapotranspiration coefﬁcients for various crops (source: Allen et al., 1998).
Crop

Lettuce
Broccoli
Cauliﬂower
CeleryEndiveMelon
WatermelonaPlant date

January/September
January/September
January/September
January/September
January/September
May
May

Growth stage (number of days)

Crop coefﬁcient (Kc)

Maximum root depth (cm)

InitialDevelopment

Mid-season

Late-season

Initial

Mid

End

35/30
35
35
25

50/40
45
50
40

45/25
40
40
45

10/10
15
15
15

a

a

a

a

25
20

35
30

40
30

20
30

0.7
0.7
0.7
0.7
0.7
0.5
0.4

1
1.05
1.05
1.05
0.95
1.05
1

0.95
0.95
0.95
1
0.90
0.75
0.75

No data available, but casual observation suggests growth stages similar to those of lettuce.

30–50
45–60
45–60
45–60
30–45
80–15080–200J. Jiménez-Martínez et al. / Journal of Hydrology 367 (2009) 138–149

143

Fig. 4. Summary of the modelled soil surface boundary conditions (I = irrigation; P = precipitation; Ep = potential evaporation; Tp = potential transpiration).

Table 4Fitted hydraulic parameter values with 95% conﬁdence intervals.
Depth (cm)a (cm1)

n (–)

0–30
30–60
60–90
90–1500.078 ± 0.01
0.046 ± 0.005
0.014 ± 0.004
0.020 ± 0.002

1.16 ± 0.01
1.23 ± 0.02
1.27 ± 0.07
1.46 ± 0.05

268-358). Figs. 6 and 7 compare the predictions with the experimental data at various depths in the soil proﬁle. Overall, good
agreement was achieved between the ﬁeld measurements andthe HYDRUS-1D predictions. The correlation (R2) between measured and predicted water contents and pressure heads was
0.82 for both of the lettuce crops. Additional measures of goodness-of-ﬁt are given in Table 5.
Recharge estimation
We next used the calibrated model simulations to evaluate theroot zone soil water balance and calculate the annual recharge rate.
The components of the water balance were water additions due toirrigation (I) and precipitation (P), losses due to evapotranspirationand drainage below the root zone, and changes in root zone water
storage. Drainage from the bottom of the soil proﬁle was assumed
to be irrigation return ﬂow, equal to the groundwater recharge
rate. Based on the simulated drainage, we could determine the recharge rate generated during cultivation of different crops at different times of the year. Fig. 8 shows the cumulative drainage(recharge) as a function of time computed for the three cropping
periods. The percentages of applied water (P + I) becoming recharge were 24.6% for the ﬁrst lettuce crop (DOY 20-136), 21.7%
for the melon crop (DOY 137-253), and 68.2% for the second lettuce
crop (DOY 254-358). The higher percentage during the last period
was due to higher precipitation and lower transpiration rates(Fig. 4). In total, recharge during the period DOY 20-358 was
492 mm.
Fig. 8 also shows for the three cropping periods the cumulative
potential and actual transpiration rates. Actual transpiration was
frequently lower than the potential rate because soil moisture atvarious times was insufﬁcient to sustain the potential uptake rate,with uptake being reduced according to Eq. (2). This aspect of theroot zone modelling methodology (i.e., the physically based calculation of the onset of water stress and subsequent reduced uptake)is crucial for calculating recharge in arid and semi-arid regions
where ET may drop below potential rates, even in irrigatedsystems.
Vegetables, fruits, and citrus are the principal crops in the Campo de Cartagena area. In addition to the melon and lettuce crops
grown in our study, common fruits and vegetables include broccoli,cauliﬂower, celery, endive, and water melon. As shown in Table 3,these crops have similar cultivation periods and water requirements as the crops grown in our ﬁeld plot, and they are grown withthe same climate and soil conditions, which are relatively uniform
across the Campo de Cartagena region (Ramírez et al., 1999).
Therefore, as a ﬁrst approximation, we assumed our results formelon and lettuce to be representative of irrigation return ﬂows
for seasonal fruits and vegetables in the Campo de Cartagena plain.
This approach allowed us to use statistics for crop acreages in theCampo de Cartagena area (CARM, 2007) to estimate recharge to the
top unconﬁned aquifer from irrigated annual fruit and vegetable
croplands. Row crops in the region cover about 128.1 km2. Approximately 19.3% of the crops are cultivated between January and
April (lettuce, broccoli, cauliﬂower, celery, and endive), 25.3% arecultivated between May and August (sweet melon and water melon), while 55.4% are cultivated between September and December(lettuce, broccoli, cauliﬂower, celery, and endive). Using our recharge values for each of these periods, we estimate that recharge
to the top unconﬁned aquifer from irrigation of annual row crops
during 2007 was 26 hm3. This value does not include recharge arising from important irrigated perennial crops such as artichoke andcitrus.
Notably, our recharge value of 26 hm3 y1 for annual row
crops is almost two times greater than the 15 hm3 y1 estimatethat was derived above using data from the Spanish Geological
Survey (IGME, 1994). Although the SGS estimate encompassed
different estimation methods depending upon available data forthe various crops and cultivated areas, the value was based inpart on calculations which assumed crop water use was at thepotential rate. As indicated by Fig. 8, that assumption tends tooverestimate water use and hence underestimates recharge.
Based on these considerations and our calculation of 26 hm3 y1

144

J. Jiménez-Martínez et al. / Journal of Hydrology 367 (2009) 138–149

Fig. 5. Pressure head and water content data (points) measured at various depths in the soil proﬁle, along with ﬁnal ﬁtted HYDRUS simulations (solid lines) for melon.

for annual fruits and vegetables, we believe that that recharge
from irrigated farmland to the top unconﬁned aquifer of the Campo de Cartagena region is likely to be larger than previously
estimated.

Uncertainty and sensitivity assessment
Our modelling approach contains several potential sources of
uncertainty or error. We can distinguish between two levels or

145

J. Jiménez-Martínez et al. / Journal of Hydrology 367 (2009) 138–149
Table 5
Goodness-of-ﬁt measures for simulations and experimental data.
Simulation

Data set

RMSEa

MAEb

Melon (DOY 137-253) (calibration)

Water content (h)Pressure head (h)Combined h and h

0.029
33.7760.024
27.971

1st Lettuce (DOY 38-136) (prediction)

Water content (h)Pressure head (h)Combined h and h

0.024
56.6300.022
41.271

2nd Lettuce (DOY 268-358) (prediction)

Water content (h)Pressure head (h)Combined h and h

0.023
40.0130.019
31.279

R2

0.90

0.82

0.82

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
PRoot mean square error, RMSE ¼ 1n ni¼1 ðxi  yi Þ2 .
P
bMean absolute error, MAE ¼ 1n ni¼1 jxi  yi j.
atypes of uncertainty. On one hand, uncertainty exists about thecomputed root zone drainage at our plot or ﬁeld. In this context,
sources of uncertainty include the values of various model parameters such as the soil hydraulic parameters in Eqs. (4) and (5), thedaily reference evapotranspiration rate ET0(t), crop coefﬁcients
Kc(t), and drought stress parameters (e.g., h3 in Eq. (3), the threshold pressure head below which uptake is reduced). Quantifying theeffect of these parameter uncertainties on drainage calculations requires knowledge of their statistical variability and correlation
structure. In the case of the hydraulic parameters a and n, we haveinformation in the form of 95% conﬁdence intervals and correlation
coefﬁcients computed by HYDRUS-1D as part of the parameter
optimization (Table 4). For other hydraulic parameters, estimatedconﬁdence intervals are available from the Rosetta model (Table2), but those estimates are unrealistically broad because they arebased on only the soil separate and bulk density data and havenot been conditioned on the water content and pressure head ﬁeld
data. To our knowledge, the literature provides little or no information on quantifying uncertainty in other model parameters such as
crop coefﬁcients.
A second type or source of uncertainty involves extrapolation of
the ﬁeld results to the larger region. Here, sources of uncertainty include the data on regional irrigation practices, cropping rotations
and acreages, and so forth—the type of uncertainty that exists inany modelling or water balance estimation of recharge. Additionally, we may consider how speciﬁc model parameters vary in the region and how those variations would affect root zone drainage
calculations. Again, however, little knowledge exists on how toquantify such variability in parameters like ET0 or Kc(t). For thehydraulic parameters, the Rosetta uncertainty estimates (Table 2)potentially provide a good starting point for quantifying the effectsof hydraulic property uncertainty. Those estimates could, for example, serve as the basis for a Monte Carlo calculation of uncertainty inthe recharge estimate. However, such a calculation is more complicated than it may ﬁrst appear. The Rosetta uncertainty estimatesare quite broad. For example, the bounds for the 95% conﬁdence
intervals for Ks span an order-of-magnitude. With that kind of variability, a problem that may be encountered is that a given irrigation regime may not be sensible over the whole range of soil
properties: a realistic regime for a medium conductivity soil may
result in water logging in a low conductivity soil (a speciﬁc example
is given in the sensitivity calculations discussed below). Thus, aquantitative analysis of hydraulic property uncertainty would haveto somehow account for the relationship between soil properties
and irrigation regimes (i.e., surface boundary conditions).
In sum, it does not seem possible at this time to put quantitative conﬁdence intervals on our recharge estimate owing to a lackof knowledge about parameter variability and correlation structure at multiple levels. However, it is possible to get some senseof the importance of parameter uncertainty by performing a sen-sitivity analysis. For this analysis, we performed a series of simulations in which individual parameters, or in one case a set ofparameters, were perturbed a ﬁxed amount while all otherparameters were held at their baseline values (that is, the values
used in our recharge calculations). The effect of the various
parameter perturbations on the calculated recharge was thenevaluated. Parameters considered for sensitivity analyses werea, n, Ks, hs, hr, l, Kc(t), h3, and ET0(t). Note in the following that Rosetta assumes a, n, and Ks are lognormally distributed, such thatthe conﬁdence intervals are not symmetric about the (geometric)
mean parameter estimate after antilog transformation. Also notethat in our calibration procedure above we did not treat thehydraulic parameter l as adjustable, instead ﬁxing its value atl = 0.5; we include it in our sensitivity analysis for completenessand because Rosetta generally estimates high levels of uncertainty for this parameter.
The ﬁrst set of sensitivity calculations involved perturbing thehydraulic parameters a, n, Ks, hs, hr, and l one at a time in individual
soil layers. The perturbations to the parameter values in this casewere large, equal to the bounds of the 95% conﬁdence intervals given in Table 2. The results for these calculations showed that among
the parameters considered, recharge calculations were least sensitive to the water contents hs and hr. Setting either of these parameter values to the bounds of the 95% conﬁdence intervals in any of
the four soil layers resulted in a change to the computed rechargeof less than 3%, a small change considering that the 95% conﬁdencebounds corresponded to parameter perturbations of about ±15% to±23% for hs (depending on the soil layer), and ±54% to ±70% for hr(Table 2). Relatively low sensitivity was also found for the remaining soil hydraulic parameters (a, n, Ks, l) in the middle two soil layers. Setting Ks or l in those layers to the 95% conﬁdence bounds
resulted in changes to the computed recharge of less than 3%, while
equivalent perturbations to a or n produced changes less than 5%.
Again, the changes in computed recharge were fairly small considering the size of the parameter perturbations, which were about±25% for n, between 70% and +400% for Ks and a (recall the asymmetric conﬁdence intervals), and greater than 1000% for l. Higher
sensitivity was found for a, n, Ks, and l in the surface soil layer. Perturbations to those parameters resulted in changes to the computed
recharge that ranged from 3% to 12% with one notable exception:
the lower bound for the surface saturated conductivity resulted in
an increase of 42% to the calculated recharge. The reason for thislarge increase was that the low surface conductivity caused the surface soil to stay very close to saturation for long periods of time,such that uptake was reduced according to Eq. (3) and hence simulated drainage increased (the Feddes et al. (1978) model, Eq. (3),speciﬁes that uptake reduction occurs both when the soil is excessively wet and when it is too dry). The highest sensitivity to hydraulic parameters was found in the bottom soil layer where theretention and conductivity functions directly impact the drainage

146

J. Jiménez-Martínez et al. / Journal of Hydrology 367 (2009) 138–149

Fig. 6. Pressure head and water content data (points) measured at various depths in the soil proﬁle, along with HYDRUS predictions (solid lines) for the ﬁrst lettuce crop.

of water out of the root zone. In this layer, setting individual parameters to the lower and upper bounds of the 95% conﬁdence region
led to changes in recharge of, respectively, 21% and +18% for Ks,
6% and 1% for n, +13% and 5% for l, and +46% and 23% for a.

Although some of those changes are quite substantial, keep in mind
the very large perturbations (e.g., the bounds for a corresponded toperturbations of about 80% and +475%) and the low likelihood ofsuch a soil parameter value.

J. Jiménez-Martínez et al. / Journal of Hydrology 367 (2009) 138–149

147

Fig. 7. Pressure head and water content data (points) measured at various depths in the soil proﬁle, along with HYDRUS predictions (solid lines) for the second lettuce crop.

When narrower, more realistic bounds for a and n in individual
soil layers were evaluated (Table 4), the changes in the calculated
recharge were ±3% or less. Unlike the bounds discussed in the pre-vious paragraph, these conﬁdence intervals incorporated our measured water content and pressure head data, and thus are likely to
be more reﬂective of the uncertainty that existed at our experi-

148

J. Jiménez-Martínez et al. / Journal of Hydrology 367 (2009) 138–149

ing both numbers. The results of our sensitivity analysis do not ruleout the possibility of that level of uncertainty in our calculations,but they do suggest it is unlikely. Future work aimed at quantifying
uncertainty in parameters such as Kc(t) or ET0(t) would greatly beneﬁt efforts to determine uncertainty in recharge calculations.

Summary and conclusions

Fig. 8. Cumulative potential transpiration, actual transpiration, and drainage(recharge) rates computed for the three crops.

mental site. We additionally considered smaller, ±10% perturbations to the remaining hydraulic parameters (Ks, hs, hr, and l) inindividual soil layers. In these calculations, the calculated recharge
was least sensitive to hr ( 1% change in recharge), while the otherparameter perturbations altered recharge by between 7% and 15%.
The simulation for a 10% perturbation of n was not completed because the small n value caused numerical difﬁculties in the simulation model (a well-known problem for n ? 1). Lastly, we
considered simultaneously perturbing the same hydraulic parameter in all soil layers. Using ±10% perturbations, computed recharge
was again insensitive to hr (about 1% change in recharge), whereaschanges in recharge for other parameters were about 7–15%. The
lower bound calculation for n was again aborted due to numericalproblems.
The ﬁnal sensitivity calculations involved ±10% perturbations to
the parameters h3 and Kc(t). Note that because Kc(t) and ET0(t) appear
as a product in Eq. (6), results for a ±10% perturbation in Kc(t) are
identical to results that would be obtained for a ±10% perturbationin ET0(t). The recharge calculation was relatively insensitive to perturbations of h3, with the calculated recharge changing by less than
1%. Greater sensitivity was found for Kc(t) or ET0(t), where the 10%
perturbations altered the recharge calculation by about 6–7%.
In sum, an accurate characterization of prediction uncertainty is
problematic because of a lack of knowledge about the underlying
parameter variability and parameter correlation structure. Our
conclusion that irrigation return ﬂow in the Campo de Cartagena
region is likely larger than previously estimated is based on our
calculated value of 26 hm3 y1 recharge from annual row crops
and a previous estimate of 15 hm3 y1. In our view, that conclusion
would remain valid unless our estimate was about 25% too large, inwhich case the two recharge estimates would not be that different
considering the various approximations and assumptions underly-

In irrigated semi-arid and arid regions, accurate knowledge ofgroundwater recharge is important for the sustainable management of scarce water resources. The Campo de Cartagena plain insoutheast Spain is a semi-arid region where irrigation return ﬂow
accounts for a substantial portion of recharge. In this study we estimated irrigation return ﬂow using a root zone modelling approachin which irrigation, evapotranspiration, and soil moisture dynamics for speciﬁc crops and irrigation regimes were simulated with
HYDRUS-1D. The model was calibrated using ﬁeld data collected
at an experimental plot in the area. Good agreement was achieved
between the HYDRUS-1D simulations and ﬁeld measurements
made under melon and lettuce crops. The simulations showed thatwater use by the crops was below potential levels, despite regular
irrigation.
Assuming that the HYDRUS-1D results for melon and lettuce
were representative of similar crops grown in the Campo de Cartagena, we estimated that the cultivation of annual fruits and vegetables on irrigated farmland produces 26 hm3 y1 groundwater
recharge to the top unconﬁned aquifer. This estimate does not include important irrigated perennial crops in the region such as artichoke and citrus. Overall, the results point toward a greater
amount of irrigation return ﬂow in the Campo de Cartagena region
than was estimated previously by the Spanish Geological Survey.
Our calculations indicated a high level of recharge late in the year
when potential evapotranspiration was lower. Improved irrigation
scheduling based on soil moisture status and crop water requirements could signiﬁcantly reduce irrigation return ﬂows.
More generally, our root zone modelling approach shows promise as a method for estimating recharge in irrigated semi-arid regions. Compared with other techniques, the modelling approachis relatively data intensive, involving several of crop- and soil-speciﬁc parameters. However, often times these parameters can beapproximated using existing databases and estimation tools. For
example, maps of soil texture and related soil physical properties
are available for many locations; soil hydraulic parameters canbe easily estimated from those data using pedotransfer functions
or related approaches. Likewise, water stress parameters for many
important crops have been tabulated (e.g., Kroes and van Damm,
2003). Perhaps more difﬁcult will be specifying representative irrigation practices, which may or may not be uniform across a region.
Future work aimed at quantifying parameter uncertainty and correlation structure is needed to establish conﬁdence intervals formodel predictions.
Acknowledgements
This work has been developed within Projects CGL-200405963-C04-01 and CGL2007-66861-C04-03, national plan I+D+I
Spanish Ministry of Education. Gratitude is expressed to Almudena
Jiménez for collecting data and working with the instrumentation,and to the Technical University of Cartagena.

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, FAO, Rome, Italy.

J. Jiménez-Martínez et al. / Journal of Hydrology 367 (2009) 138–149
Allison, G.B., Gee, G.W., Tyler, S.W., 1994. Vadose-zone techniques for
estimating groundwater recharge in arid and semiarid regions. Soil Sci. Soc.
Am. J. 58, 6–14.
Brunner, P., Bauer, P., Eugster, M., Kinzelbach, W., 2004. Using remote sensing toregionalize local precipitation recharge rates obtained from the chloride
method. J. Hydrol. 294 (4), 241–250.
Campbell, G.S., Norman, J.M., 1998. An Introduction to Environmental Biophysics,
second ed. Springer-Verlag, New York.
CARM, 2007. Consejería de Agricultura y Agua de la Región de Murcia. Agrarian
Statistics Data. <http://www.carm.es>.
Dane, J.H., Hopmans, J.W., 2002. Pressure plate extractor. In: Dane, J., Topp, C. (Eds.),
Methods of Soil Analysis, Part 4, SSSA Book Series, vol. 5. Am. Soc. Agron.,
Madison, WI, pp. 688–690.
de Vries, J.J., Simmers, I., 2002. Groundwater recharge: an overview of processes andchallenges. Hydrogeol. J. 10, 5–17.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field Water Use and Crop
Yield. John Wiley and Sons, NY.
Finch, J., 1998. Estimating direct groundwater recharge using a simple waterbalance model – sensitivity to land surface parameters. J. Hydrol. 211, 112–125.
Gartuza-Payán, J., Shuttleworth, W.J., Encinas, D., McNeil, D.D., Stewart, J.B.,
DeBruin, H., Watts, C., 1998. Measurement and modelling evaporation for
irrigated crops in northwest Mexico. Hydrol. Process. 12, 1397–1418.
Gee, G.W., Hillel, D., 1988. Groundwater recharge in arid regions: review andcritique of estimation methods. Hydrol. Process. 2, 255–266.
Gee, G.W., Or, D., 2002. Particle size analysis. In: Dane, J., Topp, C. (Eds.), Methods ofSoil Analysis, Part 4, SSSA Book Series, vol. 5. Am. Soc. Agron., Madison, WI, pp.
255–294.
Grossman, R.B., Reinsch, T.G., 2002. Bulk density and linear extensibility. In: Dane, J.,
Topp, C. (Eds.), Methods of Soil Analysis, Part 4, SSSA Book Series, vol. 5. Am. Soc.
Agron., Madison, WI, pp. 201–228.
IGME, 1994. Las aguas subterráneas del Campo de Cartagena (Murcia). Spanish
Geological Survey, Spain. 62 pp.
Kashyap, P.S., Panda, R.K., 2001. Evaluation of evapotranspiration estimation
methods and development of crop-coefﬁcients for potato crop in a sub-humid
region. Agric. Water Manage. 50, 9–25.
Kendy, E., Gérard-Marchant, P., Walter, M.T., Zhang, Y., Liu, C., Steenhuis, T.S., 2003.
A soil–water-balance approach to quantify groundwater recharge from
irrigated cropland in the North China Plain. Hydrol. Process. 17, 2011–2031.
Kroes, J.G., Van Damm, J.C., 2003. Reference Manual SWAP: Version 3.0.3. Rep. 773.
Alterra Green World Res., Wageningen, The Netherlands.
Laurent, J.P., Ruelle, P., Delage, L., Bréda, N., Chanzy, A., Chevallier, C., 2001. On theuse of the TDR trime-tube system for proﬁling water content in soil. In:
Proceedings TDR’01, Evanston-Illinois, USA, 5–2 September 2001, pp. 1–10.

149

Laurent, J.P., Ruelle, P., Delage, L., Zaïri, A., Ben Nouna, B., Adjmi, T., 2005. Monitoringsoil water content proﬁles with a commercial TDR system: comparative ﬁeld
tests and laboratory calibration. Vadose Zone J. 4, 1030–1036.
Mualem, Y., 1976. A new model predicting the hydraulic conductivity ofunsaturated porous media. Water Resour. Res. 12, 513–522.
Overman, A.R., Scholtz, R.V., 2002. Mathematical Models of Crop Growth and Yield.
Marcel Dekker, Inc., New York.
Pachepsky, Y.A., Smettem, K.R.J., Vanderborght, J., Herbst, M., Vereecken, H.,
Wosten, J.H.M., 2004. Reality and ﬁction of models and data in soil hydrology.
In: Feddes, R.A. et al. (Eds.), Unsaturated-Zone Modeling. Kluwer Academic
Publishers, Dordrecht, the Netherlands.
Ragab, R., Finch, J., Harding, R., 1997. Estimation of groundwater recharge to chalkand sandstone aquifers using simple soil models. J. Hydrol. 190, 19–41.
Ramírez, I., Vicente, M., García, J., Vaquero, A., 1999. Mapa digital de suelos de laRegión de Murcia. Consejería de Agricultura, Agua y Medio Ambiente.
Handbook and CD-ROM, p. 78.
Sánchez, M.I., López, F., Del Amor, F., Torrecillas, A., 1989. La evaporación yevapotranspiración en el Campo de Cartagena y Vega Media del Segura.
Primeros resultados. Anales Edafología Agrobiol, 1239–1251.
Scanlon, B.R., Healy, R.W., Cook, P.G., 2002. Choosing appropriate techniques forquantifying groundwater recharge. Hydrogeol. J. 10, 18–39.
Schaap, M.G., Leij, F.J., van Genuchten, M.Th., 2001. ROSETTA: a computer programfor estimating soil hydraulic parameters with hierarchical pedotransfer
functions. J. Hydrol. 251, 163–176.
SIAM, 2007. Servicio de Información Agraria de Murcia. Climatology Data. <http://
www.siam.imida.es>.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2005. The HYDRUS-1D Software
Package for Simulating the Movement of Water, Heat, and Multiple Solutes in
Variability Saturated Media, Version 3.0. Department of Environmental Sciences
University of California Riverside, Riverside, California, USA, 270 pp.
Thornthwaite, C.W., 1948. An approach toward a rational classiﬁcation of climate.
Geogr. Rev. 38, 55–94.
Topp, G.C., Ferré, P.A., 2002. Thermogravimetric method using convective ovendrying. In: Dane, J., Topp, C. (Eds.), Methods of Soil Analysis, Part 4, SSSA Book
Series, vol. 5. Am. Soc. Agron., Madison, WI, pp. 422–424.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
Young, M.H., Sisson, J.Y., 2002. Tensiometry. In: Dane, J., Topp, C. (Eds.), Methods ofSoil Analysis, Part 4, SSSA Book Series, vol. 5. Am. Soc. Agron., Madison, WI, pp.
575–609.
Zhang, L., Dawes, W.R., Hatton, T.J., Reece, P.H., Beale, G.T., Packer, I., 1999.
Estimation of soil moisture and groundwater recharge using the TOPOG_IRM
model. Water Resour. Res. 35, 149–161.