Case Study/

Groundwater Recharge at Five Representative
Sites in the Hebei Plain, China
by Xiaohui Lu1,3 , Menggui Jin1 , Martinus Th. van Genuchten2 , and Bingguo Wang1

Abstract
Accurate estimates of groundwater recharge are essential for effective management of groundwater, especiallywhen supplies are limited such as in many arid and semiarid areas. In the Hebei Plain, China, water shortageis increasingly restricting socioeconomic development, especially for agriculture, which heavily relies ongroundwater. Human activities have greatly changed groundwater recharge there during the past several decades.
To obtain better estimates of recharge in the plain, five representative sites were selected to investigate the effectsof irrigation and water table depth on groundwater recharge. At each site, a one-dimensional unsaturated flow
model (Hydrus-1D) was calibrated using field data of climate, soil moisture, and groundwater levels. A sensitivity
analysis of evapotranspirative fluxes and various soil hydraulic parameters confirmed that fine-textured surface
soils generally generate less recharge. Model calculations showed that recharge on average is about 175 mm/yearin the piedmont plain to the west, and 133 mm/year in both the central alluvial and lacustrine plains and thecoastal plain to the east. Temporal and spatial variations in the recharge processes were significant in responseto rainfall and irrigation. Peak time-lags between infiltration (rainfall plus irrigation) and recharge were 18 to
35 days in the piedmont plain and 3 to 5 days in the central alluvial and lacustrine plains, but only 1 or 2 daysin the coastal plain. This implies that different time-lags corresponding to different water table depths must be
considered when estimating or modeling groundwater recharge.

Introduction
The Hebei plain is one of the most productive agricultural areas of northern China, occupying an area of
about 73,000 km2 and having a population of approximately 50 million. Groundwater is the main source ofwater, with agricultural usage about 70% of total water
1 Corresponding author: MOE Key Laboratory of BGEG, Instituteof Water Resources and School of Environmental Studies, China
University of Geosciences, Wuhan, Hubei 430074, China; 86-2767883461; fax: 86-27-87436235; mgjin@cug.edu.cn
2 Department of Mechanical Engineering, LTTC/COPPE, Federal
University of Rio de Janeiro, UFRJ, Rio de Janeiro, RJ 21945-970,
Brazil.
3 School of Earth Sciences and Engineering, Hohai University,
Nanjing, Jiangshu, 210098, China.
Received January 2009, accepted December 2009.
Copyright © 2010 The Author(s)Journal compilation © 2010 National Ground Water Association.
doi: 10.1111/j.1745-6584.2009.00667.x

NGWA.org

use. Excessive use of both surface and groundwater inthis densely populated semiarid plain has led to a rapid
decline of groundwater levels since the 1980s, as well as
aggravated environmental problems. The fresh water table
is declining at a rate of 1 to 2 m/year, with the total decline
during the past decades being 10 to 30 m in many areas.
Especially affected is the piedmont plain to the west where
the water table is currently at depths of up to 60 m or
more. Furthermore, water levels in fresh confined aquifers
have declined to 100 m below surface in many parts of theplain, even in areas where levels were originally near (oreven above) land surface. Water shortage and associatedenvironmental problems are increasingly threatening the
local socio-economic sustainability of the area (Jin et al.
1999; Yang et al. 2006).
The long-term natural groundwater cycle in the Hebei
plain has been greatly changed over the years due tohuman activities. Vertical groundwater recharge fromGROUND WATER

1

precipitation and irrigation return flow, and evapotranspiration (ET) are the main components of the water cyclein the plain, with the unsaturated zone playing an important role. Water-saving agricultural practices since the1990s, along with the declining water tables and increased
thickness of the vadose zone, have strongly affected therecharge processes and net groundwater recharge (Jinet al. 2000; Wang et al. 2008a). As in many other arid
and semiarid areas of the world with limited groundwater
resources, knowledge of how the groundwater rechargealters with changes in natural and anthropogenic factors
is critical for effective management of available watersupplies.
Groundwater recharge in the Hebei plain is generally
estimated using empirical equations, usually by simply
multiplying the precipitation rate by an empirical rechargecoefficient; these values are used in groundwater flow
models (Wang et al. 2008b). This empirical approach
does not consider explicitly the effects of local land use,
irrigation, and water table depth on recharge.
The main objective of this study is to provide a more
process-based investigation of the effects of irrigation andwater table depth on groundwater recharge in the Hebei
plain. The Hydrus-1D variably saturated flow model(Šimůnek et al. 2005) was used at five representative sites.
Hydrus-1D has been widely used to simulate a broad
range of water flow and solute transport processes in variably saturated media (Šimůnek et al. 2008, and references
therein), although applications to recharge processes havebeen relatively limited (Scanlon et al. 2003; JimenezMartinez et al. 2009). Additional objectives were used toevaluate temporal and spatial variability in the recharge
processes in the Hebei plain. The study also provided an
opportunity for comparing recharge rates estimated with
unsaturated flow modeling to those derived using chemical
tracer studies (e.g., Wang et al. 2008a) at similar sites.

