Българска Академия на науките . Bulgarian Academy of SciencesИНЖЕНЕРНА ГЕОЛОГИЯ И ХИДРОГЕОЛОГИЯ, 29
ENGINEERING GEOLOGY AND HYDROGEOLOGY, 29София, 2015. Sofia

Estimating groundwater recharge using
HYDRUS-1D
Jirka Šimůnek1
Department of Environmental Sciences, University of California Riverside, CA 92521, USA
E-mail: jiri.simunek@ucr.edu

1-

Ю. Шимунек. 2015. Оценка на подхранването на подземните води посредством компютърен
код HYDRUS-1D. Инженерна геология и хидрогеология, 29, стр. 25 - 36.
Резюме. Моделирането на водния баланс е един от методите за количествена оценка наподхранването на подземните води. Това е от изключителна важност за определяне наводните ресурси и уязвимостта от замърсяване на водоносните хоризонти. Освен това,определянето на водния баланс близо до повърхността е от ключово значение при оценкатана въздействията на изменението на климата и промяната на земеползването върху потокана подземните води.
Целта на статията е да представи математическото описание на процесите участващи вподхранването и включени в числовия код HYDRUS-1D като се представя и пример на базатана експериментални данни от района на гр. Ломел, Белгия. В този пример се разглежда водния
поток през еднометров почвен профил при два варианта – без растителност и с тревно покритие.
Представените резултати демонстрират, че подхранването на подземните води е по-голямопри гола почва. Тревната покривка “връща” повече вода обратно в атмосферата посредством
транспирация по сравнение с физическото изпарение от почва без растителност.
Ключови думи: поток в ненаситена среда, подхранване на подземните води,моделиране, HYDRUS-1D
Abstract. Water balance modeling is one of the few available tools that can be used to predictgroundwater recharge, which is critical for assessing water resources and aquifer vulnerability tocontamination. Quantification of the near surface water balance is extremely important for evaluating land-atmosphere interactions, and the impact of climate change and land-use change on thesubsurface flow.
The objective of this paper is to provide a mathematical description of processes implemented into the numerical model HYDRUS-1D and to provide a simple example of its application
to experimental data from Lommel, Belgium. In this example, we consider transient water flow
through a one meter deep soil profile, which is assumed either bare or covered with grass. Using thisexample we demonstrate that the soil profile with plants returns more water back to the atmosphere
via transpiration than the bare soil profile via evaporation, resulting in more groundwater recharge
in a bare soil profile.
Keywords: water flow, unsaturated zone, groundwater recharge, modeling, HYDRUS-1D
Geological Institute, Bulgarian Academy of Sciences © 25

Introduction
Water balance modeling is one of the few available tools that can be used to predictgroundwater recharge, which is critical for assessing water resources and aquifer
vulnerability to contamination (Scanlon et al., 2002). Quantification of the near surface
water balance is extremely important for evaluating land-atmosphere interactions,and the impact of climate change and land-use change on subsurface flow. There aredifferent types of models that can be used to evaluate the near-surface water balance.
These models vary in complexity from simple, analytical codes and storage-routing
or bucket-type codes to more complex codes that use the Richards’ equation.
To evaluate the near-surface water balance and estimate groundwater recharge
for specific soil, vegetation, and climate conditions is one of the very common
applications of the HYDRUS-1D model (Šimůnek et al., 2008) (see for example,
Adomako et al., 2010; Assefa et al., 2013; Awan et al., 2013; Dickinson et al., 2013;
Ficklin et al., 2010; Gates et al., 2011; Huo et al., 2014; Ibrahim et al., 2014; JimenezMartinez et al., 2009; Kurtzman et al., 2011; Le Coz et al., 2013; Leterme et al., 2012;
Lu et al., 2011; Mastrocicco et al., 2010; Scanlon et al., 2002; Scott et al., 2000;
Turkeltaub et al., 2014; Wang et al., 2009).
The objective of this brief note is to provide a mathematical description ofprocesses implemented into the numerical model HYDRUS-1D and to provide a
simple example of its application to experimental data from Belgium.
Model Description
Governing Flow Equation
One-dimensional uniform water movement in a partially saturated rigid porous
medium is described in HYDRUS-1D (Šimůnek et al., 2008) by a modified form of
the Richards’ equation using the assumptions that the air phase plays an insignificant
role in the liquid flow process and that water flow due to thermal gradients can beneglected:∂θ ∂ ⎡ ⎛ ∂h ⎞ ⎤=K ⎜ + 1⎟ - S∂t ∂z ⎣⎢ ⎝ ∂z ⎠ ⎦⎥

				

