Journal of Japan Society of Civil Engineers, Ser.
B1 (HydraulicEngineering),
Vol. JSCE,
68, No.Vol.56,
4, I_211-I_216,
2012.
Annual Journalof HydraulicEngineering,
2012, February

GLOBAL SIMULATION OF GROUNDWATER
RECHARGE, WATER TABLE DEPTH, AND LOW
FLOW USING A LAND SURFACE MODEL WITHGROUNDWATER REPRESENTATION.
Sujan KOIRALA1, Hannah G. YAMADA2, Pat J.-F. YEH3, Taikan OKI4, Yukiko
HIRABAYASHI5, and Shinjiro KANAE6
1Member of JSCE, Dr. of Eng., Researcher, Department of Mechanical and Environmental Informatics, Tokyo
Institute of Technology (2-12-1-W8-4 O-okayama, Meguro-ku, Tokyo 152-8552, Japan)
2Member of JSCE, Master of Eng., Department of Mechanical and Environmental Informatics, Tokyo Institute ofTechnology (2-12-1-W8-4 O-okayama, Meguro-ku, Tokyo 152-8552, Japan)
3 Dr. of Eng., Associate Professor, Department of Civil Engineering, University of Tokyo (4-6-1 Komaba,
Meguro-ku, Tokyo 153-8505, Japan)
4Member of JSCE, Dr. of Eng., Professor, Department of Civil Engineering, University of Tokyo (4-6-1 Komaba,
Meguro-ku, Tokyo 153-8505, Japan)
5Member of JSCE, Dr. of Eng., Associate Professor, Institute of Engineering Innovation, University of Tokyo(2-11-16 Bunkyo-ku, Tokyo 113-8656, Japan)
6Member of JSCE, Dr. of Eng., Associate Professor, Department of Mechanical and Environmental Informatics,
Tokyo Institute of Technology (2-12-1-W8-4 O-okayama, Meguro-ku, Tokyo 152-8552, Japan)

With the withdrawal for anthropogenic uses in addition to climatic changes, the sustainability ofgroundwater resources is under question. Global-scale land surface models commonly used for water
resources assessment, however, simplify or completely neglect the groundwater processes making theminapplicable for groundwater resources assessment. In this study, a groundwater representation is
implemented into a global-scale LSM, the MATSIRO, enabling it to simulate the major groundwater
variables namely, groundwater recharge, water table depth, and low flow. The estimated global
groundwater recharge (29900 km3/yr) corresponds well with GSWP-2 mean baseflow (30200 km3/yr).
Global distribution of water table depth is found to be mainly controlled by climate and soil properties.
The comparison of simulated and observation-based daily flow duration curves in selected global river
basins reveals that the simulation of low flow improves significantly with the groundwater representation.

Key Words : global land surface model, groundwater resources, water table depth, low flow

1. INTRODUCTION
Land surface models (LSMs) used for the
global-scale hydrologic simulations usually neglect
or implicitly represent the groundwater process. Thebase runoff is often parameterized as a ‘free gravity
drainage’ from the soil water storage1) and hence themajor hydrologic variables related to thegroundwater process; e.g., base runoff, water table
depth (WTD) and groundwater recharge are notsimulated in a proper manner. To our knowledge,
there is no global-scale modeling study thatprovides the estimation of all major groundwater

variables. Existing global water resources
assessments are based on crude annual average due
to lack of proper estimation of seasonal variations of
runoff2). To reproduce the seasonal variation ofavailable water resources, a proper representation ofbase runoff or low flow is necessary3).
Furthermore, due to anthropogenic water uses,
the groundwater resources is depleting rapidly invarious regions of the world, e.g., central UnitedStates and northeast India4),5). Hence, the estimationof groundwater recharge is essential for global
assessment of groundwater resources as it indicatesthe naturally renewable groundwater resources6) but

I_211

Annual Journal of Hydraulic Engineering, JSCE, Vol.56, 2012, February
5),7)