Site and Data Description
The Hebei plain is a semiarid to semi-humid area witha monsoon-dominated climate. Mean annual precipitation
is 500 to 600 mm, of which more than 80% occurs in themonths of June through September, whereas mean annual
evaporation rates (of open water) are in the order of 1000
to 1300 mm (Chen et al. 2008). Geomorphologically,
the Hebei plain comprises three parts from west to east:

Figure 1. Locations of the five representative sites in theHebei plain, showing the location of the five sites (stars),rivers (in blue), boundaries of the study area (in brown),
and boundaries of geomorphological units (dashed line). The
roman numerals indicate: the piedmont plain (I), the alluvial
and lacustrine plains (II), and the coastal plain (III).

the piedmont plain, the alluvial and lacustrine plains,
and the coastal plain (Figure 1). The water table in the
piedmont plain is relatively deep, generally 10 to 40 m,
but up to 60 m at some locations, whereas the unsaturated
zone consists mainly of silt and fine sand (Table 1). The
water table depth in the alluvial and lacustrine plains is in
the order of 1 to 20 m, with the unsaturated zone mostly
consisting of silty clay and silt. By comparison, the watertable is relatively shallow in the coastal zone near theBohai Sea, generally between 1 and 4 m below the land
surface.
Five representative sites in the Hebei plain wereselected for the recharge studies (Table 1): the Luquan(LQ) and Luancheng (LC) sites in the piedmont plain
to the west, the Hengshui (HS) and Dezhou (DZ) sites

Table 1
General Features of the Representative Sites

2

Site

Geomorphology

Qijie in LQ
LC
HS
DZ
CZ

Piedmont plain
Piedmont plain
Alluvial and lacustrine plains
Alluvial and lacustrine plains
Coastal plain

X. Lu et al. GROUND WATER

Water Table Depth (m)
10–25
30–35
3–5
2–6
1–4

Lithology
Silt and fine sandSilty clay
Silty clay, clay, and silt
Silty clay
Silt and silty clay

NGWA.org

in the central alluvial and lacustrine plains, and theCangzhou (CZ) site in the coastal plain to the east (Table 1
and Figure 1). Soil water contents, pressure heads, andgroundwater levels were monitored every 5 to 10 days
from 2003 to 2005. Pressure heads were measured using
tensiometers. Water contents were measured using an
L520 neutron probe (Jiangsu Academy of Agricultural
Sciences, Nanjing, China) at the LC and HS sites,
and using MP160 moisture probes (ICT International
Pty Ltd., Armidale, Australia) at the other sites. Daily
meteorological data were obtained from local weather
stations, which were all located on site, except for the
LQ weather station, which was located about 5 km fromthe LQ study site. Water used for irrigation was estimatedby agricultural researchers or local farmers (based onelectricity use).
Mathematical Model
The Hydrus-1D software package (Šimůnek et al.
2005) was used to simulate 1D vertical flow assuming
applicability of the standard Richards equation for unsaturated flow as follows:



∂
∂h
∂θ (h)K(h)(1)=+ 1 − S(h)∂t∂z∂zsubjected to the initial and boundary conditionsh(z, 0) = h0 (z)


∂h
−K(h)
+ 1 
= q0 (t)
∂z


(2)
(3)

z=0

∂h
(L, t) = 0
∂z

(4a)

h(L, t) = hL (t)

(4b)

or

where θ is the volumetric water content [—], h is thepressure head [L], z is the vertical coordinate [L], assumed
to be 0 at the soil surface and directed upward, t is time[T], K is the unsaturated hydraulic conductivity [L T−1 ], S
is a sink term to account for root water uptake [T−1 ], h0 (z)is the initial condition, and q0 (t) is the fluid flux acrossthe soil surface boundary [L T−1 ]. As the lower boundary
condition (Equations 4a,b), we implemented either a freedrainage condition, ∂h/∂z = 0, when the water table wasrelatively deep (the LQ and LC sites), or a time variablypressure head condition, hL (t), when the water table wasrelatively close to the soil surface (the HS, DZ, and CZ
sites).
Following Feddes et al. (1978), root water uptake wasmodeled using the function
S(h) = λ(z)α(h)Tp(5)

