Journal of Hydrology 538 (2016) 339–354

Contents lists available at ScienceDirect

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

Infiltration in layered loessial deposits: Revised numerical simulations
and recharge assessment
Elad Dafny a,⇑, Jirka Šimůnek bab

National Centre for Engineering in Agriculture, University of Southern Queensland, Toowoomba, QLD 4350, Australia
Department of Environmental Sciences, University of California, Riverside, CA, United Statesa r t i c l e

i n f oArticle history:
Received 15 February 2016
Received in revised form 11 April 2016
Accepted 14 April 2016
Available online 21 April 2016
This manuscript was handled by PeterK. Kitanidis, Editor-in-Chief, with theassistance of Juan V. Giraldez, Associate
Editor
Keywords:
Infiltration test
Groundwater recharge
Arid vegetation
HYDRUS
Sandy loess
Layered profiles u m m a r y
The objective of this study is to assess recharge rates and their timing under layered loessial deposits at
the edge of arid zones. Particularly, this study is focused on the case of the coastal plain of Israel and Gaza.
First, results of a large-scale field infiltration test were used to calibrate the van Genuchten parameters ofhydraulic properties of the loessial sediments using HYDRUS (2D/3D). Second, optimized soil hydraulicparameters were used by HYDRUS-1D to simulate the water balance of the sandy-loess sediments during
a 25-year period (1990–2015) for three environmental conditions: bare soil, and soil with both sparseand dense natural vegetation.
The best inverse parameter optimization run fitted the infiltration test data with the RMSE of 0.27 d(with respect to a moisture front arrival) and R2 of 96%. The calibrated model indicates that hydraulicconductivities of the two soil horizons, namely sandy loam and sandy clay loam, are 81 cm/d and17.5 cm/d, respectively. These values are significantly lower than those previously reported, based on
numerical simulations, for the same site.
HYDRUS-1D simulation of natural recharge under bare soil resulted in recharge estimates (to the aquifer) in the range of 21–93 mm/yr, with an average recharge of 63 mm/yr. Annual precipitation in thesame period varied between 100 and 300 mm/yr, with an average of 185 mm/yr. For semi-stabilized
dunes, with 26% of the soil surface covered by local shrub (Artemisia monosperma), the mean annual
recharge was 28 mm. For the stabilized landscape, with as much as 50% vegetation coverage, it was only
2–3 mm/yr. In other words, loessial sediments can either be a source of significant recharge, or of norecharge at all, depending on the degree of vegetative cover. Additionally, the time lag between specificrainy seasons and corresponding recharge events at a depth of 22 m, increased from 2.5 to 5 years, and toabout 20 years, respectively, with an increasing vegetative cover. For this reason, and also likely due to agreat depth of loessial sediments, no correlation was found between annual recharge and annual precipitations of the same year or subsequent years. Similarly, no differences were found between summer andwinter recharge fluxes. Instead, numerical simulations indicated continuous year-round recharge of theaquifer. We conclude that the layered subsurface acts as a short-term (annual) and long-term(multi-annual) buffer to smooth sudden precipitation/infiltration events. Vegetation conditions can help
in predicting long-term recharge rates (as percentage of annual precipitation), which in turn need to beconsidered when assigning recharge characteristics in regional assessments and models.
Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction
Recharge estimations play an important role in the management and research of groundwater systems. The magnitude andtiming of groundwater recharge with respect to the corresponding
infiltration events are controlled by climatic and geological factors,
as well as the depth to the water table, and have long been of scientific and practical interest (e.g., Wu et al., 1996; Scanlon et al.,⇑ Corresponding author.
http://dx.doi.org/10.1016/j.jhydrol.2016.04.029
0022-1694/Ó 2016 Elsevier B.V. All rights reserved.

2006). Recharge can be estimated by various methods, whichshould be selected based on the climate zone (arid/humid),
expected fluxes, the spatial scale of the aquifer, and the studied
time scale (Scanlon et al., 2002a). Such a task should also consider
the data availability in terms of soil hydraulic properties and climatic data (e.g., precipitation and evaporation time series).
Physically based models, such as those solving the Richards’
equation for water flow in the vadose zone and water balance ofsurficial sediments, have often been used for estimations ofgroundwater recharge under various conditions (Wu et al., 1996;

340

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

Scanlon et al., 2002b; Kurtzman and Scanlon, 2011; Leterme et al.,2012; Turkeltaub et al., 2015). These models require the knowledge of soil water retention and unsaturated hydraulic conductivity functions for all soil horizons. Nevertheless, as such knowledge
exists (or is measured) only in rare cases, inverse modeling, inwhich numerical models are calibrated using measured flowrelated variables, has gained popularity for estimating the hydraulic functions characterizing the unsaturated zone (Hopmans et al.,
2002; Vrugt et al., 2008). Inverse models are then tested by theirability to reproduce independent measurements (Jacques et al.,
2002; Ritter et al., 2003; Turkeltaub et al., 2015). Inverse modelinghas several important advantages over other methods of estimating soil hydraulic properties. These mainly include (a) the ability
to simultaneously estimate both the soil water retention andunsaturated hydraulic conductivity function from a single transient experiment, and (b) the fact that similar numerical models
are used both for estimation of soil hydraulic properties (calibration) and for predictive forward modeling (assessment of recharge)(e.g., Hopmans et al., 2002).
Loessial soils are regarded by many as transmissive andrechargeable (Lin and Wei, 2006; Seiler and Gat, 2007; Aish,2014). Nevertheless under some vegetation and climate conditions
they might ‘‘act as an impervious barrier to recharge” (Weinthal
et al., 2005). The objective of this paper is to provide a physically
based assessment of recharge rates under loessial soils for arid climate and vegetation. We use the case of the coastal aquifer bounding Israel and the Gaza strip. Large parts of this aquifer are coveredby loess and loessial sediments.
The thorough data set collected during a large infiltration test in
this area, as reported by Gvirtzman et al. (2008), is first used toquantify the hydraulic conductivities and van Genuchten retention
parameters of the loessial sediments of the studied area. Second,
calibrated soil hydraulic properties and long-term meteorological
data are used to assess the long-term natural recharge to theunderlying sandy aquifer. These two modeling stages are sequentially discussed in separate sections. Finally, insights obtainedusing this recharge model are discussed in relation to previously
used models and empirical recharge estimations of the same aquifer. The conclusions obtained in this study for the coastal plain
underlying Israel and Gaza Strip are relevant for other phreaticaquifers overlain by a thick, layered unsaturated-zone, especiallythose within a similar arid climate.

2. Study area
While it is also relevant for other areas on the edge ofarid-climate zones, the study focuses on the area covered by sandy
loess soils south of the Gaza Strip (Fig. 1). In this area, a relatively
thick (>20 m) layer of loessial deposits covers permeable sandy
layers of the coastal aquifer (Fig. 1b). The water table is at a depthof about 25–35 m, and the hydraulic gradient is toward theMediterranean coast, i.e., to the NW.
The term ‘loessial deposits’ is used here to describe the sandy
grain-size loess, according to the Yaalon and Dan (1974) definition,and to distinguish it from silty grain-size loess. The sediments
accumulated during the Late Pleistocene period when desert
storms brought sands from the northern Sinai Desert to the area
(Yaalon and Dan, 1974; Crouvi et al., 2010). Since the distance fromthe parent material to the accumulation site is relatively short,these wind-blown sediments maintained a fairly coarse-graincomposition, contrary to the classic silt-size loess sediments(Crouvi et al., 2010). Nevertheless, the grain size distribution in
the loessial sediments is highly variable and depends on the windspeed (Crouvi et al., 2010). During periods of intense winds, active