its estimations on global scale are relatively few.
In previous model-based studies, recharge wassimulated either as a fraction of total runoff basedon slope relief, soil texture, hydrogeology,
permafrost, and precipitation intensity, while theupward capillary flux from shallow groundwater
was not considered7), or as the flux between thelowermost soil layer and groundwater storage, while
the dynamic interaction between them was not
explicitly modeled5). Due to these limitations, animprovement in the model representation ofgroundwater recharge is desirable.
Hence, in this study, a dynamic groundwater
representation is integrated into a global-scale LSM,
the Minimal Advanced Treatments of Surface
Integration and runoff (MATSIRO)8). With thegroundwater representation, the base runoff isgenerated from groundwater reservoir, WTD is
prognostically updated, and groundwater recharge isdirectly simulated based on Richards’ equation.

Fig.1 Structure of the soil column in the (a) original MATSIRO
and (b) MATSIRO with groundwater representationa grid, and Ls [L] is the length of a conceptual
hillslope and it is inversely proportional to tanβs,
and dgw [L] is WTD.
(2) MATSIRO with groundwater (MAT-GW)
A simple unconfined aquifer model10) is coupled
to the soil model of MAT-ORI. One-dimensional
lumped water balance equation for the unconfined
aquifer can be expressed as,

2. MODEL DESCRIPTIONSy

MATSIRO is the land surface scheme of an
Atmospheric Ocean General Circulation Model, theModel for Interdisciplinary Research On Climate(MIROC), jointly developed by the Atmosphere andOcean Research Institute at the University of Tokyo,
the National Institute of Environmental Studies, andthe Frontier Research Center for Global Change inJapan. Although the majority of hydrologic
processes are physically represented in MATSIRO,
it lacks a proper representation of groundwater
processes. The structures of soil column and baserunoff calculation in the original MATSIRO and
MATSIRO with groundwater representation arebriefly explained in the following.

K 0tan s
exp(1  f atn d gw )f atn Lst

 I gw  Qgw(2)

where Sy [-] is specific yield, Δdgw [L] is WTD, Δt[T] is time step, Igw [LT-1] is groundwater recharge
to (when positive) or capillary flux from (whennegative) the groundwater reservoir, and Qgw[LT-1] is baseflow. The fluctuation of dgw isgoverned by the difference of Igw and Qgw. Igw
depends on the degree of saturation of lowermost
soil layer and soil properties. It is calculated as,

I gw  G f  C f
Gf  k(3)

 d 
C f  k  dz 

(1) Original MATSIRO (MAT-ORI)
Soil column is divided into five layers (Fig. 1a)and the temperature and moisture (liquid andfrozen) are calculated for each layer. The thicknessof soil layers is 5, 20, 75, 100 and 200 cm
respectively.
A simplified version of TOPMODEL9) is
adopted to represent runoff process. Therefore,
WTD is implicitly calculated and does not directly
represent the boundary between unsaturated andsaturated soil zones. Baseflow is calculated as,

Qgw 

d gw

where Gf [LT-1] is gravity drainage flux (downward
to groundwater reservoir), Cf [LT-1] is capillary flux
(upward from groundwater reservoir), k [LT-1] is theunsaturated hydraulic conductivity of soil, dψ [L]and dz [L] is difference in matric potentials andelevation of saturated and unsaturated zones.
The size of the groundwater reservoir isdynamic with time and exact location of the watertable determines the number of soil layers inunsaturated zone for which the Richards’ equation issolved at each modeling time step. In order toaccommodate the variable WTD and accurately
locate its position, soil column is extended to ten
meters below the ground surface - in total 12 layerswith the depth of 5, 20, and 75 cm for the first threelayers, and 1 m each for the remaining nine layers.
A schematic representation of the linkage of the(1)

where Qgw [LT-1] is the baseflow, K0 [LT-1] is thesaturated hydraulic conductivity at ground surface,fatn [L-1] is the attenuation coefficient of K0 withdepth, tanβs [-] is the mean topographic slope within

2

I_212

Annual Journal of Hydraulic Engineering, JSCE, Vol.56, 2012, February

Fig.2 Global distribution map of long-term mean (1985-1999) of (a) groundwater recharge in mm/yr and (b) water table depth in m

soil-groundwater model in MAT-GW is presented inFig. 1b. Unlike MAT-ORI, the soil column now has
an explicit dynamic representation of unsaturated
and saturated zones separated by the water table,
which fluctuates with time.
The TOPMODEL-based baseflow in MAT-ORI(Eq.1) is replaced by the following threshold
relationship developed based on observations in
Illinois12),

Qgw  K  d 0  d gw  if
Qgw  0if