where λ(z) is the relative root distribution function[—], Tp is the potential transpiration rate [L T−1 ], andNGWA.orgα(h) is a dimensionless water stress response function(0 ≤ α ≤ 1) to account for reductions in uptake due todrought stress [—].
For α(h), we adopted the functional form ofFeddes et al. (1978) using Hydrus-1D default parametervalues for wheat and maize based on studies by Taylorand Ashcroft (1972) and Wesseling (1991). This function assumes that within some pressure head range h3 ≤h ≤ h2 root water uptake is optimal. Below h3 , root
water uptake declines from 1.0 to 0 at h4 (the wilting point) due to water stress. Above h2 , root water
uptake similarly declines linearly from 1.0 to 0 at hl
due to insufficient aeration. The critical pressure head h3
increases for higher potential transpiration rates, Tp . The
Hydrus-1D default parameters of wheat used in thisstudy were: h1 = 0 cm, h2 = −1 cm, h3l = −500 cm,
h3h = −900 cm, h4 = −16,000 cm. The default parameters of maize were: h1 = −15 cm, h2 = −30 cm, h3l =−325 cm, h3h = −600 cm, h4 = −8000 cm (Wesseling1991).
For λ(z), we used a linear root distribution varying
between 0 at the bottom of the root zone and 1.0 at
the soil surface (Feddes et al. 2001; Mertens et al. 2006).
Observed plant root data in the Hebei plain indicated that,on average, the rooting depth of winter wheat increased
from approximately 15 cm in March to 75 cm in May,and the rooting depth of summer maize from 0 cm onJune 6 to 15 cm on July 2, and to 40 cm on August 1,
and to a final 45 cm on September 1.
The surface boundary Equation 3 was implemented
as an atmospheric condition without surface ponding inwhich q0 (t) equals precipitation plus irrigation minus
potential evaporation as long as the pressure head determined at the soil surface exceeds some minimum negativevalue (−10, 000 cm in our study). However, the boundary flux becomes a time-variable flux during periods ofhigh evaporation as determined by the pressure head distribution near the soil surface with h(0, t) = −10,000 cm.
Surface runoff was assumed to occur when the surface
becomes saturated, in which q0 (t) in Equation 3 decreases
in value. Equation 3 was implemented using daily precipitation, irrigation, and ET rates.
Simulations based on Equation 1 were applied to thefive sites using a self-adjusting numerical time stepping
scheme, and with finite element discretizations in which
nodal spacings generally ranged from 1 cm near the soil
surface up to a maximum of 10 cm at the bottom of
the flow domain for deep water tables (the DZ site).
The simulation period for model calibration was fromJanuary 1, 2003 to August 31, 2005, whereas groundwater
recharge rates at the sites were calculated from January
1, 2004 to December 31, 2004, using initial conditionsof January 1, 2003. Annual (2004) rainfall and irrigation
rates at the sites are given in Table 2.
Unsaturated Soil Hydraulic Properties
The solution of Equations 1 through 4 requires
estimates of the soil water retention, θ (h), and unsaturated
hydraulic conductivity, K(h), relationships. We used the
X. Lu et al. GROUND WATER

3

reference crop evapotranspiration rate, ET0 (mm/d), as
follows:

Table 2
Rainfall and Irrigation in 2004 at the Five Sites
Items

LQ

LC

HS

DZ

CZ

Rainfall (mm)Irrigation (mm)

434.6
375

456.7
300

434.6
225

534.5
300

447.8
225

van Genuchten (1980) soil hydraulic functions

θ (h) =⎧⎪
⎨ θr +⎪⎩θ −θs[1 + |αh|n ]m

(6)
h≥0

θs
1/2

K(h) =

h < 0 (m = 1−1/n)

1/m m 2

Ks Se [1 − (1 − Se

) ]

h<0
h≥0

Ks

(7)

where θr and θs are the residual and saturated watercontents [—], respectively, α [L−1 ] and n [—] areempirical shape factors that depend on soil type, Ks is
the saturated hydraulic conductivity [L T−1 ], and Se iseffective saturation given by
Se (h) =θ (h) − θrθs − θr(8)

Soil water retention data were measured in the laboratoryon undisturbed soil cores taken from the field sites(Lu et al. 2006). Measurements were made for pressure
heads between 0 and about −800 cm, whereas watercontents (including θs ) were measured gravimetrically.
The parameters θr , α, and n were subsequently fitted tothe observed retention data using the RETC program ofvan Genuchten et al. (1991).
Measurements of the vertical saturated hydraulic conductivity, Ks , were obtained in the laboratory using fallinghead experiments on undisturbed soil cores taken fromseveral depths at the LC and HS sites. For the LC site,
we obtained Ks values of 42.0, 28.6, and 41.4 cm/d for the
0 to 20, 35 to 50, and 280 to 310 cm soil layers, respectively. For the HS site, measured Ks values were 143.6,
35.3, and 84.8 cm/d, respectively, for cores taken from
depths of 0 to 30, 50 to 70, and 200 to 220 cm, respectively. Saturated hydraulic conductivities of the othersites were estimated from measured particle distributions
and soil textural averages using the Rosetta pedotransfer functions of Schaap et al. (1998) as implemented in
Hydrus-1D. All the laboratory-measured hydraulic parameters were considered to be initial estimates subjected tofurther calibration against observed soil water content timeseries measured at the field sites.
Evapotranspiration
In this study, we used the Penman-Monteith equation(Allen et al. 1998; Liu et al. 1997) to estimate the4