eolian abrasion of the sand grains increased, silt-size sediments
were imported, and vice versa. Gradually, a layered sandy loessdeposit was accumulated.
The area is dominated by the Mediterranean climate, with dry
and hot summers and a rainy season between October and April.
The average annual precipitation is about 200 mm/yr. At En Habesor meteorological station (#144872, Israel’s Meteorological Service (IMS) database, see location in (Fig. 1a)), an average
precipitation of 185 mm/yr was recorded during the last 25 years(1990–2015). The cumulative annual potential evaporation
exceeds 1800 mm.
The dominant species in the area is a desert dwarf shrub, Artemisia monosperma (Bar Kutiel et al., 2016). It has an important rolein the dune stabilization process, and covers up to 16%, 16–36%,
and 36–65% of the total area of mobile, semi-stabilized and stabilized dunes, respectively. A. monosperma develops a unique rootsystem in response to the accumulating or dispersing sands around
it. Similar to other species of dune vegetation, its roots seldom
exceed a depth of 3 m, but extend laterally to accommodate any
available water from the ‘open matrix’ in-between neighboring
shrubs (Bar Kutiel et al., 2016).
The observed run-off from sandy loess soils in the region is verysmall (Yair, 1990; Givati and Atzmon, 2009; Eshtawi et al., 2015).
For the entire catchment of the ephemeral Besor River, which iscovered in large parts by sandy loess sediments, the run-off coefficient for 1985/6–2009/10, defined as % of annual precipitation, wasonly 1.3% (Givati and Atzmon, 2009). For comparison, the average
run-off coefficient over the entire western drainage system ofIsrael in the same period was 4.4%.

2.1. 2004 Infiltration test
A large trench infiltration test was conducted within the study
area in 2004 by Gvirtzman et al. (2008). We only provide here ashort description of the infiltration test, its set-up, the measurements made, and the obtained results that are relevant to the current study. Detailed information about the experiment can befound in Gvirtzman et al. (2008).
An elongated trench (3 m wide by 17 m long at the base) was
excavated in the loessial section, and its sloping edges were covered by PVC sheets. Four bores (marked ‘A’ through ‘D’) were
drilled in the vicinity of the trench midline (Fig. 2). Each was
equipped with 7–8 TDR probes at different depths, which monitored the temporal changes in the water saturation. The probes
were numbered from top to bottom (i.e., probe C2 is deeper thanprobe C1, etc.). Due to technical issues, only data of 22 sensors
was collected, as noted in Fig. 2.
Soil and sediment samples taken from these bores allow for
detailed characterization of the profile. Soil moisture, bulk density,
porosity, content of fines (<0.075 mm), and other parameters were
measured every 2 m (Gvirtzman et al., 2008; Hatzor et al., 2009).
Overall, the loessial sediment profile is stratified, with alternating
horizons classified as silty sand, clayey sand, and low plasticity clay(‘SM’, ‘SC’, and ‘CL’, respectively) according to the Unified Soil Classification System (USCS). The exact description of the sedimentprofile, as well as soil physical properties, slightly varied from bore
to bore.
In March 2004, following the end of the winter season, thetrench was flooded with 1 m of water for 17 days. TDR readings
were taken in decreasing frequencies during 20 days of the experiment. The analysis of water content measurements shows that thewetting front propagates below the trench in an ‘onion-shape’ pattern. The data further indicates that the wetting front reached thelowest probe (C8), located 20 m below the trench base, between

341

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

+100

A’
SSE
Israel

Shoreline

Gaza Strip

A
NNW

Elevation [m ASL]

+50

0

A

Coastal aquifer

B1
B2

-50

Saqiye aquiclude

C
-100
0

2

4

6

8

10

12

14

16 km

Fig. 1. Location map and geological cross-section (after Tulmatz, 1977).

days 17 and 20 (i.e., at 8:00 a.m. on the 17th day, the sensor was
still dry, and by 8:00 a.m. on the 20th day, it was already wet).
The overall water volume used to fill the trench in this experimentwas 1400 m3. Inevitably, a small portion of water, which filled thetrench, directly evaporated.
A second experiment was carried out about three months later
in July 2004. For this experiment, the trench was lengthened to48 m. The sides of the trench were re-sheeted with PVC sheets.
The trench was then flooded with 1.25 m of water for 25 days.
TDR readings were taken less and less frequently during the 4th,
7th, 10th and 14th days of the experiment. The wetting front propagation during the second experiment was noted to be substantially faster than its propagation during the first experiment,
though the results were not reported in the paper. The overallwater usage in this experiment was 10,300 m3.

3. Inverse calibration of the first infiltration test
The results of the first infiltration test served to calibrate thehydraulic properties of the loessial sediments. The high uncertainties and generally less meticulous data for the second infiltration
test prevented it from serving as a validation test. Among otherissues, the water use (and percolation rate) in the second infiltration experiment was substantially higher than expected. Had theinfiltration rate been the same as in the first test, the water use
would
have
been
5813 m3(1400 m3  (25 d/17 d) 
(48 m/17 m)). Furthermore, had the infiltration rate increased
due to the higher head in the trench, the water use would havebeen 7266 m3 (5813  1.25/1.0). However, the reported water
use was 41–71% higher than expected. This difference cannot becaused by variations in the initial conditions at the onset of the342

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354
Table 1
Vertical discretization of the 2D domain.

Trench
base

C1

B1
A1
A2
A3
A4
A5

D1

C2

B2
B3
B4

C3

D3

C4

D4

B5

D5

B6
C6

D7

C8
x
z

Fig. 2. 2D domain showing the bores (black lines), finite element mesh (blue lines),and observation points (red). The diagram is 30 m high and 35 m wide. (Forinterpretation of the references to color in this figure legend, the reader is referredto the web version of this article.)

two experiments (as also demonstrated by our numerical results,
simulations not included here). Direct evaporation from the trench,in the order of 7.6 mm/d during July (see Section 4.1.2) could
explain only a fraction of the difference (only 47 m3). We hypothesize (partly based on a comparison of photos included in
Gvirtzman et al. (2008) and Hatzor et al. (2009)) that the trenchcover was less durable (hence less effective) during the second test,
thus increasing the effective percolation area. Alternatively, there
may have been technical issues with the water level or water use
metering.

3.1. Methods
3.1.1. Model settings
Numerical simulations of the infiltration experiment were conducted using version 2.03.0450 of the HYDRUS (2D/3D) code(Šimůnek et al., 2008). The code solves the two-dimensional (2D)
form of the Richards’ equation while allowing for time-variableboundary conditions. In order for the results to be comparable tothe previous study, a 2D domain resembling the one describedby Gvirtzman et al. (2008) was considered in numerical simulations. The domain was set perpendicular to the flooded trench,
assuming that a 3D flow pattern can be approximated using theflow field in the 2D x–z plane (where x and z are the horizontal
and vertical coordinates). This approximation seems to be justified,as the length of the trench (17 m, in the y direction) is significantly
longer than its width (3 m; in the x direction). Regardless, the 2D
domain is adequate to simulate the flow field in the central partof the trench, and it represents recharge through a one length unit.
The 2D transport domain is 35 m wide and 30 m deep, with a
2.5 m-deep trench located in the middle of its upper surface(Fig. 2). The total area of the transport domain is 1032.5 m2. The
profile is comprised of two soil units that form a sequence of fivehorizons (Table 1). Due to the heterogeneity between different
bores, no attempt was made to further sub-divide the subsurface
into distinct units. The locations of TDR probes were explicitly
introduced into the domain prior to generating the unstructured
finite element grid. The average cell size for the entire transport

HorizonElevation [m]Soil classification

No. of observation points

1
2
3
4
5

30–24.5
24.5–18.5
18.5–16.5
16.5–14.5
14.5–0.0

Sandy loam
Sandy clay loam
Sandy loam
Sandy clay loam
Sandy loam

3
11
3
1
4

domain was set to 20 cm, with point refinements (of 5 cm) at bothedges of the trench bottom and coarsening (of 50 cm) along thebottom part of the domain. The final mesh consists of 11,783 finite
elements, substantially less than the 24,174 elements of the structured grid used by Gvirtzman et al. (2008) in their model.
The trench base (3 m wide) was set as a variable head boundary(HYDRUS specific technical terms are in italics) to represent ponding. The trench was assumed to fill-up within 6 h to a final pressure
head of 100 cm. Following the end of the infiltration test, it drained
at a rate of the effective saturated hydraulic conductivity, i.e., forabout three days. The bottom of the domain was assumed to beunder a unit gradient condition and was set accordingly as a freedrainage boundary, as it is located far above the water table. The
left- and right-hand sides, as well as the top of the domain (except
for the trench topography), were assumed to be no-flow
boundaries.
The initial pressure head profile was set according to measured
moisture contents of the profile using the soil’s saturation-pressurefunctions. Overall, the pressure head alternates between 490 cm
and 410 cm, for the more permeable and less permeable horizons, respectively. Note that higher values of the pressure head(300 cm) were used in the previous model (Gvirtzman et al.,2008) throughout the entire profile.
3.1.2. Calibration strategy
The soil hydraulic properties were calibrated in a series of
inverse simulations for the first infiltration test (days 0–21). The
van Genuchten soil hydraulic parameters of both soil layers wereinversely optimized by the HYDRUS-2D software using the internal
Marquardt–Levenberg-type (Marquardt, 1963) parameter optimization algorithm. For the inverse calibration, time series ofwater contents at each monitoring point (sensor) were reconstructed based on the available data from the first infiltration
experiment. These time series included the timing of the first reading when the probe was found to be wet and the last reading when
it was still dry. Each time series thus contained, at least, threeobservations. The initial, low water content, which was assumed
to be constant until the arrival of the moisture front, was used initially (t = 0 d) and at the last time that the data logging suggested
that the sensor was dry. A relatively high water content (90% of
the saturated water content of a particular soil) was used at thenext data logging, after the arrival of the wetting front. For example, for the C8 sensor, h = 0.054 was set at t = 0 and 16.8 d andh = 0.33 was set at t = 19.8 d. For three probes in the uppermost
part of the profile (B1, C1, and C2), the data included discrete readings of water contents 4–5 times during the experiment. Thesesensors showed a gradual transition from dry (h = 0.05) to wet(h = 0.40) conditions, over 5 h to 4 days.
Additionally, the net cumulative volume of water that infiltrated to the subsurface during the experiment (total cumulative
influx minus evaporation) was used to constrain the solution.
The water volume used in the experiment was first divided by
the length of the trench, and then potential evaporation during
the experimental period (Section 4.1.2) was subtracted to give
cumulative infiltration flux (81.6 m3/m length), which was usedduring inverse optimization.

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

The inverse parameter estimation was conducted by minimizing the objective function (Šimůnek et al., 2009):

u ¼ vhnX2
wi ½hobs ðti Þ  hsim ðt i Þ2 þ v cf wcf ½cfobs  cfsim 

ð1Þ

i¼1

where n is the number of water content observations, hobs and hsim
are measured and simulated volumetric water contents at times ti,
wi is the weight for individual observations, wcf is the weight for the
cumulative infiltration, and cfobs and cfsim are measured and simulated cumulative infiltrations at the end of the experiment. Due torelatively large uncertainties involved, the observations weights(wi and wcf) were set to 0.1, to allow for a wider confidence range.
vh and vcf are the model-assigned weighting coefficients, which
minimize differences in weighting between the water content data
set and the cumulative flux data set.
The inverse calibration was carried out to estimate at once six
soil hydraulic parameters (three for each soil unit): the saturated
hydraulic conductivity Ks, and the a and n parameters (which control the shape of the retention curve). The residual and saturated
water contents (hr and hs, respectively) were set to representativevalues obtained in the field and lab (Hatzor et al., 2009) (Table 2).
The pore connectivity and tortuosity factor l, which showed to be
unimportant and/or insensitive in numerical experiments (e.g.,
Wang et al., 2003; Turkeltaub et al., 2015), was fixed at 0.5 for bothsoil units (Mualem, 1976).
Uniqueness and stability of the inverse solution were tested byrepeatedly solving the inverse problem using six different combinations of the observations while setting various initial estimates
for the calibrated parameters. The list of optimized parametersand different combinations of the observations is given in Table 4,
which also gives the final optimized values (discussed below).
3.2. Results
The simulated instantaneous infiltration rate through thetrench floor at the beginning of the first test was very highTable 2
Fixed (not calibrated) soil hydraulic parameters.
Soil unit

hr [ ]

hs [ ]l []Sandy loam
Sandy clay loam

0.035
0.085

0.38
0.36

0.5
0.5

343

(13.2 mm/min). The infiltration rates exponentially fell until stabilizing at a much lower level of about 0.9 mm/min, which is equivalent to 3.8 m3/m length/d. The cumulative infiltration during thefirst test was 81.6 m3/m length. Integration of the simulated infiltration rate during the first flood over the 17 m trench resulted
in a total of 1387 m3, with a perfect match of the measured influx(accounting for the calculated direct evaporation of 13 m3). During
this period, no recharge occurred under the experiment site, i.e.,
almost all infiltrated water remained (temporarily) in the unsaturated zone.
The simulated wetting front propagation for the first infiltration
test satisfactorily matched the conceptualized (Gvirtzman et al.,2008) wetting front propagation, as demonstrated in Fig. 3. ThoughFig. 3a may suggest some mismatch in the early stages of the simulation, careful examination of the first arrival times shows thatthe simulated results matched at least 3 out of 6 early observationsof saturation (i.e., C1, C3, D1). Using the considered conceptual setup, it was not possible to further reduce discrepancies between
measured and simulated arrivals of the wetting front at the otherthree observation points (B1, B2, and C2), as well as at five otherpoints (Fig. 4). Similar discrepancies were described in earlier studies (Wierenga et al., 1991; Zou et al., 2001; Gvirtzman et al., 2008;
Turkeltaub et al., 2015) and were attributed to the inherent heterogeneity of the natural soil material, including the spatial variabilityof hydraulic conductivities and the uncertain thickness of different
soil horizons. Thus, a better fit for each and every measurement
point was not expected nor pursued. Regardless, a very good match
can be seen from day 13 onward (Fig. 3b and c).
It should be noted that all simulated wetting front arrivals atbore B were earlier than observed. For that reason, we conducted
6 runs using and/or excluding various parts of the dataset. For
example, information collected in bore B was excluded in run #2,
and information collected in the top 10 m of the profile was
excluded in run #6. Table 3 indicates that all calibration runs produced the same narrow range of hydraulic conductivities for bothsoil units, regardless of what observation sets were used. Table 4summarizes the statistical results of the runs. Excluding bore B
observations did not change the calibration results (compare runs
#1 and #2 in Table 3), although it slightly improved the statisticalfit. This means that it is unlikely that the calibration is biased bydata from one or more bores (particularly bore B) or from the depth
at which they were taken.

Fig. 3. Simulated (colored) and conceptualized (background) water content profiles at three times during the first infiltration experiment: (a) day 4, (b) day 13, and (c) day 21.
Each diagram is 30 m high and 35 m wide. The colored contours depict water content of 19% (blue) and 29% (red). Background image is Fig. 6 of Gvirtzman et al. (2008). (Forinterpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)344

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354
Table 4
‘Goodness of calibration’ statistics for different inverse runs.
Run

Observation sets

Observation
nodes

Min
error
[d]

Max
error
[d]

RMSE
[d]

R2a
(%)

1
2

All bores
All bores except
B
Bores A and D
Bores B and C
Only probes in
top 10 m
Only probes
below top 10 m

22
16

0
0

3.2
0.7

1.14
0.27

91
96

10
12
10

0
0
0

1.2
3
3

0.60
1.43
1.06

92
93
71

12

0

1.5

0.95

90

3
4
5
6
a

Fig. 4. Simulated and observed (ranges) arrival times of the wetting front at various
probes (bores (A–D and probe numbers are shown in the upper row).

Based on these inverse optimization runs, the hydraulic conductivities of the two soil units were set to K1 = 80.2 cm/d andK2 = 17.3 cm/d. These calibrated values are substantially lower
than those previously proposed based on numerical simulations(i.e., 1103 cm/d and 525 cm/d, respectively, Gvirtzman et al.
(2008)). The previous, higher values would generally characterize
sandy soils, rather than loamy soils (Radcliffe, 2002). For example,
calibrated hydraulic conductivities for five sandy soils at a nearbysite range between 1069 and 1756 cm/d, while that of a sandy
loam soil was only 41 cm/d (Turkeltaub et al., 2015).
3.3. Soil classification
The calibrated saturated hydraulic conductivities for the twosoil units were K1 = 80.2 cm/d and K2 = 17.3 cm/d. These values
are in the range of conductivity values for ‘sandy loam’ and ‘sandy
clay loam’ soils, respectively (Radcliffe and Šimůnek, 2010). A further analysis of the grain-size distribution using the USDA soil textural triangle (Hillel, 2004) rather than the USCS’, justifies theabove soil classifications. Soil samples taken from the highconductivity layers contain 16–38% fines (i.e., grain smaller than
0.075 mm, as reported by Gvirtzman et al., 2008), with a medianR2 varies for each run, depending on which data were considered in calibration.

of 31% fines. Therefore, the proportion of sands (including the finesands category) in this soil unit is in excess of 70% and can thus beclassified as either sandy loam or sandy clay loam. Similarly, thelow-conductivity layers containing 57–68% fines, with a medianof 60% fines, can be classified as silt loam, loam, clay loam, or sandy
clay loam. Due to a lack of a further analysis of grain-size distribution, besides a general knowledge of the relatively low content ofclay in the area’s loess (Eshtawi et al., 2015), the final classificationof the two soils was made based on the measured and calibrated
van Genuchten parameters. The high sand content layer was
denoted ‘sandy loam’ and the low sand content unit was denoted‘sandy clay loam’ (Tables 1 and 2). It should be noted that we foundthe use of the USCS term ‘low plasticity clay’ for the latter soil unit
misleading.
The behavior of the entire system can be perceived as slightly
anisotropic, though not due to an initial explicit parameterization.
Simple calculations (Fetter, 2001) for the upper 22 m section givean effective vertical hydraulic conductivity (Kv,eff) of about
35 cm/d and an effective horizontal conductivity (Kh,eff) of
62 cm/d. The anisotropy is thus about 2 (i.e., the effective horizontal hydraulic conductivity is twice as large as the effective verticalone), merely due to soil layering. Similar low anisotropy values insandy loam soils, especially following long periods of no-tillage,
were described in previous studies (e.g., Petersen et al., 2008).
The low implicit anisotropy of the loessial sediments is in contrast
with the much higher, conventional, 1:10 explicit anisotropy of the
sandy layers of the Israeli coastal aquifer (Amir et al., 2013).
4. Natural recharge simulations
At the second stage, HYDRUS-1D (Šimůnek et al., 2009) wasused to assess natural recharge over a 25-year period throughloessial sediments similar to those found at the infiltration testsite. These recharge simulations aim to explore the overall water

Table 3
Calibrated van Genuchten (1980) parameters and their 95% confidence intervals for two soil materials (subscripts 1 and 2) and different inverse runs.
Run

Observation sets

K1 [cm/d]n1 [–]

a1 [1/cm]

K2 [cm/d]n2 [–]

a2 [1/cm]

1

All bores2

All bores except B

3

Bores A and D

4

Bores B and C

5

Only probes in top 10 m

6

Only probes below top 10 m

80.13
(79.98–80.27)
81.34(80.27–81.40)
80.11(80.11)
81.34(81.23–81.45)
80.38(79.99–80.78)
81.34(79.93–82.75)
80.77(79.93–82.75)

2.68(2.68)
2.72(2.70–2.74)
2.68(2.41–2.95)
2.72(2.71–2.73)
2.73(2.73)
2.72(2.57–2.88)
2.71(2.41–2.95)

0.0079(0.0079)
0.0078(0.00779–0.00785)
0.0079(0.0072–0.0086)
0.0078(0.0077–0.0078)
0.0078(0.0078)
0.0078(0.0076–0.0080)
0.0078(0.0072–0.0086)

17.35(17.35)17.14(17.13–17.15)
17.35(17.35)17.14(17.14)17.46(17.20–17.73)
17.14(17.14)17.26(17.13–17.73)

3.06(3.06)
2.82(2.78–2.86)
3.05(3.05)
2.82(2.77–2.86)
2.84(2.84)
2.82(2.82)
2.90(2.77–3.06)

0.0058(0.0058)
0.0070(0.0070)0.0058(0.0049–0.0067)
0.0070(0.0066–0.0073)
0.0068(0.0055–0.0082)
0.0070(0.0058–0.0082)
0.0066(0.0049–0.0082)

Statistics: Average range

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

balance of the vadose zone, including its main six components –precipitation (rainfall), evaporation, infiltration, runoff, recharge,
and the change in soil water storage. The simulations have been
carried out for three scenarios involving three different natural settings, namely, bare soil, sparse vegetation upon semi-stabilized
dunes, and denser vegetation upon a stabilized landscape. These
three settings represent successive steps in the stabilization of atypical loessial landscape at the fringe between the arid Sinai
Peninsula and the Mediterranean climate (Bar Kutiel et al., 2016).
4.1. Methods
4.1.1. Model settings
The sediment profile was considered to be 22 m deep, i.e., tomatch the maximal depth from which soil data was obtained during the infiltration test. It was discretized using 150 finite elementsof varying sizes, from 0.29 cm at the top to 29 cm at the bottom of
the profile. The sediment profile had the same layout (depths ofsoil horizons and their hydraulic parameters) as the 2D model ofthe infiltration site. Eight additional realizations, with a combination of ±10% of hydraulic conductivities for the two soil units, werealso run to test the sensitivity of the results. The ±10% variations
exceed the 95% confidence intervals for the hydraulic conductivities of the two soil materials, as obtained in the simulations of
the infiltration test (±3% for both). Five observation points thatenabled the recording of changing pressure heads, water fluxes,
and water contents, were defined at the bottom of every soil
horizon.
An atmospheric boundary condition with surface run-off was
assigned to the upper boundary of the transport domain. Thisboundary condition allows water to infiltrate, evaporate, and/or
discharge as surface run-off depending on prescribed precipitationand potential evaporation rates, as well as on the transient soil
conditions (moisture). This boundary condition requires the setting of a single empirical parameter, hCritA, which determines
the minimum allowed pressure head at the soil surface and separates the first and second stage of evaporation. This parameter
was assigned a value of 15,000 cm, which is the default value
for relatively coarse-grain soils. Multiple runs with lower values,
down to 100,000 cm, produced similar results, indicating thatmodeling results (in the scope of this paper) are not sensitive tovariation in hCritA. Precipitation and potential evaporation fluxes
are described in the following section. The bottom of the sedimentprofile was assumed to be under a unit gradient condition, and wasset accordingly as a free drainage boundary because it is still located
far above the water table.
The initial conditions for each of the three simulations wereobtained following initial runs with the proper vegetation parameters (Section 4.1.3) and 100 repetitions of the 1989/90 meteorological conditions (Section 4.1.2).
4.1.2. Meteorological data and temporal variables
Two meteorological data sets were obtained from the Israeli
Meteorological Service (IMS) website in order to represent local
conditions in the simulations of the recharge to the aquifer: (1)time series of daily precipitation from 1/10/1990 until
30/9/2014; and (2) time series of daily pan evaporation from
1/10/1990 to 31/12/1999 (no evaporation data were recorded
thereafter). Both sets were recorded at the En Habesor station,
which is located in the middle of the study area (Fig. 1).
The precipitation time series (Fig. 5a) was comprised of daily
values. To approximate higher intensity rainfalls than the daily values that occur in this area, each time record was further divided
into a short period (0.1 d), which was given the entire daily precipitation volume and a longer period (0.9 d) without precipitation.
For example, a day with a total precipitation of 12 mm would be345

divided into a 2.4-h period with 120 mm/d precipitation(120 mm/d⁄ 0.1 d = 12 mm) and a 21.6-h period without precipitation. This is a more realistic conceptualization of the local weather
patterns, and in turn, resulted in higher and more realistic rain
intensities (for the example above it is 5 mm/h vs. 0.5 mm/h,
respectively) and run-off.
The pan evaporation dataset included some records (about 15%of the records each year) that were obtained over a period of several days rather than for every single day. For this reason, theobtained data could not be used directly. Rather, a synthetic annual
time series was prepared based on the collected data, as follows.
Overall, daily evaporation varied between about 1 and 9 mm/day,without any significant differences between years (Fig. 5b). The
lowest values were recorded each year during January and thehighest values during July. Accordingly, a sinusoidal function was
fitted to observed daily evaporation values:

EV d ¼ sin

p



DOY
 EV fluc þ EV min
365

ð2Þ

where EVd is daily pan evaporation, DOY is the day within the calendar year (e.g., 1st of October of each year is DOY = 273), EVmin is minimal daily evaporation on January 1st, set to be 0.8 mm/d, and EVflucis the amplitude, at which the daily evaporation fluctuates throughout the year. EVfluc was optimized to fit the average annual potential
evaporation during a hydrological year (1850 mm/yr) and wasfound to be 6.8 mm/d. Overall, there is a good match between thesynthetic evaporation series, the daily measured values, and a 7day running average (Fig. 5b).
Finally, a time series of rainfall and potential (pan) evaporationfor 25 years (October 1990 to May 2015) was reconstructed byrepeating the existing records; this served as input to the 1Dmodel.
4.1.3. Natural vegetation parameterization
The natural vegetation of A. monosperma was represented in the
numerical model by a set of simplifications. First, its soil cover fraction (SCF) was assumed to be 26% for semi-stabilized dunes and
50% for stabilized landscape, that is, an average of reported values(Bar Kutiel et al., 2016). Its leaf area index (LAI) index was accordingly calculated using Eq. (3) (Ritchie, 1972):

1
LAI ¼   lnð1  SCFÞa

ð3Þwhere a is an empirical parameter, representing the radiation
extinction by the canopy, which equals a default value of 0.463.
An LAI of 0.65 and 1.49 were calculated for semi-stabilized and stabilized landscapes, respectively.
Second, the vertical root distribution was simplified as follows.
For both vegetation settings, no roots were set in the 0–5 cm depth
interval, and linear growth was assumed in the 5–30 cm interval.
For the semi-stabilized scheme, a logarithmic decrease in root density toward 0 was set in the 30–100 cm depth interval. For the stabilized landscape, vegetation was assumed to be moreacclimatized, with a maximal density in the 30–50 cm and a logarithmic decrease toward 0 in the 50–250 cm interval (Fig. 6).
Finally, for the root water uptake, we used the Feddes model(Feddes et al., 1978) with the default HYDRUS values for wheat(Wesseling et al., 1991). The latter was used in light of a lack ofmore appropriate values in the literature. Importantly, the wilting
point pressure head was P3 = 16,000 cm and the anaerobiosis
point pressure head was P0 = 0, i.e. extraction of water from thesoil can occur, virtually, under any saturation level.
To verify that the pre-set rooting depth and parameters of thewater stress response function in the Feddes model do not havea substantial effect on the results, we ran additional simulations
for the semi-stabilized scenario, with a higher (less negative) and346

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

Fig. 5. Daily precipitation (a) and daily evaporation (b) at the En Habesor meteorological station for years 1990/91–1998/99.

Fig. 6. Root distribution (a) and initial water content (b) in the upper 3 m of the soil for three scenarios. (For interpretation of the references to color in this figure legend, thereader is referred to the web version of this article.)

347

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

Fig. 7. Changes in cumulative surface flux, root uptake, recharge and evaporation for the amended semi-vegetated simulations during the 1999/2000–2008/9 period. Index
lists the changes made, see text for further explanation.

Table 5
Simulated water balance components for a 25-year period (mm and % of precipitation).
Bare soil

Semi-stabilized

Stabilized landscapePrecipitation (P)5328(100%)

5328(100%)

5328(100%)

Actual Evaporation (AE)

3624
(68%)
4
(0.1%)
1700(31.9%)
0(0%)+14(0.3%)
1686(31.6%)
+162
(3%)
1544
(29%)
20
(0.3%)

2498
(47%)
5
(0.1%)
2825(52.9%)
2170(40.7%)+5(0%)
650(12.2%)+6(0.1%)
683
(13%)
39
(0.7%)

2219
(42%)
4
(0.1%)
3105(57.9%)
2961(55.6%)+6(0%)
138(2.6%)+24(0.4%)
63(1%)+51(0.9%)

Runoff (RO)Net Infiltration (I)(I = P  AE  RO)
Transpiration (T)Change in storage in the root zone (M1)
Deep drainage (DD)(DD = P  AE  RO  T  M1)
Change in storage below the root zone (M2)
Recharge (R)(R = DD  M2 + ER)
Simulated mass balance error (an error of the numerical solution) (ER)a lower wilting point pressure head, P3 = 8000 cm andP3 = 32,000 cm, respectively), a lower anaerobiosis point pressure head (P0 = 30 cm), an extended range at which roots can
extract water at the maximum rates (by doubling P2H = 1000 cm
and P2L = 1800 cm), shallower and deeper root zones (70 and
130 cm, respectively), and a linear decrease in root density in the
30–100 cm interval. All in all, these runs demonstrate that theresults are virtually insensitive to the selected parameters(Fig. 7). For example, recharge rates (P/R) during the 1999/2000 –2008/9 period were 12.2% for all the simulations testing Feddes
parameters, and slightly lower (11.5–12.0%) for the three simulations with various rooting schemes.
4.1.4. Analysis of simulation results
The water balance was calculated for the root zone and the sediments below the root zone. Deep drainage (flux at the bottom of
the root zone), which is not included among the reported HYDRUS
results, was obtained by subtracting actual transpiration and thechange of storage in the root zone from net infiltration. Bottom flux
was serving as a proxy to recharge. Cumulative precipitation, evap-

otranspiration, net infiltration, soil moisture changes, and recharge
for the 25-year simulation were divided into annual values. The
annual values were tested for their probability distribution using
the quantile method (Helsel and Hirsch, 2002, p. 27). The temporalvariability of annual precipitation, annual deep drainage, andannual recharge were statistically analyzed. Then, annual values
were cross-correlated to demarcate or refute possible interrelations.
To detect the wetting front propagation in the sediment profile,simulated vertical fluxes as recorded at the bottom of each soil
horizon were used. The wetting front was assumed to propagate
downward following seasonal rains, with a decreasing amplitude
and an increasing time-lag (Wu et al., 1996). The merging of a
few annual peaks was expected to occur at larger depths, especiallyin scenarios with low annual deep drainage. Fortunately, the data
record contained two rainy years (1994/1995 and 2004/5) with
precipitation of +50% above average. These two years produced
fairly high net infiltration and deep drainage. The resulted flux
peaks could be traced down to the depth of the simulated profileand served to delineate the precipitation–recharge time lag. Next,

348

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

Fig. 8. Simulated annual precipitation (a), deep drainage (b), and recharge (c) for the three scenarios. (For interpretation of the references to color in this figure legend, thereader is referred to the web version of this article.)a few other peaks and lows in the flux time-series were correlated
to rain events with the same time lags. For the stabilized landscapescenario this was conducted only for the first peak, due to a highdetected time lag.
4.2. Results
The following section discusses water balances and different
fluxes for the three evaluated scenarios. For consistency, all lengthvariables reported in this section are in mm.
4.2.1. Water balance
Total recharge and five other components of the water balance,
as simulated for the 25-year period and the three scenarios, aresummarized in Table 5. Net infiltration (I) is the result ofsubtracting actual evaporation (AE) and runoff (RO) from precipitation (P) at the same time interval. Deep drainage (D), which is theflux at the bottom of the root zone, is obtained by subtracting theactual transpiration (T) and the change of storage in the root zone(M1) from net infiltration. And finally, recharge (R) is the result offurther subtracting changes in soil storage below the root zone(M2) from deep drainage. One should note that simulated massbalance errors are kept very low (less than 1% of the overall water
budget).
Annual time series of precipitation, deep drainage, and recharge
for the three scenarios are shown in Fig. 8. Simulated annual precipitation perfectly matched observed precipitation. Statisticalanalysis shows, with high confidence, that annual precipitation
and recharge are normally distributed (i.e., the null hypothesis
could not be rejected). Run-off was rather constant and small inE. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

349

Fig. 9. Simulated surface infiltration (a), and vertical fluxes (a and b), at different depths for the bare soil scenario. The arrows indicate the flux peaks that are discussed in the
text. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 6
Simulated recharge for a 25-year period using mean estimated hydraulic conductivities and their ±10% variations (total mm, and % of precipitation).

Ksl + 10%
Ks1
Ksl  10%

Ks2 + 10%

Ks2

Ks2  10%1574 mm, 29.5%
1551 mm, 29.1%
1533 mm, 28.6%

1567 mm, 29.4%
1544 mm, 29.0%
1517 mm, 28.5%

1559 mm, 29.3%
1533 mm, 28.8%
1508 mm, 28.3%

all scenarios, comparable to measurements and run-off estimates
for loessial sediments (Yair, 1990; Givati and Atzmon, 2009;
Eshtawi et al., 2015). The soil water content showed pronounced
annual fluctuations, which were superimposed on smooth, multiannual cycles. Soil water content, deep drainage, and recharge patterns are discussed in the following sections.

4.2.2. Bare soil
The initial water storage of 2076 mm (9.5% of pore volume) in
the 22-m sediment profile of the bare soil scenario was the highestof all three scenarios. During winter periods, the upper part of thesoil profile quickly gained moisture (up to 180 mm in some years),
while during the following months it gradually dried. Thesechanges can also be seen in the plots for the vertical fluxes (Fig. 9).

In the upper sandy layer (less than 4.3 m deep), large fluctuations in vertical fluxes can be observed closely resembling precipitation/infiltration events (Fig. 9a). Almost every annual cycle iscomprised of an instantaneous rise in the flux followed by its gradual decrease. The flux usually starts to increase 40–60 days after
the first major infiltration event. In some years, when infiltrationis small, the rise in the flux is absent (e.g., 1993/94, 1998/99,
2010/11) or relatively small (e.g., 1995/96, 1999/2000, 2007/8,
2008/9). At the bottom of the sandy layer (4.3 m deep), due to surface evaporation and redistribution, the range of simulated fluxes
is two orders of magnitude smaller than the range of surface infiltration fluxes.
Deeper in the sediment profile, the flux (and water content)
fluctuations are much attenuated and tend to smooth out individual annual cycles. Individual flux peaks and dips at greater depths
appear at larger time lags (Fig. 9b). For example, the peak that
occurred at a depth of 4.3 m in January 1995 reached the depthof 10.2 m in Nov 1995, the depth of 12.1 m in Feb 1996, the depthof 12.4 m in July 1996, and the depth of 21.7 m in June 1997
(downward arrows in Fig. 9b). The vertical flux was correspondingly attenuated from 1.26 mm/d to 0.29 mm/d, 0.28 mm/d,
0.26 mm/d, and finally 0.25 mm/d at increasing depths, respectively. That reflects a time lag of about one year (November
1994–November 1995) at the bottom of the first clay horizon

350

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

Fig. 10. Simulated surface infiltration (a), and vertical fluxes (a and b), at different depths for the semi-stabilized dunes scenario. The arrows indicate the flux peaks that are
discussed in the text. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 11. Deep drainage vs. precipitation for the semi-stabilized and stabilized scenarios.

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

and up to a time lag of 2.5 years at the bottom of the 22-m sediment column. These trends (attenuations of annual signals andincreasing time lags at deeper layers) are consistent with both theoretical considerations and field observations (e.g., Turkeltaub
et al., 2015).
The overall recharge over the 25-year period is equivalent to
29% of the total precipitation. This value is very close to the rangeof estimated recharge rates for sandy loess soils (31–40%, Seiler
and Gat, 2007). One should note that this is substantially less thanrecharge from sand dunes, which was lately reported to be about
72% of annual precipitation (Turkeltaub et al., 2015) and substantially more than an estimated range of recharge for silty loess soils(8–10%, Seiler and Gat, 2007).
Recharge rates were found to be relatively insensitive to smallchanges in hydraulic conductivities of the two soils. In further simulations, a ±10% change in soil hydraulic conductivities (Ks1 andKs2) resulted in a change of only ±2% in recharge (Table 6).

4.2.3. Semi-stabilized dunes
In the semi-stabilized dunes scenario, the initial water storage
in the sediment profile, as well as the average water storage during
the 25-year simulation, was 1901 mm. This value is about 8% lower

351

than that for the bare soil scenario and reflects the effects of plantroots in removing water from the soil profile.
The fluctuations in the flux (and water content) resemble thepattern described for the bare soil scenario. However, the magnitude of fluxes is smaller, and their downward propagation isslower. For example, the peak of the winter 1995 flux is
0.57 mm/d for the semi-stabilized dunes scenario (Fig. 10a), while
it is 1.25 mm/d for the bare soil scenario. The peak reached thedepth of 10.2 m in April 1996, the depth of 12.1 m in November
1996, the depth of 14.2 m in May 1997, and the depth of 21.7 m
in August 1999 (downward arrows in Fig. 10b). That corresponds
to a time lag of about 1.5 years at the bottom of the first clay horizon and almost five years at the bottom of the 22-m sedimentcolumn.
It can be concluded that consideration of root water uptake inthe upper 1 m of the soil profile reduced the overall soil water storage by about 8%, the peak flux by about 50%, and slowed the progress of the wetting front by about 50%. Additionally, the presenceof plants reduced the overall recharge to 13% of total precipitation,with an average recharge of 28 mm/yr. In fact, during eight years
there effectively was no deep drainage down from the root zone.
During the other four years, only minimal (<20 mm/y) deep drainage occurred. As suggested in Fig. 11, a linear relationship may

Fig. 12. Simulated surface infiltration (a), and vertical fluxes (a and b), at different depths for the stabilized landscape scenario. The arrows indicate the flux peaks that are
discussed in the text. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)352

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

be interpreted between annual deep drainage and annual precipitation when the latter exceeds 210 mm/yr. In other words, forsemi-stabilized loessial dunes, there is a threshold of 210 mm/yr
before deep drainage occurs.
4.2.4. Stabilized landscape
For the stabilized landscape, the initial water storage in the sediment profile was 1544 mm, and the 25-year average was1578 mm; both values are about 25% lower than for the bare soil
scenario.
The effects of the vegetation (roots) on the amplitude and temporal fluctuations in the flux (and water contents) discussed above
increased for this scenario, with roots extending down to a depthof 3 m. The infiltration peak of winter 1994/1995 reached thedepth of 4.3 m in April 1995 (6 months after the first major infiltration event) with a peak of 0.14 mm/d (Fig. 12a). It reached thedepth of 10.2 m in April 2001, with a flux of 0.0097 mm/d, andthe depth of 21.7 m at the end of the simulation, with a flux of
0.008 mm/d (downward arrows in Fig. 12b). That reflects a time
lag of about 6.5 years at the bottom of the first clay horizon andmore than 20 years at the bottom of the 22-m sediment column.
It can be concluded that consideration of plant roots in theupper 3 m of the soil profile reduced the overall soil water storage
by about 26%, the water flux peaks by 88–96%, and slowed the progress of the wetting fronts by a factor of 5–7. Overall recharge was
reduced in this scenario to only about 1% of total precipitation,with an average annual recharge of 2.5 mm/yr. In fact, during mostyears, there effectively was no deep drainage down from the rootzone. During only four years was the deep drainage substantial(>20 mm). When annual deep drainage is plotted against annual
precipitations, it appears that deep drainage occurs only when precipitation is higher than 270 mm/y, yet its magnitude is unpredictable (Fig. 11).

Fig. 13. Annual recharge vs. annual precipitation (a) of the same year and (b)shifted according to the wetting front propagation, i.e., +2 yr for bare soil, +5 yr forsemi-stabilized dunes, and +20 yr for the stabilized landscape. Linear trend lines
and R2 are shown for each scenario.

4.3. Discussion
The following section utilizes the simulated results to discuss
two prominent issues: (a) the relationship between annual precipitation and recharge, and (b) the effect of vegetation on the overall
recharge rates. These two factors have an important impact on ourability to predict recharge from precipitation and its representationin numerical models.
4.3.1. Recharge–precipitation relationshipAs noted before and demonstrated in Fig. 13a, there is no clear
correlation between annual recharge in loessial sediments andannual precipitation for any of the evaluated scenarios. Also, there
is no correlation between annual recharge and annual precipitationof the preceding year. These two correlations have been successfully tested to predict recharge for non-vegetated sandy soil profiles (Turkeltaub et al., 2015). There is no strong linear
correlation even between annual precipitation and annual
recharge shifted by 2.5, 5, and 20 years (Fig. 13b), reflecting the
propagation of the moisture front in the three evaluated scenarios,
respectively.
These results are not surprising and demonstrate the importance of the capacity of the soil profile in buffering short-term
fluxes (i.e., rain events). On average, only 7.0–9.5% of the pore volume is filled with water, allowing substantial fluctuations in thewater storage to occur. This buffering capacity is strengthened by
the layered structure of the sediment profile. In the studied case,vertical flow is accompanied by water movement from sandy clay
loam horizons with water contents of 13–15% to sandy loam horizons with water contents of 6–8% and back. That is, recharge, contrary to net infiltration, cannot be predicted on the basis of annual

Fig. 14. Overall simulated recharge (as %R/P) as a function of SCF and LAI.

precipitation, as the soil acts as a multi-annual storage to ‘buffer’
infiltration.
In agreement with previous studies (e.g., Wu et al., 1996), thisresearch demonstrates that recharge estimations, as a constantfraction of annual precipitation, may be an over-simplified
approach. Rather, dynamic models that account for the soil capacity to buffer water fluxes (e.g., Wu et al., 1997; Sheffer et al., 2010)
need to be used.
4.3.2. Effects of dune stabilization and natural vegetation on recharge
Similar to the recent literature (e.g., Hugenholtz and Koenig,
2014; Bar Kutiel et al., 2016), the results of this study suggest thatdune stabilization, through the establishment of local vegetation,
reduces the water content in the sediment profile and subsequently the groundwater recharge. Bare sandy loess may serve as
an important source of recharge and may provide long-term

353

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354
Table 7
Calculated travel time for matrix flow in the loessial sediments.

Bare soil
Semi-stabilized dunes
Stabilized landscapeMean recharge[mm/yr]Volume of water in thesoil profile [mm]Travel time to 22 m
depth [yr]Interpolated travel time to
30 m depth [yr]

63.5
28.0
2.6

2076
1901
1544

33
68
>500

45
92
808

recharge of up to 30% of precipitation. Semi-stabilized dunes with
scattered vegetation and a limited root depth of up to 1 m reduce
recharge rates by nearly half, whereas stabilized loess landscapeswith mature vegetation may prevent recharge altogether or allow
only minimal recharge of about 1% of precipitation. Simulation
results show that an increase in SCF and/or LAI results in a (nearlylinear) reduction of recharge as a fraction of precipitation (R/P)(Fig. 14). This trend has many physical and biological reasons,the main one being the presence of the rooting zone and thus
the greater ability to capture infiltrated water and return it to
the atmosphere in the form of transpiration (Kurtzman and
Scanlon, 2011). In other words, while recharge may not be accurately predicted on an annual basis, the long-term R/P can be determined solely from the vegetation coverage factors (LAI or SCF).

The HYDRUS simulations also show that vegetation maturity is
an important factor in determining the long-term recharge capacity of porous sediments. In general, recharge rates decrease with anincrease in the vegetation cover. For the evaluated layered sandy
loess, the recharge rate decreased from nearly 30% for bare soils
to only about 1% of precipitation under stabilized landscapes having a 50% vegetation cover.
These findings suggest that the long-term recharge (as R/P) canbe determined by accounting for the vegetation cover. Dynamic
models, which can account for the soil capacity to buffer water
fluxes (e.g., this study; Sheffer et al., 2010), should still be used
to estimate recharge time series.

4.3.3. Travel time
As a by-product of the recharge simulations, one can calculate
the travel time of water from the soil surface to the groundwatertable. A travel time of 33 years to a depth of 22 m can be calculated
for the bare soil scenario assuming a ‘piston flow’ through the soil
matrix, a recharge rate of 63.5 mm/yr, and water storage of2076 mm. Interpreting this value for a water table depth of 30 m
would result in a travel time of 45 years. Corresponding values
for the semi-stabilized dunes and stabilized landscape scenarios
would produce travel times of 92 years and more than 800 years,
respectively (Table 7). This time range is aligned with findings ofother arid zones loess recharge studies (Huang et al., 2013)

The authors thank Dr. Onn crouvi (Geological Survey of Israel)for valuable input regarding the loess origin, and Talia Hurowitz(IMS) for clarification of the meteorological data inventory. The
manuscript was improved following comments of three anonymous reviewers.

5. Summary
The transient nature of recharge through layered loessial sediments at the edge of arid zones was investigated using a Richards’
equation-based model. First, the sediment profile was conceptualized, and then the van Genuchten parameters of the soil hydraulicproperties were quantified using the dataset collected during thelarge infiltration test by Gvirtzman et al. (2008) and the inverse
method in the HYDRUS-2D numerical code. Simulation results
and the grain-size distribution classified the two soils, which
formed different horizons of the sediment profile, as sandy loam
and sandy clay loam. The calibrated saturated hydraulic conductivities of these two soils were found to be lower than previously
reported, with an average hydraulic conductivity of 80.2 cm/d forthe sandy loam and 17.3 cm/d for the sandy clay loam.
Next, diffuse recharge originating from natural precipitation
over layered loessial sediments was investigated using HYDRUS1D for three geomorphological conditions of bare soils, semivegetated dunes, and a stabilized landscape. The buffering effectof the layered sediment profile attenuated the vertical fluxes andmasked the signal of individual years, resulting in a multi-annual
recharge pattern. Thus, it is no surprise that no correlation was
found between annual recharge and annual precipitation. On thecontrary, groundwater recharge occurs year-round with relatively
small multi-annual changes. This research demonstrates, in agreement with related literature, that recharge estimations as a constant fraction of annual precipitation may be an over-simplified
approach.

Acknowledgements

References
Aish, A.M., 2014. Estimation of water balance components in the Gaza Strip with
GIS based WetSpass model. Civ. Environ. Res. 6 (11), 77–84.
Amir, N., Kafri, U., Herut, B., Shalev, E., 2013. Numerical simulation of submarine
groundwater flow in the coastal aquifer at the Palmahim area, theMediterranean coast of Israel. Water Resour. Manage. 27 (11), 4005–4020.
http://dx.doi.org/10.1007/s11269-013-0392-2.
Bar Kutiel, P., Katz, O., Ziso-Cohen, V., Divinsky, I., Katra, I., 2016. Water availability
in sand dunes and its implications for the distribution of Artemisia
monosperma.
Catena137,
144–151.
http://dx.doi.org/10.1016/
j.catena.2015.09.007.
Crouvi, O., Amit, R., Enzel, Y., Gillespie, A.R., 2010. Active sand seas and theformation of desert loess. Quatern. Sci. Rev. 29 (17–18), 2087–2098. http://dx.
doi.org/10.1016/j.quascirev.2010.04.026.
Eshtawi, T., Evers, M., Tischbein, B., 2015. Quantifying the impact of urban area
expansion on groundwater recharge and surface runoff. Hydrol. Sci. J. http://dx.
doi.org/10.1080/02626667.2014.1000916.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of Field Water Use and Crop
Yield. John Wiley & Sons, New York.
Fetter, C.W., 2001. Applied Hydrogeology, fourth ed. Prentice-Hall, New Jersey.
Givati, A., Atzmon, B., 2009. Model to predict flood-flow volume in the catchments:
Yarqon, Soreq, Lakhish, and Besor. Rep. Hydro/2009/01. Israeli Hydrological
Service (IHS), Jerusalem, 28 pp (Hebrew).
Gvirtzman, H., Shalev, E., Dahan, O., Hatzor, Y.H., 2008. Large-scale infiltration
experiments into unsaturated stratified loess sediments: monitoring andmodeling. J. Hydrol. 349, 214–229.
Hatzor, Y.H., Gvirtzman, H., Wainshtein, I., Orian, I., 2009. Induced liquefaction
experiment in relatively dense, clay-rich sand deposits. J. Geophys. Res. 114,
B02311. http://dx.doi.org/10.1029/2008JB005943.
Helsel, D.R., Hirsch, R.M., 2002. Statistical methods in water resources. In:
Techniques of Water-Resources Investigations of the United States Geological
Survey – Book 4, Hydrologic Analysis and Interpretation. USGS.
Hillel, D., 2004. Introduction to Environmental Soil Physics. Elsevier Academic Press,
Amsterdam.
Hopmans, J.W., Šimůnek, J., Romano, N., Durner, W., 2002. Inverse modeling oftransient water flow. In: Dane, J.H., Topp, G.C. (Eds.), Methods of Soil Analysis –Part 4, Physical Methods, third ed. SSSA, Madison, WI, pp. 963–1008.
Huang, T., Pang, Z., Edmunds, W.M., 2013. Soil profile evolution following land-use
change: implications for groundwater quantity and quality. Hydrol. Process. 27(8), 1238–1252. http://dx.doi.org/10.1002/hyp.9302.
Hugenholtz, C.H., Koenig, D.K., 2014. Sand dune stabilization reduces infiltrationand soil moisture: a case study from the northern Great Plains. Ecohydrology 7(4), 1135–1146. http://dx.doi.org/10.1002/eco.1445.
Jacques, D., Šimůnek, J., Timmerman, A., Feyen, J., 2002. Calibration of Richards’ andconvective–dispersion equations to field-scale water flow and solute transport

354

E. Dafny, J. Šimůnek / Journal of Hydrology 538 (2016) 339–354

under rainfall conditions. J. Hydrol. 259, 15–31. http://dx.doi.org/10.1016/
S0022-1694(01)00591-1.
Kurtzman, D., Scanlon, B.R., 2011. Groundwater recharge through vertisols:
irrigated cropland vs. natural land, Israel. Vadose Zone J. 10, 662–674.
Leterme, B., Mallants, D., Jacques, D., 2012. Sensitivity of groundwater rechargeusing climatic analogues and HYDRUS-1D. Hydrol. Earth Syst. Sci. 16 (8), 2485–
2497. http://dx.doi.org/10.5194/hess-16-2485-2012.
Lin, R., Wei, K., 2006. Tritium profiles of pore water in the Chinese loess unsaturated
zone: implications for estimation of groundwater recharge. J. Hydrol.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear
parameters. J. Soc. Ind. Appl. Math. 11 (1), 431–441. http://dx.doi.org/10.1137/
0111030.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity ofunsaturated porous media. Water Resour. Res. 12 (3), 513–522.
Petersen, C.T., Trautner, A., Hansen, S., 2008. Spatio-temporal variation ofanisotropy of saturated hydraulic conductivity in a tilled sandy loam soil. Soil
Tillage Res. 100 (1–2), 108–113. http://dx.doi.org/10.1016/j.still.2008.05.004.
Radcliffe, D.E., Šimůnek, J., 2010. Soil Physics with HYDRUS – Modeling andApplications. CRC Press, 373 pp.
Radcliffe, J.C., 2002. Pesticide Use in Australia. Australian Academy of Technological
Sciences and Engineering, Parkville, Victoria, 310 pp.
Ritchie, J.T., 1972. Model for predicting evaporation from a row crop withincomplete cover. Water Resour. Res. 8 (5), 1204–1213.
Ritter, A., Hupet, F., Muñoz-Carpena, R., Lambot, S., Van-clooster, M., 2003. Using
inverse methods for estimating soil hydraulic properties from field data as analternative to direct methods. Agric. Water Manage. 59, 77–96.
Scanlon, B., Healy, R., Cook, P., 2002a. Choosing appropriate techniques forquantifying groundwater recharge. Hydrogeol. J. 10 (1), 18–39. http://dx.doi.
org/10.1007/s10040-001-0176-2.
Scanlon, B.R., Christman, M., Reedy, R.C., Porro, I., Šimůnek, J., Flerchinger, G.F.,
2002b. Intercode comparisons for simulating water balance of surficial
sediments in semiarid regions. Water Resour. Res. 38 (12), 1323. http://dx.
doi.org/10.1029/2001WR001233.
Scanlon, B.R., Keese, K.E., Flint, A.L., Flint, L.E., Gaye, C.B., Edmunds, W.M., Simmers,I., 2006. Global synthesis of groundwater recharge in semiarid and arid regions.
Hydrol. Process. 20, 3335–3370.
Seiler, K.P., Gat, J.R., 2007. Groundwater Recharge from Run-off, Infiltration andPercolation. Springer.
Sheffer, N., Dafny, E., Gvirtzman, H., Navon, S., Frumkin, A., Morin, E., 2010.
Hydrometeorological daily recharge assessment model (DREAM) for theWestern Mountain Aquifer, Israel: model application and effects of temporal
patterns. Water Resour. Res. 46, W05510. http://dx.doi.org/10.1029/
2008WR007607.

Šimůnek, J., van Genuchten, M.T., Šejna, M., 2008. Development and applications of
the HYDRUS and STANMOD software packages and related codes. Vadose Zone
J. 7 (2), 587–600. http://dx.doi.org/10.2136/VZJ2007.0077.
Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M.T., 2009. The HYDRUS1D Software Package for Simulating the One-dimensional Movement of Water,
Heat, and Multiple Solutes in Variably-saturated Media Version 4.08. Dep. ofEnviron. Sci., Univ. of California, Riverside.
Tulmatz, Y., 1977. Hydrogeological Atlas of Israel – The Coastal Basin. The
Hydrological Service, Israel Water Authority, Jerusalem.
Turkeltaub, T., Kurtzman, D., Bel, G., Dahan, O., 2015. Examination of groundwater
recharge with a calibrated/validated flow model of the deep vadose zone. J.
Hydrol. 522, 618–627. http://dx.doi.org/10.1016/j.jhydrol.2015.01.026.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898.
Vrugt, J.A., Stauffer, P.H., Wöhling, T., Robinson, B.A., Vesselinov, V.V., 2008. Inverse
modeling of subsurface flow and transport properties: a review with newdevelopments. Vadose Zone J. 7 (2), 843–864.
Wang, W., Neuman, S.P., Yao, T., Wierenga, P.J., 2003. Simulation of large-scale field
infiltration experiments using a hierarchy of models based on public, genericand site data. Vadose Zone J. 2, 297–312.
Weinthal, E., Vengosh, A., Marei, A., Gutierrez, A., Kloppmann, W., 2005. The watercrisis in the Gaza Strip: prospects for resolution. Ground Water 43 (5), 653–660.
Wesseling, J.G., Elbers, J.A., Kabat, P., van den Broek, B.J., 1991. SWATRE: Instructions
for Input (Internal Notes). Winand Staring Centre, Wageningen, The
Netherlands.
Wierenga, P.J., Hills, R.G., Hudson, D.B., 1991. The Las Cruces trench site:characterization, experimental results, and one-dimensional flow predictions.
Water Resour. Res. 27 (10), 2695–2705.
Wu, J., Zhang, R., Yang, J., 1996. Analysis of rainfall-recharge relationships. J. Hydrol.
177 (1–2), 143–160. http://dx.doi.org/10.1016/0022-1694(95)02935-4.
Wu, J., Zhang, R., Yang, J., 1997. Estimating infiltration recharge using a response
function model. J. Hydrol. 198 (1–4), 124–139. http://dx.doi.org/10.1016/
S0022-1694(96)03309-4.
Yaalon, D.H., Dan, J., 1974. Accumulation and distribution of loess-derived deposits
in the semi-desert and desert fringe areas of Israel. Zeitschrift fur
Geomorphologie Supplementband 20, 91–105.
Yair, A., 1990. The role of topography and surface cover upon soil formation along
hillslopes in arid climates. Geomorphology 3, 287–299.
Zou, Z.Y., Young, M.H., Li, Z., Wierenga, P.J., 2001. Estimation of depth averaged
unsaturated soil hydraulic properties from infiltration experiments. J. Hydrol.
242 (1–2), 26–42. http://dx.doi.org/10.1016/S0022-1694(00)00385-1.