(1)

where h is the water pressure head [L], θ is the volumetric water content [L3 L-3], t is
time [T], z is the spatial coordinate [L] (positive upward), S is the sink term [L3 L-3 T-1],
and K is the unsaturated hydraulic conductivity function [L T-1] given by

K (h, x) = K s ( z ) K r (h, z )

			

(2)

where Kr is the relative hydraulic conductivity [-] and Ks the saturated hydraulicconductivity [L T-1].
Root Water Uptake
The sink term, S, is defined as the volume of water removed from a unit volume ofsoil per unit time due to plant water uptake. Feddes et al. (1978) defined S as

S ( h) = α ( h) S p

				

(3)

where the root-water uptake water stress response function α(h) is a prescribed
dimensionless function of the soil water pressure head (0 ≤ α ≤ 1), and Sp the potential
water uptake rate [T‑1]. Water uptake is assumed to be zero close to saturation (i.e.,wetter than some arbitrary “anaerobiosis point”, h1). For h<h4 (the wilting point
26 Geological Institute, Bulgarian Academy of Sciences ©

pressure head), water uptake is also assumed to be zero. Water uptake is considered
optimal between pressure heads h2 and h3, whereas for pressure head between h3 andh4 (or h1 and h2), water uptake decreases (or increases) linearly with h. The variable
Sp in (3) is equal to the water uptake rate during periods of no water stress whenα(h)=1.
When the potential water uptake rate is non‑uniformly distributed over theroot zone, Sp becomes

S p = b( z )Tp

					

(4)

where Tp is the potential transpiration rate [L T-1] and b(z) is a normalized water
uptake distribution [L-1]. This function describes the spatial variation of the potential
extraction term, Sp, over the root zone, and is obtained by normalizing any arbitrarily
measured or prescribed root distribution function.
The actual water uptake distribution is obtained by substituting (4) into (3):S (h, z ) = α (h, z ) b( z )Tp

			

(5)

whereas the actual transpiration rate, Ta, is obtained by integrating (5) over the
rooting depth, LR, as follows:

Ta = ∫ S (h, z ) dz = Tp ∫ a ( h, z ) b( z ) dz
LR

(6)

LR

Unsaturated Soil Hydraulic Properties
The unsaturated soil hydraulic properties, θ(h) and K(h), in (1) are in general highly
nonlinear functions of the pressure head. HYDRUS permits the use of five different
analytical models for the hydraulic properties. In this work we use the soil-hydraulicfunctions of van Genuchten (1980) who used the statistical pore-size distribution
model of Mualem (1976) to obtain a predictive equation for the unsaturated hydraulicconductivity function in terms of soil water retention parameters. The expressions of
van Genuchten (1980) are given byθs - θr⎧⎪ θr +[1+ |α h|n ]m
θ (h) = ⎨⎪θs⎩
K (h) = K s Sel [1 - (1 - Se1/ m ) m]

h<0

		

(7)

h≥0
2

		

		

(8)

where

m = 1-1 / n , n > 1

			

(9)

and the effective saturation Se is defined as follows:
Geological Institute, Bulgarian Academy of Sciences ©

27

Se =

θ - θrθs -θr

					

(10)