0  d gw  d0

was provided by the International Satellite Land
Surface Climatology Project- Initiative II(ISLSCP-2). The global distribution of vegetation
wasprovidedbytheInternational
Geosphere-Biosphere Programme, and vegetation
properties were provided by the University of Wales.
The parameters related to the groundwater model(d0, K, and Sy in Eq.(2) and Eq.(4)) were estimatedat global grid scale using climatic characteristics13).

4. RESULTS(4)d gw  d 0
In this section, the simulations of groundwater
recharge and WTD by MAT-GW are discussed first,
followed by a comparison of flow duration curves
simulated by the MAT-ORI, MAT-GW, andobservation by Global Runoff Data Centre.
Globally, mean (1985-1999) groundwater
recharge by MAT-GW is 29900 km3/yr. The
simulated global recharge volume is larger thanprevious estimates of 12700 km3/yr7) and 15200
km3/yr5). In the long term (without human
influences), however, the recharge should bebalanced by baseflow. The global groundwater
recharge by MAT-GW is closer to the GSWP-2
multimodel average baseflow (30200 km3/yr)12) thanprevious estimates5),7). The global distributions oflong-term mean groundwater recharge and WTD by
MAT-GW are presented in Fig. 2. Humid regions
like the Amazon, Congo and southeastern Asia havethe largest groundwater recharge (>1000 mm/yr),
while arid and semi-arid regions have small
recharge (Fig. 2(a)). In hot (semi-)arid regions, the-1where K [T ] is the outflow constant, and d0 [L] is
the threshold WTD above which baseflow is
initialized. Both d0 and dgw are taken as positive
values during calculation.

3. DATA DESCRIPTION
Simulation was driven using the global NCC
forcing dataset with 1ox1o spatial and 6-hourly
temporal resolutions11). Simulation is first carriedout for a 15-year (1985-1999) period, and theobtained climatologies of hydrologic states are used
to initialize another 15-year (1985-1999) simulation
for analysis. This procedure is deemed necessary
since it usually takes more than 10 years for thesimulated WTD to reach equilibrium in the arid
regions. The model time step is 1 hour. In addition
to forcing data, leaf area index and forest floor
albedo were taken from the Global Soil Wetness
Project (GSWP-2)12). The global distribution of soil
3

I_213

Annual Journal of Hydraulic Engineering, JSCE, Vol.56, 2012, February

Fig.3 Difference of the MAT-GW simulated groundwater recharge (mm/yr) and (a) FAO Aquastat and (b) previous
model-based estimate

7)

recharge from South America, while Congo river
basin contributes 54% of recharge from African
continent. Due to limited amount of observation
data for calibration of model in these regions, thereliability of previous recharge estimation7) is lowerand hence the uncertainty in prediction is relatively
higher. Nonetheless, the simulation seems to agree
fairly well in semi-arid and arid regions where thegroundwater recharge is limited and have large risk
to be under water stress in future with increasing
population as well as the climate change conditions.
Generally, spatial distribution of WTD iscontrolled by climate, soil and topographicalproperties15). A shallow WTD may reflect either
large infiltration, governed by climatic condition, or
poor drainage condition, governed by soil andtopographical characteristics. In Fig. 2(b), theglobal pattern of MAT-GW simulated WTD isfound to be mainly controlled by recharge,
baseflow, which are controlled by climate, as well
as soil properties.
The spatial distribution of the WTD correspondswith that of climatic characteristics indicated by
Budyko dryness index (ratio of net radiation tomean precipitation and has high values for dry
regions and low for humid regions) shown in Fig.
4(a). The WTD is deeper in dry regions whereas it
is shallower in humid regions. Further, thesimulated WTD pattern does not bear any
relationship with topography as topographic factors
are not considered in calculating the WTD in