X. Lu et al. GROUND WATER

ET0 =

900
γ u2 (ea − ed )T + 273
 + γ (1 + 0.34u2 )0.408(Rn − G) +(9)

where Rn is net solar radiation (J/m2 /d), G is soil heat flux
(J/m2 /d), ea is the saturation water vapor pressure at mean
air temperature (kPa), ed is the mean actual vapor pressure (kPa),  is the slope of the saturation vapor pressure(kPa/◦ C), γ is the psychrometric constant (kPa/◦ C), andu2 is wind speed at a height of 2 m (m/s). Potential ETp(mm/d) was calculated using the formula:
ETp = Kc ·ET0(10)where Kc is a dimensionless crop coefficient. As crop
coefficients for both the HS and DZ sites of the alluvial
and lacustrine plains we used values estimated by Chenand Guo (1993) at a geomorphologically similar site near
Gaocheng in Hebei Province. The crop coefficients forthe experimental sites at LC and LQ were based ondata obtained by Zhang (2002) for a site near LC in the
piedmont plain (Table S1, Supporting Information).
The potential evaporation rate of a soil under astanding crop was derived from the Penman-Monteithequation by neglecting the aerodynamic term. The only
source of energy for soil evaporation is then net radiationthat reaches the soil surface (Ritchie 1972). Assumingthat net radiation inside a canopy decreases according to
an exponential function, and that the soil heat flux can beneglected, we can derive (Goudriaan 1977; Belmans et al.
1983):
E = Ep e−Kgr LAI(11)

where E is the soil evaporation [—], Kgr is the extinction
coefficient for solar radiation [—], and LAI is the Area
Index of wheat or maize [—]. Ritchie (1972) and Feddes(1978) both used Kgr = 0.39.
Daily potential transpiration rates (Tp ) were calculated by subtracting evaporation from total ET. Estimatedpotential evaporation and transpiration rates at the fivesites are shown in Figure S1, Supporting Information.

Results and Discussion
Simulations of the Field Soil Moisture Dynamics
Initial simulations using the laboratory estimatedunsaturated soil hydraulic parameters, or those estimated using pedotransfer functions, were found to berelatively poor and in need of improvement. For thesereasons, we used the Hydrus-1D inverse module (Šimůnek
et al. 2005) at all the sites to improve estimates of
the hydraulic parameters in terms of providing optimal
descriptions of the field data. Because of the large
number of parameters involved for the layered profiles,
many of them correlated, values of selected parameters(notably the parameter n) were adjusted only when the
NGWA.org

Table 3
Calibrated Soil Hydraulic Parameters of the Five Sites
Representative Sites

Depth (cm)Materials

LQ
LC

0–500
0–42.542.5–90
90–3300
0–50
50–150
150–400
0–140
140–5100–300

Silt
Silty clay
Silty clay
Silty clay and silt
Silt
Clay
Silty clay
Silt
Silt
Silty clay

HS

DZ

soil moisture

CZ

0.45

θr

θs

α

n

Ks (cm/d)

0.068
0.086
0.089
0.099
0.068
0.110
0.080
0.090
0.068
0.086

0.400
0.360
0.361
0.365
0.400
0.410
0.390
0.360
0.400
0.360

0.00100.01340.03150.01510.00200.01100.00700.00600.00100.0134

2.400
1.574
1.179
1.747
2.320
1.480
2.060
1.660
2.400
1.574

69.0
33.5
28.7
23.3
93.3
30.5
82.0
50.0
69.0
33.5

40cm

0.35
0.25
time(d)
0.15
0

50

100

150

200

soil moisture

0.45

250

300

350

400

450

80cm

0.35
0.25
time(d)
0.15
0

50

100

150

200

250

300

350

400

450

soil moisture

0.45
100cm
0.35
0.25
time(d)
0.15
0

50

100

150

200

250

300

350

400

450

soil moisture

0.45
120cm
0.35
0.25
time(d)
0.15
0

50

100

150

200

250

300

350

400

450

Figure 2. Comparison between measured (dots) and simulated (lines) volumetric water contents of different depths at the
Luancheng site.

differences between the calculated and field-measured
water contents were larger than the standard deviation.
The final identified parameters for the five sites (Table 3)
provided a close fit between simulated and measured watercontents, as illustrated for the LC site in Figure 2 using
data from four depths.
Temporal and Spatial Changes in Recharge
Temporal variations in the vertical groundwater
recharge rate depend mainly on the temporal distributions
NGWA.orgof water input (precipitation and irrigation) and ET.
Figure 3 shows the calculated daily recharge rate as wellas rainfall plus irrigation at the LC site from January 1,
2003 to August 31, 2005. The recharge rate from January to March was less than 1 mm/d due to low rainfall
and the absence of any irrigation. Recharge then gradually
increased from April to September because of irrigationin April and May, and because of more rainfall from
June to September. Starting at the end of September, therecharge rate decreased again. Figure 3 also indicates that
X. Lu et al. GROUND WATER