The above equations contain five independent parameters: θr, θs, α, n, and
Ks. The pore-connectivity parameter l in the hydraulic conductivity function was
estimated (Mualem, 1976) to be about 0.5 as an average for many soils.
Initial and Boundary Conditions
The solution of Eq. (1) requires knowledge of the initial distribution of the pressure
head within the flow domain:h( z , t ) = hi ( z )t = t0				(11)where hi [L] is a prescribed function of z, and to is the time when the simulation
begins.
The surface boundary is exposed to atmospheric conditions. The potential fluid
flux across this interface is controlled exclusively by external conditions. However,
the actual flux depends also on the prevailing (transient) soil moisture conditions near
the surface. The soil surface boundary condition may change from a prescribed flux
to a prescribed head type condition (and vice-versa). The numerical solution of (1)is obtained by limiting the absolute value of the surface flux by the following two
conditions:⎛ ∂h ⎞- K ⎜ + 1⎟ ≤ E⎝ ∂z ⎠

			

(12)

			

(13)

at z = LandhA ≤ h ≤ hS

at z = L

where E is the maximum potential rate of infiltration or evaporation under the current
atmospheric conditions [L T-1], and hA and hS are, respectively, minimum and maximum
pressure heads at the soil surface allowed under the prevailing soil conditions [L].
The value for hA is determined from the equilibrium conditions between soil waterand atmospheric water vapor, whereas hS is usually set equal to zero; if positive, hS
represents a small layer of water ponded which can form on top of the soil surface
during heavy rains before initiation of runoff. One option in HYDRUS is to assume
that any excess water on the soil surface above zero will be immediately removed.
When one of the end points of (13) is reached, a prescribed head boundary condition
will be used to calculate the actual surface flux.
Another option in HYDRUS is to permit water to build up on the surface. If
surface ponding is expected, a “surface reservoir” boundary condition of the typedh
⎛ ∂h ⎞- K ⎜ + 1⎟ = q0 (t ) dt
⎝ ∂z ⎠

at z = L		(14)

may be applied. The flux q0 in this equation is the net infiltration rate, i.e., thedifference between precipitation and evaporation. Equation (14) shows that the heighth (L, t) of the surface water layer increases due to precipitation, and reduces because
of infiltration and evaporation.
28 Geological Institute, Bulgarian Academy of Sciences ©

A free drainage (zero pressure head gradient) boundary condition is usuallyspecified at the bottom of the soil profile. This boundary condition assumes that the
groundwater table is deep below the considered soil profile and that the bottom flux
is driven only by gravity.

Evaluation of Potential Evapotranspiration
Hargreaves Formula
The potential evapotranspiration can be evaluated using the Hargreaves formula (e.g.,
Jensen et al., 1997):

ETp = 0.0023Ra (Tm + 17.8 ) ΔTm

			

(15)

where Ra is the extraterrestrial radiation in the same units as ETp [e.g., mm d-1 orJ m-2 s-1], Tm is the daily mean air temperature, computed as an average of maximum
and minimum air temperatures [oC], ∆Tm is the temperature range between mean
daily maximum and minimum air temperatures [oC]. The extraterrestrial radiation,
Ra [J m-2 s-1], can be calculated as follows:

Ra =

Gsc
d r (ωs sin ϕ sin δ + cos ϕ cos δ sin ωs )π

(16)where Gsc is the solar constant [J m-2 s-1] (1360 W m-2), ϕ is the site latitude [rad], ωs isthe sunset hour angle [rad], dr is the relative distance between Earth and Sun [-], and δ
is the solar declination [rad]. The last three variables are calculated as follows:ωs = arccos(- tan ϕ tan δ )

				

(17)

⎛ 2π ⎞
d r = 1 + 0.033 cos ⎜J⎟⎝ 365 ⎠

				

(18)

		

(19)

⎛ 2π
⎞
J -1.39 ⎟⎝ 365⎠

δ = 0.409 sin ⎜where J is the number of the day in the year [-].
The minimum pressure head at the soil surface allowed under the prevailingsoil conditions, hA [L], is calculated from the air humidity, Hr [-], as follows:⎡ h Mg ⎤H r = exp ⎢ − A⎣ RT ⎥⎦RTln( H r )hA = −
Mg

				

(20)