precipitation is relatively low (0-50 mm/mon), of
which 70-85% is found to be evaporated, resulting
in small infiltration to soil and lower groundwater
recharge. Also, high latitude regions, with frozen ice
content in soil > 0.1 m/m, have relatively lower
recharge. If the frozen ice content in the root zone is0.10 m/m, the unsaturated hydraulic conductivity (kin Eq.3) is found to be decreased by 20-25%,
depending upon soil type, resulting in lower
recharge.
The validation of groundwater recharge and
WTD simulation by MAT-GW in the Illinois region
has been previously carried out14). On the
global-scale, however, there is no observation data
available. As MAT-ORI cannot simulate therecharge, the differences of only MAT-GW
recharge simulation and previous model-based
country-wise estimate7) and statistics-based Foodand Agriculture Organization AQUASTAT (FAO)
database is presented in Fig. 3. The model-based
estimate7) was constrained using observed riverdischarge and is reliable in regions with sufficient
observations5). Further, the FAO estimate is basedon statistical data collected from many countries
around the world. The difference between
MAT-GW and previous estimates is found to be thelargest in the humid countries while it is relatively
small in relatively drier as well as cold countries.
The largest difference can be seen in the countries
located within the humid river basins. In MAT-GW
simulation, Amazon basin contributes about 57% of

4

I_214

Annual Journal of Hydraulic Engineering, JSCE, Vol.56, 2012, February(Fig. 2(b)) are caused by the difference in soil
properties; for example, in the Amazon basin thegrid cells with a loamy clay soil usually have a
deeper WTD than the clay grid cells (Fig. 4(b)).
Under a similar climate, regions with clay or clayeysoils may have a shallower WTD compared to sandy
soils, as drainage in latter case is more efficient.
Finally, a comparison of flow duration curves(FDCs) generated from daily flows is presented. A
FDC is a plot that shows percentage of the time thatflow in a river is likely to equal or exceed somespecified value of interest. Daily FDC looks steeper
than monthly FDC and it clearly displays theextreme flows. Analysis of flow duration curves isneeded to investigate the availability of riverdischarge; especially in dry season, i.e., low flow. A
daily time series of one year was generated from
GRDC data available for the study period. Long
term mean was used to replace the missing data.
The daily FDCs for target river basins arepresented in Fig. 5. MAT-ORI cannot simulate low
flow correctly for majority of basins. The difference
between MAT-ORI and MAT-GW is small for highlatitude river basins (e.g., Lena and Kolyma river
basins). For humid basins, the MAT-GW simulationof 90th percentile flow (the flow available in 9 out of
10 days) is much closer to observation than the
MAT-ORI simulation. Similar improvement can be

Fig.4 Global distribution of (a) Bukyo dryness index (-) and(b) soil types

MAT-GW (Eq.2). Streamflow network (drainage
density) and other hillslope terrain attributes also
have critical influences on baseflow and WTD.
However, these effects work at smaller spatial scales
and their control at the grid-size scale of the global
simulation may be secondary. As in most LSMs,
topography and lateral flow between grid cells are
not explicitly considered in MATSIRO. Additional
smaller-scale heterogeneities in WTD distribution

Fig.5 Comparison of daily flow duration curves for MAT-GW simulation (dashed lines), MAT-ORI simulation(dotted lines), and observation (solid lines)
5

I_215

Annual Journal of Hydraulic Engineering, JSCE, Vol.56, 2012, February

seen in dry basins (e.g., Darling, Orange, andZambezi) where surface runoff dominated
MAT-ORI simulation cannot produce sufficientbase runoff in dry season as there is no
parameterization of groundwater reservoir. For
Congo basin, none of the parameterizations canreproduce the long term mean observed flow
duration characteristic correctly.