5

P+I

8

Recharge

75

6

50

4

25

2

0
1-Jan-03

Recharge(mm)

P+I(mm)100

031-May-03

28-Oct-03

26-Mar-04

23-Aug-04

20-Jan-05

19-Jun-05

Figure 3. Daily recharge and rainfall plus irrigation rates for the Luancheng site (from January 1, 2003 to August 31, 2005).

Table 4
Summary of 2004 Simulation Results Obtained for the Five Study Sites
Geomorphology

Water Table Depth (m)Site

R (mm)Rr (mm/d)

Rc (%)

Piedmont plain

10-40

Alluvial and lacustrine plains

2–20

Coastal plain

<5

LQ
LC
HS
DZ
CZ

169
180
140
155
102

0.46
0.49
0.38
0.43
0.28

21
24
21
19
15

Note: R is the total 2004 recharge, Rr is the average daily groundwater recharge rate, Rc is the recharge coefficient (%).

the recharge peaks are delayed relative to the maxima
in precipitation and irrigation. Antecedent moisture conditions and the time period over which rainfall occurs
appear equally or more important than the total quantityof rainfall. For example, the intense rain storms generated more recharge during the summer monsoon season
than did the totality of smaller precipitation events andirrigation applications during the remainder of the year.
Figure S2 (Supporting Information) shows daily
recharge rates vs. precipitation plus irrigation (P + I ) at
all the five sites for 2004. The LQ and LC sites in the
piedmont plain, both having deep water tables (Table 1),
showed the longest time-lags between the rainfall and irrigation events and recharge. The HS and DZ sites in the
alluvial and lacustrine plains, where water tables wereintermediate in depth, showed much shorter time-lags
between the rainfall and irrigation events and recharge.
Recharge at the CZ site in the coastal plain showed no
obvious time-lag after a rainfall or an irrigation event.
Recharge of groundwater in the more humid piedmont
plain was always larger than that in the alluvial and lacustrine plains and the coastal plain.
The simulation results for 2004 are further summarized in Table 4. They show that both the average recharge
rate Rr and the recharge coefficient Rc reduce graduallyfrom the piedmont plain to the coastal plain.
Relation Between Precipitation and Recharge
in the Hebei Plain
The results in Figures 3 and S2 (Supporting Information) indicate that the total amount of recharge throughthe unsaturated zone is mostly determined by the totalityof the local precipitation, irrigation, and ET rates, whereas
water table depth and the lithology of the unsaturated
6

X. Lu et al. GROUND WATER

zone (and hence soil moisture storage in the unsaturatedzone) determine the time delays and the degree of smoothing in the recharge rates relative to the precipitation andirrigation events. For example, the water table in CZ isrelatively shallow, which causes the infiltrated water toreach the water table relatively quickly so that most of
the individual rainfall/irrigation corresponds to isolated
infiltration recharge events with very small time-lags (generally only about 1 or 2 days).
Following Wu et al. (1996), we used the peak timelag here as the time elapsed from the beginning of
an individual rainfall event to the time at which therecharge rate reaches its maximum, and the cessation timelag as the time elapsed from the beginning of rainfall
to the cessation of recharge. As the depth to watertable increased, the correspondence between rainfall andrecharge decreased. The water table depths at the DZ
and HS sites were in the order of 2 to 5 m, whichcaused the recharge events produced by several individual
rainfalls to merge into one single process, although a
few smaller recharge peaks still corresponded to some of
the larger rainfall events or concentrated rainfall clusters.
A peak in the annual recharge process represents an
integrated recharge event produced by a set of individual
rains clustered around a large rainfall event corresponding
to this peak. The peak time-lags at DZ and HS weregenerally between about 3 and 5 days.
By comparison, the water tables at the LQ and LC
sites are relatively deep, thus requiring much more time
for the infiltrated water to reach the water table. The peaktime-lags at these sites were between 18 and 35 days ormore. The peak at LQ in late March was different fromthose at the other sites due to early irrigation and having
higher initial soil water contents. Recharge after the peak
NGWA.org

Table 5Comparisons of Calculated 2004 Recharge Rates Obtained with and Without Accounting for Irrigation

Study Sites
LQ
LC
HS
DZ
CZ

A (Rainfall + Irrigation)
Rr (mm/year)

B (Rainfall)Rr (mm/year)

C = A − B (Irrigation Only)
Rr (mm/year)

C /A (%)

169
180
140
155
102

86
103
83
126
75

83
77
57
29
27

49
43
41
20
28

decreased smoothly, which is typical for locations with adeep water table.

40

A

Simulated