where M is the molecular weight of water [g mol-1] (=0.018015 kg mol-1), g is thegravitational acceleration [L T-2], (=9.81 m s-2), and R is the universal gas constant[J mol-1 K-1] (=8.314 J mol-1 K-1) [M L2 T-2 mol-1 K-1].
Geological Institute, Bulgarian Academy of Sciences ©

29

Penman-Monteith Combination Formula for Reference Evapotranspiration
Defining reference evapotranspiration (ETo) as the rate of evapotranspiration froma hypothetic crop with an assumed crop height of 12 cm, a fixed canopy resistanceof 70 s m-1 and an albedo of 0.23, closely resembling the evapotranspiration froman extensive surface of green grass of uniform height, actively growing, completely
shading the ground and not short of water, the estimation of the ETo can be determinedwith the combination formula based on the Penman-Monteith approach. When
combining the aerodynamic and radiation terms, the combination formula can benoted as (FAO, 1990):

ET0 =

900
U 2 (ea - ed )T + 273
Δ + γ (1 + 0.34U 2 )0.408 Δ ( Rn - G ) + γ(21)

where ET0 is the reference crop evapotranspiration [mm d-1], Rn is the net radiation at
crop surface [MJ m-2d-1], G is the soil heat flux [MJ m-2 d-1], T is the average temperature [oC], U2 is the wind speed measured at 2 m height [m s-1], (ea-ed) is the vapour
pressure deficit [kPa], Δ is the slope vapour pressure curve [kPa oC-1], γ is the psychrometric constant [kPa oC-1], and 900 is the conversion factor.
When no measured radiation data are available, the net radiation can be estimated as follows:

Rn = Rns - Rnl

					

n⎞
⎛
Rns = 0.77 ⎜ 0.25 + 0.50 ⎟ Ra
N⎠
⎝

			

n
⎛
⎞
Rnl = 2.45 .10−9 ⎜ 0.9 + 0.1⎟ 0.34 - 0.14 ed
N
⎝
⎠

(

G = 0.14(Tmonth , n - Tmonth , n −1 ) ≈ 0) (T + T )
4
kx

4
kn

			

(22)

(23)

(24)
(25)

where Rn is net radiation [MJ m-2 d-1], Rns is the net shortwave radiation [MJ m-2 d-1],
Rnl is the net longwave radiation [MJ m-2 d-1], Ra is the extraterrestrial radiation[MJ m-2 d-1], n/N is the relative sunshine fraction [-], Tkx and Tkn are the maximum andminimum temperatures [K], respectively; ed is the actual vapour pressure [kPa], andG is the soil heat flux [MJ m-2 d-1].
Partitioning of ETp into Ep and Tp
Potential evaporation and transpiration fluxes can also be calculated from potential
evapotranspiration using Beer’s law that partitions the solar radiation component of
the energy budget via interception by the canopy (Ritchie, 1972) as follows:

30 Geological Institute, Bulgarian Academy of Sciences ©

Tp = ETp (1 - e -k *laI )= ETp SCF(26)

E p = ETp e -k *laI = ETp (1 - SCF )where ETp, Tp, and Ep are potential evapotranspiration, transpiration and evaporation
fluxes [L T-1], respectively, LAI is the leaf area index [-], SCF is the soil cover
fraction [-], and k is a constant governing the radiation extinction by the canopy [-]as a function of sun angle, the distribution of plants, and the arrangement of leaves(between 0.5-0.75).

Demonstrative Example
To demonstrate the use of the HYDRUS-1D software package to evaluate thenear-surface mass balance and to estimate groundwater recharge, we have adopted
an example from the HYDRUS short course, which uses experimental data from
Seuntjens (2000). In this example we consider transient water flow through a 1-m deep
multi-layered Podzol soil profile. Transient flow is induced by atmospheric boundary
conditions. Either a fully plant-covered (grass) or bare soil profile is considered, thus
restricting the atmospheric boundary conditions to daily values of precipitation, andtranspiration or evaporation, respectively.
Input Data
Soil hydraulic and physical parameters (Table 1; Fig. 1) of the dry Spodosol locatedat the “Kattenbos” site near Lommel, Belgium were taken from Seuntjens (2000,
Tables 3.1 and 7.1).
The initial water content conditions (Fig. 2) are defined assuming a constantflux of 0.12 cm day-1 and a free-drainage lower boundary condition. The flux
corresponds to the long-term (1972-1981) actual infiltration rate (precipitation actual evapotranspiration). The upper boundary condition involves daily values ofprecipitation and potential evaporation (for a bare profile) or potential transpiration(for a planted profile) fluxes defined using meteorological data from the Brogel station
weather (Belgium) for 1972 (Fig. 3). In a simulation with plants, it is assumed thatroots are uniformly distributed in the top 30 cm.