REFERENCES
1) Lohmann, D., Lettenmaier, D. P., Liang, X., Wood, E. F,
Boone, A. and Chang, S., et al.: The Project forIntercomparison of Land-surface Parameterization Schemes(PILPS) phase 2(c) Red-Arkansas River basin experiment::
3. Spatial and temporal analysis of water fluxes, Global andPlanetary Change, pp.161-179, 1998.
2) Oki, T. and Kanae, S.: Global Hydrologic Cycle and World
Water Resources, Science, Vol. 313, pp.1068-1072, 2006.
3) Maxwell, R. M. and Miller, N. L.: Development of a
Coupled Land Surface and Groundwater Model, Journal of
Hydrometeorology, Vol.6, pp.233-247, 2005.
4) Rodell, M., Velicogna, I. and Famiglietti, J. S.:
Satellite-based estimates of groundwater depletion in India,
Nature, Vol. 460, pp. 999-1002, 2009.
5) Wada, Y., van Beek, L. P. H., van Kempen, C. M.,
Reckman, J. W. T. M., Vasak, S. and Bierkens, M. F. P.:
Global depletion of groundwater resources, GeophysicalResearch Letters, Vol. 37(20), pp.L20402, 2010.
6) Doll, P.: Vulnerability to the impact of climate change on
renewable groundwater resources: a global-scale assessment,
Environmental Research Letters, Vol.4, pp.035006, 2009.
7) Doll, P. and Fiedler, K.: Global-scale modeling ofgroundwater recharge, Hydrology and Earth System
Sciences, Vol. 12, pp.863-865, 2008.
8) Takata, K., Emori, S. and Watanabe, T.: Development ofminimal advanced treatments of surface interaction andrunoff, Global and Planetary Change, pp.209-22, 2003.
9) Stieglitz, M., Rind, D., Famiglietti, J. S. and Rosenzweig,
C.: An Efficient Approach to Modeling the Topographic
Control of Surface Hydrology for Regional and Global
Climate Modeling, Journal of Climate, Vol.10(1), pp.
118-137, 1997.
10) Yeh, P. J.-F. and Eltahir, E. A. B.: Representation of watertable dynamics in a land surface scheme. Part I: Modeldevelopment, Journal of Climate, Vol.18, pp.1861-1880,
200511) Ngo-Duc, T., Polcher, J. and Laval, K.: A 53-year forcing
data for land surface models, Journal of GeophysicalResearch, Vol.110 (D06116), 2005.
12) Dirmeyer, P. A., Gao, X., Zhao, M., Zhichang, G., Oki, T.
and Hanasaki, N.: GSWP-2: Multimodel Analysis and
Implications for Our Perception of the Land Surface,
Bulletin of American Meteorological Society, Vol. 87(10),pp.1391-1397, 2006.
13) Koirala, S., Yeh, P. J.-F., Oki, T. and Kanae, S.: Parameter
estimation of a groundwater representation applicable in a
global-scale land surface model, Annual Conference ofJapan Society of Hydrology and Water Resources, 2010.
14) Koirala, S., Yeh, P. J.-F., Oki, T. and Kanae, S.: Fully
Dynamic Groundwater Representation in the MATSIROLand Surface Model, Annual Journal of HydraulicEngineering, Japan Society of Civil Engineers, Vol.54,
2010.
15) York, J.P., Person, M. and Gutowski, W.J.: Putting aquifers
into atmospheric simulation models: an example from theMill Creek Watershed, northeastern Kansas, Advances in
Water Resources, Vol. 25(2), 2002.
(Received September 30, 2011)

5. SUMMARY AND DISCUSSIONS
To our knowledge, a first comprehensive global
land surface hydrologic simulations of all the major
groundwater related variables were presented.
Groundwater recharge, which indicates the potential
groundwater resources that can be used without
depleting the source, was estimated to be 29900
km3/yr globally. Humid regions have the largest
recharge while the dry and cold regions have thelowest. Compared to previous estimates, therecharge was higher in the humid region for
MAT-GW simulation but the amount was similar inother regions suggesting acceptable simulations.
This implies that the physically-based calculation ofgroundwater recharge using Richards’ equation can
reproduce the previous estimates using a calibrated
conceptual model. The conceptual models are
considered relatively poorer for assessment of water
resources under climate change, when observation is
not available for model calibration. Hence,
estimation of groundwater variables using aphysically-based land surface model enhances theability for future assessment of water resources.
Further, global WTD distribution was found to bemainly controlled by the climate while secondary
control was provided by the soil properties. The
result indicates that the effect of topography might
be more pronounced in local scale and cannot be
observed directly in grid scale of 100’s of km.
Finally, the simulation of river discharge was used
to calculate the daily flow duration curves. The
comparison of flow duration curves for MAT-ORI
and MAT-GW with the observation reveals asignificant improvement in the simulation of low
flow (90th percentile value), especially in the humid
regions where base runoff is the significant runoff
generation mechanism, and dry regions, where mostof the dry season flow is from groundwater storage.
ACKNOWLEDGMENT: The research was
supported by CREST project of Japan Science andTechnology Agency (JST) and also partially
supported by JSPS- KAKENHI, Grant-in-Aid for
Scientific Research (S)(19106008), and Funding
program for next generation world-leading
researchers, JSPS.

6

I_216