TracedComparison of Numerical Results with Chemical Tracer
Studies
Wang et al. (2008a) previously used vadose zone
chemical (tritium and bromide) tracing to estimate groundwater recharge at similar locations in the Hebei plain asa function of soil type, land use, and irrigation. Figure 4
compares their results for the recharge coefficient Rc , with
our numerical simulations for different areas in the Hebei
plain. The two approaches give approximately the same
results for flood-irrigated croplands, with the locationsin the piedmont plain (LC and DZ) showing the highest recharge coefficients. The nonirrigated uncultivatedareas showed much lower Rc values, not only in the
coastal plain to the east (CZ), but also DZ in the central alluvial and lacustrine plains. It is to be noted thatthe tracer method gave somewhat higher estimates for therecharge coefficient as compared to the numerical modeling approach.
As compared to the tracer method, numerical simulation is a far more convenient and cost-effective approach,but it requires relatively long-term weather data, information about the local irrigation practices, and soil physicalparameters. However, if the data are available, numerical simulation may well be a very attractive alternative tochemical tracer methods.
NGWA.org

20
10
0
LQ

LC

DZ

CZ

B

Simulated

Traced

LC

DZ

40
30
Rc(%)

Analysis of Groundwater Recharge from Irrigation
Return FlowTo estimate groundwater recharge from irrigation
return flow, we reran the previously calibrated models forthe five sites assuming the same climate conditions, but
now without irrigation. We assumed that the recharge simulated with both irrigation and rainfall minus the recharge
simulated with rainfall only is the recharge stemming from
irrigation return flow. The simulation results (Table 5)show that groundwater recharge decreased significantly
and that recharge caused by irrigation in the Hebei plain
accounted for 27% to 49% of the total precipitation andirrigation. Irrigation-induced recharge at LQ was largest
among the five sites (49% of the total precipitation andirrigation). Recharge caused by irrigation in the piedmont
plain was slightly larger than that in the alluvial and lacustrine plains.

Rc(%)

30

20
10
0
CZ

Figure 4. Comparison of the recharge coefficient (Rc ) at different locations in the Hebei plain obtained with numericalunsaturated flow modeling (this study) and chemical tracing:(A) flood-irrigated cropland and (B) nonirrigated uncultivated land.

Uncertainty and Sensitivity Assessment
The vadose zone modeling approach used in thisstudy contains several potential sources of uncertainty.
Uncertainty exists about the computed infiltration fluxes,
ET fluxes, and recharge fluxes. Possible sources of uncertainty are the values of the many model parameters
needed to run Hydrus-1D, including the soil hydraulicparameters, parameters in the reference ET rate equation,
crop coefficients, and drought stress parameters in the
uptake model of Feddes et al. (1978). For the soil
hydraulic parameters, information is available in theform of 95% confidence intervals (CI) and correlation coefficients computed using the inverse solution
module of Hydrus-1D. Unfortunately, no such information is available about uncertainty in the other model
parameters.
X. Lu et al. GROUND WATER

7

Hence, as previously indicated by Jimenez-Martinez
et al. (2009) it does not seem possible at this time to placereliable quantitative CI on the computed recharge rates
due to a lack of knowledge about parameter variability andcorrelation structure. Still, some notion of the importanceof parameter uncertainty may be obtained by carrying out
a sensitivity analysis of various model parameters. Here,
we limit the analysis to the soil hydraulic parameters inEquations 6 and 7. Our calculations showed that the orderof sensitivity of the parameters for groundwater recharge
were: n, θs , Ks , α, and θr . The residual water content,
θr , was found to be the least sensitive, consistent withfindings by Jimenez-Martinez et al. (2009). Taking the
LC site as an example, allowing individual parameters
for all layers to decrease or increase by 50% from thecalibrated values in Table 3 caused groundwater recharge
to change by −18% and +15% for Ks , −8% and +6%
for α, and +2% for θr , respectively. When decreasing
the value of n by 10% or increasing by 10% and 50%,
groundwater recharge changed to −24%, +26%, and+59%, respectively. Except those for n, these changes are
relatively small and somewhat expected for simulations
in which rainfall, irrigation and evaporation rates aremostly imposed as boundary conditions (flux-controlled
infiltration), with still important but lesser effects due towater stress. The combined effects of these parametersreflect the fact that fine-textured soils (i.e., those having
lower n, Ks , and α values) retain more moisture near thesoil surface for subsequent root water uptake and ET, thus
generally leading to lower recharge rates (e.g., KennethSmith et al. 1994; Nolan et al. 2003; Jimenez-Martinez
et al. 2009).
The sensitivity of ET was calculated by perturbing
an input variable by a fixed percentage at all thetimes during the calibration period. ET in this case wasindividually increased or decreased by 10%. The model
was subsequently run using the perturbed input to obtain
the resulting groundwater recharge values. Results showthat decreasing or increasing the value of ET by 10%
caused the groundwater recharge to change by +5% and−4%, respectively.
We note that accurate characterization of the prediction uncertainty is problematic because of a lack of knowledge about parameter variability and parameter correlation
structure. Future study aimed at quantifying uncertainty
in parameters would greatly benefit from more comprehensive analyses of uncertainty in groundwater recharge
calculations of the type provided in this study.