Table 1. Soil hydraulic and physical properties of seven soil horizonsТаблица 1. Хидравлични свойства и физични параметри за седем почвени хоризонтаHorizon

depth

Bulk
density, ρorganic
carbonresidual
watercontent, θr

Saturated
water
content, θsempirical
parameter, αempirical
parameter, nSaturatedhydraulicconductivity, Ks
cm d-1cm

g cm-3

%

-

-

cm-1

-

A

0–7

1.31

2.75

0.065

0.48

0.016

1.94

95

e

7 – 19

1.59

0.75

0.035

0.42

0.015

3.21

311

Bh1

19 – 24

1.30

4.92

0.042

0.47

0.016

1.52

39

Bh2

24 – 28

1.38

3.77

0.044

0.46

0.028

2.01

864

Bc

28 – 50

1.41

0.89

0.039

0.46

0.023

2.99

1210

c1

50 – 75

1.52

0.12

0.030

0.42

0.021

2.99

1210

c2

75 – 100

1.56

0.08

0.021

0.39

0.021

2.99

1210

Geological Institute, Bulgarian Academy of Sciences ©

31

Fig. 1. Soil hydraulic properties (retention curves – left; hydraulic conductivity functions – right)of the seven soil horizons of the Podzol soil profile
Фиг. 1. Хидравлични свойства (криви на влагозадържащата способност на почвата θ = f(h)- вляво; зависимост на коефициента на влагопренасяне от всмукващия потенциал K = f(h) вдясно) на седемте почвени хоризонти на подзолистта почва

Fig. 2. Initial water contents
Фиг. 2. Начална обемна влажност

Fig. 3. Daily values of precipitation and potential evapotranspiration fluxes defined using meteorological data from the Brogel weather station (Belgium) for 1972Фиг. 3. Ежедневни стойности на валежите и потенциалната евапотранспирация
определени по метеорологични данни за станция Brogel (Белгия) за 1972 г.
32 Geological Institute, Bulgarian Academy of Sciences ©

Output Data
Figure 4 shows water contents and water fluxes at depths of 50 and 100 cm during the
one-year long simulation for a bare soil profile. Both water contents and water fluxes
reflect the surface drivers of water flow in the soil profile, i.e., variable precipitation.

Fig. 4. Water contents (left) and water fluxes (right) at depths of 50 and 100 cm
Фиг. 4. Обемна влажност (вляво) и водни потоци (вдясно) на дълбочини 50 и 100 cm(вдясно)Figure 5 shows main cumulative water fluxes for both bare and planted soil
profile. For the bare soil profile, Figure 5 shows potential and actual cumulative
surface fluxes, cumulative infiltration and evaporation, and cumulative recharge. The
potential surface flux (positive downwards in Fig. 5) is the sum of precipitation andpotential evaporation. Since the saturated hydraulic conductivity of the surface layer
is significantly higher than the highest values of precipitation and thus all precipitationwater infiltrates into the soil profile, any difference between the potential and actual
surface fluxes is the result of the reduction of evaporation from its potential to actual
values. Note that this reduction is about 10 cm (less water leaving the soil profile,and thus overall more water moving downwards through the soil surface). Also notethat since the soil profile has little retention capacity, the cumulative recharge follows
closely the actual surface flux. Any difference between these two fluxes is the change

Fig. 5. Main cumulative water fluxes for the bare soil profile (left) and the soil profile with grass(right)Фиг. 5. Основни кумулативни потоци за гол почвен профил (вляво) и затревен почвен
профил (вдясно)Geological Institute, Bulgarian Academy of Sciences ©

33

in water storage in the soil profile.
For the planted soil profile, Figure 5 shows potential and actual cumulativetranspiration, cumulative infiltration, and cumulative recharge. Infiltration is equal
to both potential and actual surface fluxes, since evaporation is considered to be
zero for cropped soil profile, and again, since the saturated hydraulic conductivity ofthe surface layer is significantly higher than the highest values of precipitation, allprecipitation water infiltrates into the soil profile. Note that the difference between
potential and actual transpiration (about 3 cm) (for a cropped profile) is significantly
smaller than the difference between potential and actual evaporation (about 10 cm)(for a bare soil profile). This indicates that the soil profile with plants returns morewater back to the atmosphere than the bare soil profile.
Figure 6 shows cumulative recharge for the bare and planted soil profiles,
reflecting the processes discussed above, i.e., that the soil profile with plants returns
more water back to the atmosphere via transpiration than the bare soil profile via
evaporation, resulting in more recharge in a bare soil profile.

Fig. 6. Cumulative recharge for the bare and planted soil profiles
Фиг. 6. Кумулативно подхранване на подземните води за гол почвен профил и за затревен
почвен профил (вдясно)Summary
In this brief note we have provided mathematical description of processes implemented
into the numerical model HYDRUS-1D and a simple example of its application toexperimental data from Belgium. In this simple example, we considered transientwater flow through a 1-m deep multi-layered Podzol soil profile, which was assumed
either bare or covered with grass. Using this example we demonstrated that the soil
profile with plants returns more water back to the atmosphere via transpiration thanthe bare soil profile via evaporation, resulting in more groundwater recharge in a baresoil profile.
References
Adomako, D., P. Maloszewski, C. Stumpp, S. Osae, T. T. Akiti. 2010. Estimating
groundwater recharge from water isotope (δ2H, δ18O) depth profiles in the
Densu River basin, Ghana. - Hydrological Sciences Journal, 55, 8, 14051416.
34 Geological Institute, Bulgarian Academy of Sciences ©