Conclusions
Hydrus-1D was found to be a very useful tool for
simulating groundwater recharge in the Hebei plain characterized by relatively little topographic relief, relatively
deep water tables, and insignificant snowmelt, and withavailable data limited to the lithology of the unsaturated
zone and laboratory-measured soil physical parameters.
Our study of recharge at the five representative zones permits several conclusions.
8

X. Lu et al. GROUND WATER

1. Temporal variations in the groundwater recharge
rate at the five sites were found to be significant.
Groundwater recharge was significantly delayed relative
to rainfall or irrigation in the piedmont plain (the LQ
and LC sites) and the alluvial and lacustrine plains (DZ
and HS), but not obviously so in the coastal plain (CZ).
The highest recharge rates to groundwater at the DZ,
HS, and LC sites occurred in August, whereas at LQ thehighest recharge rates occurred in late March and at CZ
in July.
2. The recharge rate was on average 175 mm/yearin the piedmont plain, and 133 mm/year in both the
alluvial and lacustrine plains and the coastal plain. The
average infiltration recharge coefficient in the piedmont
plain (0.22) was higher than those of the alluvial andlacustrine plains and the coastal plain (0.18 in both theplains).
3. The peak time-lag between an infiltration event(rainfall or irrigation) in the coastal plain with its shallow
water tables was generally at most 1 or 2 days. Peak timelags for sites in the central alluvial and lacustrine plains
were usually about 3 to 5 days, with the recharge peaks
often generated by a series of individual rainfalls clusteredaround some large rainfall event. By comparison, thepeak time-lags in the piedmont zone with its much
deeper water tables (generally >10 m) were about 18 to
35 days. Recharge rates in this zone were much more
gradual with time because of buffering by the deep
unsaturated zone. This implies that different time-lagscorresponding to different water table depths may need tobe considered when estimating or modeling groundwater
recharge.
4. Groundwater recharge from irrigation in the Hebei
plain accounted for about 27% to 49% of the total
precipitation plus irrigation. Of the five sites, LQ showed
the highest rate stemming from irrigation (49% of the total
precipitation plus irrigation).
5. As compared with the tracer methods, numericalsimulation seems a much more convenient and costeffective approach for estimating groundwater recharge. A
disadvantage of unsaturated flow modeling is the need forlong-term data about climate, irrigation practices, and soil
physical parameters. Still, the numerical approach may be
preferred if the required data are available.

Acknowledgments
This study was financed in part by the 863 Program of China (2007AA06Z337), the 973 Programof China (2010CB428802), the China Geology Survey(200310400035-1), the National Natural Science Foundation of China (40772155, 40472123), and CAPES, Brazil.
We would also like to thank all the members of the projectgroup for their fruitful collaboration and the China Scholarship Council. We gratefully acknowledge the reviewers (Dr. Phillip Stauffer and two anonymous reviewers)
for their constructive comments and suggestions on thearticle.
NGWA.org

Supporting Information
Additional Supporting Information may be found in
the online version of this article:
Table S1. Monthly crop coefficients for the Luancheng and Gaocheng experimental sites
Figure S1. Daily potential evaporation and transpiration rates of the five sites
Figure S2. Daily recharge and precipitation plus
irrigation (P + I ) rates in 2004 at the five sites
Please note: Wiley-Blackwell are not responsible forthe content or functionality of any supporting materials
supplied by the authors. Any queries (other than missing
material) should be directed to the corresponding author
for the article.

References
Allen, R.G., L.S. Pereira, D. Raes, and M. Smith. 1998. Cropevapotranspiration. Guidelines for computing crop waterrequirements. Irrigation and Drainage Paper No. 56. FAO,
Rome.
Belmans, C., J. Wesseling, and R.A. Feddes. 1983. Simulation
function of the water balance of a soil: SWATRE. Journalof Hydrology 63, 271–286.
Chen, M., Q. Yin, and T.Y. Zhao. 2008. Evaluation of evaporation in the Hai river basin (in Chinese). Design ofWater Resources and Hydraulic Engineering 27, no. 2:
11–12.
Chen, Y.M., and G.S. Guo. 1993. Study on Water Requirement
Isoline of the Main Crop of China. Beijing, China:
Agricultural Scientech Press (In Chinese).
Feddes, R.A., H. Hoff, M. Bruen, T. Dawson, P. de Rosnay,
R. Dirmeyer, R.B. Jackson, P. Kabat, A. Kleidon, A. Lilly,
and A.J. Pitman. 2001. Modeling root water uptake inhydrological and climate models. Bulletin of the American
Meteorological Society 82, no. 12: 2797–2809.
Feddes, R.A., P.J. Kowalik, and H. Zarandny. 1978. Simulationof Field Water Use and Crop Yield . New York: John Wiley
Sons.
Goudriaan, J. 1977. Crop meteorology: a simulation study.
Simulation monographs, Pudoc, Wageningen.
Jimenez-Martinez, J., T.H. Skaggs, M.Th. van Genuchten, and
L. Candela. 2009. A root zone modeling approach toestimating groundwater recharge from irrigated areas.
Journal of Hydrology 367, no. 1–2: 138–149.
Jin, M.G., R.Q. Zhang, L.F. Sun, and Y.F. Gao. 1999. Temporaland spatial soil water management: a case study in
the Heilonggang region, PR China. Agricultural Water
Management 42, no. 2: 174–187.
Jin, M.G., X. Liang, I. Simmers, Y.F. Gao, and R.Q. Zhang.
2000. Estimation of groundwater recharge using artificial
tritium tracing. Proceedings of the International Symposiumon Hydrogeology and the Environment, Wuhan, China,
17–21 October 2000. China: Environmental Science Press,
340–345.
Kenneth-Smith, A., P.G. Cook, and G.R. Walker. 1994. Factor
affecting groundwater recharge following clearing in thesouth western Murray Bason. Journal of Hydrology 154,
no. 1: 85–105.
Liu, Y., L.S. Pereira, and J.L. Teieirab. 1997. Update definitionand computation of reference evapotranspiration comparison with former method (In Chinese). Journal HydraulicEngineering 6, 27–33.