Assefa, K.A., A. D. 2013. Woodbur. Transient, spatially varied groundwater
recharge modeling. - Water Resources Research, 49, 4593–4606, doi:10.1002/
wrcr.20332.
Awan, U. K., B. Tischbein, C. Martius. 2013. Combining hydrological modeling and
GIS approaches to determine the spatial distribution of groundwater recharge
in an arid irrigation scheme. - Irrig. Sci., 31, 793–806, doi: 10.1007/s00271012-0362-0.
Dickinson, J. E., T. P. A. Ferré, M. Bakker, B. Crompton. 2014. A screening tool fordelineating subregions of steady recharge within groundwater models, Vadose
Zone Journal, 13, 6, 15 pp., doi: 10.2136/vzj2013.10.0184,
Food and Agriculture Organization of the United Nations. 1990. Expert consultation
on revision of FAO methodologies for crop water requirements, ANNEX V,
FAO Penman-Monteith Formula, Rome Italy.
Feddes, R. A., P. J. Kowalik, H. Zaradny. 1978. Simulation of Field Water Use andCrop Yield, John Wiley & Sons, New York, NY.
Ficklin, D. L., E. Luedeling, M. Zhang. 2010. Sensitivity of groundwater recharge
under irrigated agriculture to changes in climate, CO2 concentrations andcanopy structure. - Agricultural Water Management, 97, 7, 1039-1050.
Gates, J. B., B. R. Scanlon, X. Mu, L. Zhang. 2011. Impacts of soil conservation on
groundwater recharge in the semi-arid Loess Plateau, China. - Hydrogeology
Journal, 19, 865–875.
Huo, S., J. Menggui, X. Liang, D. Lin. 2014. Changes of vertical groundwater recharge
with increase in thickness of vadose zone simulated by one-dimensionalvariably saturated flow model, Journal of Earth Science, 25, 6, 1043-1050.
Ibrahim, M., G. Favreau , B. R. Scanlon, J. L. Seidel, M. Le Coz, J. Demarty, B.
Cappelaere. 2014. Long-term increase in diffuse groundwater recharge
following expansion of rainfed cultivation in the Sahel, West Africa. Hydrogeology Journal, 22, 6, 1293-1305.
Jensen, D. T., G. H. Hargreaves, B. Temesgen, R. G. Allen. 1997. Computation of Eto
under nonideal conditions. - J. of Irrig. Drainage, 123, 5, 394-400.
Jimenez-Martinez, J., T. H. Skaggs, M. Th. van Genuchten, L. Candela. 2009. A rootzone modelling approach to estimating groundwater recharge from irrigatedareas. - J. Hydrol., 367, 1-2, 138-149.
Kurtzman, D., B. R. Scanlon. 2011. Groundwater recharge through Vertisols: Irrigatedcropland vs. natural land, Israel. - Vadose Zone Journal, 10, 2, 662-674.
Le Coz, M., G. Favreau, S. D. Ousmane. 2013. Modeling increased groundwater
recharge due to change from rainfed to irrigated cropping in a semiarid region.
- Vadose Zone Journal, 12, 2, 12 pp., doi:10.2136/vzj2012.0148.
Leterme, B., D. Mallants, D. Jacques. 2012. Sensitivity of groundwater rechargeusing climatic analogues and HYDRUS-1D. - Hydrol. Earth Syst. Sci., 16,2485–2497, doi:10.5194/hess-16-2485-2012,
Lu, X., M. Jin, M. Th. van Genuchten, B. Wang. 2011. Ground water recharge at five
representative sites in the Hebei Plain of China: Case study. - Ground Water,
49, 2, 286-294.
Mastrocicco, M., N. Colombani, E. Salemi, G. Castaldelli. 2010. Numericalassessment of effective evapotranspiration from maize plots to estimate
groundwater recharge in lowlands. - Agricultural Water Management, 97, 9,1389-1398.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous media. - Water Resour. Res., 12, 3, 513-522.
Ritchie, J. T. 1972. Model for predicting evaporation from a row crop with incomplete
cover. - Water Resour. Res., 8, 5, 1204-1213.
Scanlon, B. R., M. Christman, R. Reedy, I. Porro, J. Šimůnek, G. Flerchinger.
2002. Intercode comparisons for simulating water balance of surficial
sediments in semiarid regions. - Water Resour. Res., 38, 12, 13-23,
doi:10.1029/2001WR001233, 59.1-59.16.
Scott, R. L., W. J. Shuttleworth, T. O. Keefer, A. W. Warrick. 2000. Modeling multiyear
observations of soil moisture recharge in the semiarid American Southwest. Geological Institute, Bulgarian Academy of Sciences ©

35

Water Resour. Rer., 36, 8, 2233-2247.
Seuntjens, P., Reactive solute transport in heterogeneous porous media. Cadmium
leaching in acid sandy soil. PhD, University of Antwerp, 236 p., 2000.
Šimůnek, J., M. Th. van Genuchten, M. Šejna. 2008. Development and applications of
the HYDRUS and STANMOD software packages and related codes. - Vadose
Zone Journal, 7, 2, 587-600, doi:10.2136/VZJ2007.0077.
Turkeltaub, T., O. Dahan D. Kurtzman. 2014. Investigation of groundwater recharge
under agricultural fields using transient deep vadose zone data. - Vadose Zone
Journal, 13, 4, 13 pp., doi:10.2136/vzj2013.10.0176.
van Genuchten, M. Th. 1980. A closed‑form equation for predicting the hydraulicconductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892‑898.
Wang, T., V. A. Zlotnik, J. Šimůnek, M. G. Schaap. 2009. Using pedotransfer functions
in vadose zone models for estimating groundwater recharge in semiarid regions.
- Water Resour. Res., 45, 12 pp., doi:10.1029/2008WR006903 W04412.

Постъпила: 02.11.2015				

Received: 02 November 2015Приета: 07.12.2015				

Accepted: 07 December 201536 Geological Institute, Bulgarian Academy of Sciences ©