NGWA.org

Lu, X.H., M.G. Jin, and B.G. Wang. 2006. Discussion on thesoil water characteristic curve of the agricultural ecosystem
experiment station in Luancheng. Hebei China (In Chinese).
Rural Water and Hydropower 12, 30–32.
Mertens, J., R. Stenger, and G.F. Barkle. 2006. Multiobjective
inverse modeling for soil parameter estimation and model
verification. Vadose Zone Journal 5, no. 3: 917–933.
Nolan, B.T., A.L. Baehr, and L.J. Kauffman. 2003. Spatial
variability of groundwater recharge and its effect on shallow
groundwater quality in southern New Jersey. Vadose Zone
Journal 2, 677–691.
Ritchie, J.T. 1972. A model for predicting evaporation from a
row crop with incomplete cover. Water Resource Research8, no. 5: 1204–1213.
Scanlon, B., K. Keese, R.C. Reedy, J. Šimůnek, and B. Andraski.
2003. Variations in flow and transport in thick desert vadosezones in response to paleoclimatic forcing (0–90 kyr):monitoring, modeling, and uncertainties. Water Resource
Research 39, no. 7: 1179. DOI: 10.1029/2002WR001604.
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 1998.
Neural network analysis for hierarchical prediction ofsoil hydraulic properties. Soil Science Society of America
Journal 62, 847–855.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2008. Development and applications of the HYDRUS and STANMOD
software packages and related codes. Vadose Zone Journal
7, 587–600.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2005.
The HYDRUS-1D software package for simulating the
one-dimensional movement of water, heat, and multiple
solutes in variably-saturated media, Version 3.0, HYDRUS
Software Series 1, Department of Environmental Sciences,
University of California, Riverside, CA, 270p.
Taylor, S.A., and G.M. Ashcroft. 1972. Physical Edaphology,
434–435. San Francisco, California: Freeman and Co.
van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil
Science Society of America Journal 44, 892–898.
van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The
RETC code for quantifying the hydraulic functions of
unsaturated soils. Report No. EPA/600/2-91/065. Ada,
OK: R.S. Kerr Environmental Research Laboratory, U.S.
Environmental Protection Agency, 85p.
Wang, B., M. Jin, J.R. Nimmo, L. Yang, and W. Wang. 2008a.
Estimating groundwater recharge in Hebei plain, China
under varying land use practices using tritium and bromide
tracers. Journal of Hydrology 356, no. 1–2: 209–222.
Wang, S.Q., J.L. Shao, X.F. Song, Y.B. Zhang, Z.B. Huo, and
X.Y. Zhou. 2008b. Application of MODFLOW and geographic information system to groundwater flow simulationin North China Plain, China. Environmental Geology 55,
no. 7: 1449–1462.
Wesseling, J.G. 1991. Meerjarige simulates van groundwateronttrekking voor verschillende bodemprofielen, grondwatertrappen en gewassen met het model SWATRE. Report 152.
Winand Starring Centre, Wageningen.
Wu, J., R. Zhang, and J. Yang. 1996. Analysis of rainfallrecharge relationships. Journal of Hydrology 177, no. 1–2:
143–160.
Yang, Y., M. Watanabe, X. Zhang, X. Hao, and J. Zhang. 2006.
Estimation of groundwater use by crop production simulated by DSSAT-wheat and DSSAT-maize models in the
piedmont. Hydrological Processes 20, 2787–2802.
Zhang, X. 2002. Index system for irrigation scheduling of
winter wheat and maize in the Piedmont of mountain
Taihang (In Chinese). Transactions of the Chinese Societyof Agricultural Engineering 6, 35–39.

X. Lu et al. GROUND WATER

9