Published June 24, 2019

Reviews and Analyses
Core Ideas• Land surface models (LSMs) showa large variety in describing andupscaling infiltration.
• Soil structural effects on infiltrationin LSMs are mostly neglected.
• New soil databases may help toparameterize infiltration processes
in LSMs.
Received 22 Oct. 2018.
Accepted 16 Feb. 2019.
Citation: Vereecken, H., L. Weihermüller, S.
Assouline, J. Šimůnek, A. Verhoef, M. Herbst,
N. Archer, B. Mohanty, C. Montzka, J. Vanderborght, G. Balsamo, M. Bechtold, A. Boone,
S. Chadburn, M. Cuntz, B. Decharme, A.
Ducharne, M. Ek, S. Garrigues, K. Goergen, J.
Ingwersen, S. Kollet, D.M. Lawrence, Q. Li, D.
Or, S. Swenson, P. de Vrese, R. Walko, Y. Wu,
and Y. Xue. 2019. Infiltration from the pedon
to global grid scales: An overview and outlook for land surface modeling. Vadose Zone
J. 18:180191. doi:10.2136/vzj2018.10.0191

Affiliations: H. Vereecken, L. Weihermüller,
M. Herbst, C. Montzka, J. Vanderborght, K.
Goergen, and S. Kollet, Agrosphere InstituteIBG-3, Forschungszentrum Jülich GmbH,52425 Jülich, Germany; H. Vereecken, International Soil Modeling Consortium (ISMC),
Forschungszentrum Jülich GmbH, 52425
Jülich, Germany; S. Assouline, Dep. of Environmental Physics and Irrigation, Institute ofSoils, Water and Environment Sciences, A.R.O.
Volcani Center, Bet Dagan, Israel; J. Šimůnek,Dep. of Environmental Sciences, Univ. of California, Riverside, CA 92521; A. Verhoef, Dep.
of Geography and Environmental Science,
Univ. of Reading, Reading, RG6 6AB, UK; N.
Archer, British Geological Survey, The Lyell
Centre, Edinburgh, EH14 4AP, Scotland; B.
Mohanty, Dep. of Biological and AgriculturalEngineering, Texas A&M Univ., College Station, TX 77843; G. Balsamo, European Centre
for Medium-Range Weather Forecasts, Reading, RG2 9AX, UK; M. Bechtold, Dep. of Earthand Environmental Sciences and Dep. of
Computer Science, KU Leuven, 3001 Heverlee, Belgium; A. Boone and B. Decharme,
CNRM, UMR3589 (Météo-France, CNRS), Toulouse, France; S. Chadburn, School of Earthand Environment, Univ. of Leeds, Leeds,
LS2 9JT, UK; S. Chadburn, Univ. of Exeter,
Exeter, EX4 4QE, UK; M. Cuntz, Univ. de Lorraine, INRA, AgroParisTech, UMR Silva, 54000
Nancy, France; A. Ducharne, Sorbonne Univ.,
CNRS, EPHE, UMR 7619 METIS, Paris, France;(continued on next page)© 2019 The Author(s). This is an open access article
distributed under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-ncnd/4.0/).

Infiltration from the Pedon to Global
Grid Scales: An Overview and Outlook
for Land Surface Modeling
Harry Vereecken, Lutz Weihermüller, Shmuel Assouline,
Jirka Šimůnek, Anne Verhoef, Michael Herbst, Nicole Archer,
Binayak Mohanty, Carsten Montzka, Jan Vanderborght,
Gianpaolo Balsamo, Michel Bechtold, Aaron Boone,
Sarah Chadburn, Matthias Cuntz, Bertrand Decharme,
Agnès Ducharne, Michael Ek, Sebastien Garrigues,
Klaus Goergen, Joachim Ingwersen, Stefan Kollet,
David M. Lawrence, Qian Li, Dani Or, Sean Swenson,
Philipp de Vrese, Robert Walko, Yihua Wu,
and Yongkang Xue
Infiltration in soils is a key process that partitions precipitation at the land surface into surface runoff and water that enters the soil profile. We reviewed thebasic principles of water infiltration in soils and we analyzed approaches commonly used in land surface models (LSMs) to quantify infiltration as well as its
numerical implementation and sensitivity to model parameters. We reviewed
methods to upscale infiltration from the point to the field, hillslope, and grid cell
scales of LSMs. Despite the progress that has been made, upscaling of local-scale
infiltration processes to the grid scale used in LSMs is still far from being treatedrigorously. We still lack a consistent theoretical framework to predict effective
fluxes and parameters that control infiltration in LSMs. Our analysis shows that
there is a large variety of approaches used to estimate soil hydraulic properties.
Novel, highly resolved soil information at higher resolutions than the grid scale of
LSMs may help in better quantifying subgrid variability of key infiltration parameters. Currently, only a few LSMs consider the impact of soil structure on soil
hydraulic properties. Finally, we identified several processes not yet consideredin LSMs that are known to strongly influence infiltration. Especially, the impact
of soil structure on infiltration requires further research. To tackle these challenges and integrate current knowledge on soil processes affecting infiltration
processes into LSMs, we advocate a stronger exchange and scientific interaction
between the soil and the land surface modeling communities.
Abbreviations: BC, boundary condition; CLM, common land model; FTC, freeze–thaw cycle; GCM, global
climate model; LSM, land surface model; pdf, probability density function; PDM, probability-distributedmodel; PTF, pedotransfer function; SWR, soil water repellency.

Infiltration or water entry into the soil profile is a key process in the hydrological cycle.
Its rate and dynamics affect the partitioning of precipitation at the land surface and determines the onset of ponding and, consequently, the formation of overland flow and runoff.
Infiltration affects irrigation efficiency on managed lands and the resulting stored soil
water available to vegetation (e.g., Verhoef and Egea-Cegarra, 2013), overland flow and soil
erosion processes (e.g., Assouline and Ben-Hur, 2006; Garrote and Bras, 1995; Poesen et al.,2003), groundwater recharge (e.g., Dahan et al., 2008), the exchange of water and energy
between the soil and atmosphere by controlling soil water content at the surface (Kim et al.,2017; MacDonald et al., 2017), and with this the flux partitioning into latent and sensible
heat flux with multiple atmospheric feedbacks (e.g., Keune et al., 2016; Seneviratne et al.,2010), stream flow and flooding events (Garrote and Bras, 1995), and various soil physical processes such as the onset of landslides (Lehmann and Or, 2012) and soil mechanical
Vadose Zone Journal | Advancing Critical Zone ScienceM. Ek, Environmental Modeling Center, National Centers for Environmental Prediction,
College Park, MD 20740; S. Garrigues, French
National Institute for Agricultural Research(INRA) CESBIO, Toulouse, France, and Centre
for Ecology and Hydrology, CEH, Wallingford,
UK; K. Goergen, Centre for High-Performance
Scientific Computing in Terrestrial Systems,
Geoverbund ABC/J, 52425 Jülich, Germany; J.
Ingwersen, Institute of Soil Science and Land
Evaluation, Biogeophysics, Univ. Hohenheim,
70599 Stuttgart, Germany; D.M. Lawrence,
National Center for Atmospheric Research,
Boulder, CO 80305; Q. Li, Center for Monsoon
Research, Institute of Atmospheric Physics,
Chinese Academy of Sciences, Beijing 100029,
China; D. Or, ETH Zurich, Institut für Biogeochemie und Schadstoffdynamik, 8092 Zürich,
Switzerland; S. Swenson, Climate and Global
Dynamics Lab., Terrestrial Sciences, NCAR,
Boulder, CO 80307; P. de Vrese, Max-Planck
Institut for Meteorology, Land in the Earth
System, 20146 Hamburg, Germany; R. Walko,
Rosenstiel School of Marine and Atmospheric
Science, Univ. of Miami, Miami, FL 33149; Y.
Wu, Environmental Modeling Center, National Centers for Environmental Prediction and
I.M. System Group at NCEP, College Park, MD
20740; Y. Xue, Dep. of Atmospheric and Oceanic Sciences and Dep. of Geography, Univ. of
California, Los Angeles, CA 90095. *Corresponding author (h.vereecken@fz-juelich.de).

stress–strain behavior. The spatial distribution of infiltration rates is affected by soil type,
local topography, and attributes of surface cover. Infiltration feedbacks have been shown todrive the formation of vegetation patterns. Infiltration rates are crucial inputs to the designof any irrigation system and many soil and water conservation practices. The significanceof infiltration has made it a subject of studies in many domains ranging from hydrology
to agricultural, environmental, and civil engineering.
Even at the single profile scale, infiltration exhibits strong dynamics that are dependent on soil properties, rainfall characteristics, wetting rates, vegetation cover and type, soil
and crop management, and initial and boundary conditions within the soil flow domain.
Based on the definition of Hillel (1980) and Brutsaert (2005), infiltration is defined as“the entry of water into the soil surface and its subsequent movement through the soil
profile.” The sources of liquid water for infiltration include direct precipitation (rainfall,
dewfall, and snowmelt), leaf drip and stem flow, irrigation, or runoff that was routed overthe land surface and reinfiltrates (a process termed run-on). A detailed understanding of
the primary controls on infiltration rates and the onset of ponding with subsequent runoff,and their translation into model equations, is of great importance at all scales. The accurate process representation of infiltration is also essential for crop water use studies, thedesign of irrigation systems, and the optimal management of water resources. Different
approaches have thus been developed over the past decades to provide quantitative tools
able to describe and predict infiltration into porous media such as soils, ranging fromempirical expressions to analytical and numerical solutions of the basic flow equations.
Infiltration dynamics are determined by soil properties (hydraulic conductivity, sorptivity), the hydraulic gradients that drive flow, and initial and boundary conditions (Philip,
1957). Depending on the initial soil water content in the soil profile, the water supply rate,
and the corresponding soil wetting dynamics, all the available water can infiltrate in different amounts into the soil. Hence, these factors will influence the infiltration curve (thechange with time of the infiltration rate during a wetting event) (Mein and Larson, 1973).
The infiltration capacity or potential infiltration rate of a soil, qcap(t), is the maximum
rate at which the soil surface can take up water for given initial conditions (Horton, 1940).
The actual infiltration capacity is also affected by the initial soil water content of the soil,
but for practical considerations it may be considered a time-dependent soil property wherewater inputs in excess of this maximum infiltration rate will pond and likely generate
runoff. For surface fluxes at rates lower than the soil’s infiltration capacity, the realized
infiltration rate will depend on the state of the soil (as shaped by the temporal history of theapplication rates and the consecutive sequences of wetting and drying). Two infiltration
regimes in unsaturated soils can be distinguished and lead to different occurrence timesof ponding and thus runoff generation: constant water supply (occurring during irrigation or simulated rainfall) and variable-rate supply (during natural rainfall), as is shownin Fig. 1 and 2 in terms of rates (top panel) or cumulative depths (bottom panel; in thiscase uppercase symbols are used).
The dashed curves in Fig. 1 and 2 represent the infiltration capacity, qcap(t) [L T−1]
or Qcap(t) [L], of a given soil profile, the dotted line depicts the water supply rate, r(t)[L T−1] or R(t) [L], for the constant precipitation rate (Fig. 1) or for the variable precipitation rate, as is often the case during natural rainfall (Fig. 2), and the solid line shows theactual infiltration rate, q(t) [L T−1] or Q(t) [L], during these events. In the case of constantr(t), all the applied water can infiltrate in the first stage of wetting (Fig. 1), and q(t) = r(t).
Compared with qcap(t) corresponding to an unlimited water application rate, q(t) can behigher than qcap(t) because the hydraulic gradients resulting from the unsaturated condition generated by the limited wetting rate, r(t), are larger than the ones that result fromthe saturated condition inherent to the unlimited wetting rate. However, both curves
tend toward a similar “quasi-steady” infiltration rate corresponding to a gradient of unity
but approach it at different rates resulting from the difference in the rate of decrease of
the hydraulic gradients corresponding to each wetting condition. As wetting progresses,

Vadose Zone Journal | Advancing Critical Zone ScienceFig. 2. (A) Illustrative representation of the infiltration capacity rate,qcap(t) (dashed curve) and the actual infiltration rate, q(t) (solid line) ofa soil profile exposed to a rainfall event with variable intensity r(t) (dotted line)—for low r(t) values below qcap(t), q(t) = r(t), and the result
is that tc is close to tp; and (B) cumulative infiltration capacity Icap(t)(dashed line), cumulative actual infiltration, I(t) (solid line), and cumulative rainfall R(t) (dotted line)—in this case, tp cannot be estimated bymeans of Icap(tp − t0). Adapted from Assouline et al., 2007.
Fig. 1. (A) Illustrative representation of the infiltration capacity rate,qcap(t) (dashed curve) and the actual infiltration rate, q(t) (solidline) of a soil profile exposed to a rainfall event with constant intensity r(t) (dotted line), where the time tc denotes the moment whereqcap(t) = r(t), while the time tp denotes the ponding time whereq(t) = r(t), and shifting the representation of qcap(t) by t0 = (tp − tc)(dashed-dotted line) allows estimation of tp ; (B) illustration of tp in
terms of the cumulative infiltration, with I(tp) = R(tp); and (C) cumulative infiltration capacity Icap(t) (dashed line) and cumulative rainfall
R(t) (dotted line); the ponding time tp can be estimated by means ofIcap(t − t0) (dashed-dotted line).

q(t) begins to decrease, and at the time of ponding, tp, where theinfiltration rate, q(tp) is smaller than the wetting rate, r(tp), ponding occurs at the soil surface.
In terms of cumulative infiltrating volumes of water (or water
depths, which are volumes per unit area) (Fig. 1B), Q(t) is always
smaller than Qcap(t) even when q(t) is higher than qcap(t) (Fig. 1A),
indicating that the concept of infiltration capacity corresponds tocumulative infiltration rather than to infiltration rates. Because
the rainfall rate is constant, the picture in terms of infiltration

VZJ | Advancing Critical Zone Science

rates (Fig. 1A) and cumulative depths (Fig. 1B and 1C) provide
a similar condition for tp, corresponding to the time when Q(t)departs from R(t). It is interesting to note that if it is assumed thatponding occurs when r(t) exceeds qcap(t), tp is underestimated (Fig.
1A and 1C). The infiltration regime under variable r(t) and theprediction of tp is much more complicated because q(t) depends
on the pattern of r(t) (Fig. 2A). The picture is much simpler when
expressed in terms of cumulative depths (Fig. 2B). Assouline etal. (2007) have proposed a simple method to estimate tp undervariable r(t).
The cases described in Fig. 1 and 2 represent one mechanism of runoff formation termed infiltration excess or Hortonian
overland flow (Horton, 1933). Another important mechanism corresponds to the formation of a saturated topsoil layer that enhances
runoff formation. This mechanism is termed saturation excess orDunne overland flow (Dunne, 1978; Freeze, 1980). Dunne overland flow occurs when the soil reaches saturation from above, via
direct precipitation, or below and no additional water supplied
p. 3 of 53

from the top can be stored in the soil profile. This can happen,
for example, when the local water table is high or a hydraulically
impeding layer close to the soil surface exists (causing a so-called“perched” water table). In general, the Dunne mechanism occurs inareas close to the channels of the drainage network of a catchment
or in areas in the low-lying parts of a catchment where the depth to
the water table is shallow. It is therefore more common in humid
climates. These saturated runoff source areas vary in size, seasonally, and during individual storm events. Therefore, they are often
referred to as variable source areas (Dunne and Black, 1970) andthe runoff generated on them as variable source area runoff (as, e.g.,
in TOPMODEL [Beven and Kirkby, 1979], which is implemented[sometimes in modified form] in many LSMs). On the other hand,
Hortonian overland flow is more common in semiarid climates(Entekhabi and Eagleson, 1989).
Infiltration theory, and related numerical and analytical solutions, was originally the domain of soil physicists. Several reviews
on different aspects of infiltration into soils have been published(Assouline, 2013; Barry et al., 2007; Clothier, 2001; Gardner,
1960; Hopmans et al., 2007; Parlange, 1980; Parlange et al., 1999;
Philip, 1969; Philip and Knight, 1974; Raats, 2001; Raats et al.,
2002; Skaggs, 1982) and also constitute key chapters in textbooks(e.g., Bear, 1972; Brutsaert, 2005; Childs, 1969; Chow et al., 1988;
Delleur, 2006; Hillel, 1998; Warrick, 2003). The fundamental
concepts of infiltration have been applied in hydrology to deal withthe prediction of infiltration at the field scale (Corradini et al.,2011; Govindaraju et al., 2012; Morbidelli et al., 2006), on hillslopes (e.g., Morbidelli et al., 2018), for heterogeneous soil systems(Govindaraju et al., 2001), and to handle the impact of complex
precipitation events and patterns on infiltration (Corradini et al.,
1994, 1997).
With the advent of efforts to model the global water andenergy cycles at the large scale, infiltration theory was taken up by
the climate and hydrological modeling community. Largely due
to the (originally) limited computing power and the difficulty
in defining spatially distributed and appropriately upscaled soil
parameters, simplifications and approximations of the infiltration process were introduced into the LSMs embedded in weather
and climate models. The main role of a LSM is to compute theenergy partitioning at the interface between land surface andatmosphere. At the land surface, net radiation is converted into
latent heat, sensible heat, and ground heat flux, where the latent
heat flux is the equivalent of the evapotranspiration flux in thewater balance but in this case expressed in energy units. The energy
partitioning at the land surface directly affects the state of theatmosphere. For example, the relative magnitude of the latent andsensible heat fluxes will modify atmospheric state variables suchas relative humidity and the height of the atmospheric boundary layer, which in turn will affect cloud-forming processes andultimately rainfall. Infiltration acts on this energy partitioning
indirectly via its control on soil moisture content. Near-surface soil
moisture content is an important state variable in both the waterand energy balance. For example, it affects net radiation due to its

VZJ | Advancing Critical Zone Science

effects on land surface radiative properties, albedo, and emissivity.
Also, both soil evaporation and transpiration depend strongly on
soil moisture content. A reduction in soil moisture content will
lower evaporation via a reduced soil–atmosphere vapor gradientand decreased replenishment of water to the evaporation front
due to reduced hydraulic conductivity, whereas transpiration isaffected via a decrease in root water uptake under drought conditions. Finally, because soil thermal conductivity and heat capacity
are functions of soil moisture, infiltration also indirectly affects
the soil heat regime. Anwar et al. (2018) showed that the choiceof the infiltration scheme had a significant effect on the simulated
regional climate. An infiltration scheme with a lower soil infiltration rate yielded lower topsoil soil moisture, which led to a lower
latent heat flux and a higher sensible heat flux, resulting in a net
warming effect within the simulation domain.
With regard to the simulation of infiltration, a range ofapproaches at different levels of complexity currently exists in thesemodels. Paniconi and Putti (2015) reviewed the last five decades ofphysically based numerical models in hydrology and addressed thetreatment of infiltration from local via hillslope to catchment scale.
They focused on the flow path heterogeneity, where analysis has
been focused on the nonlinearity and upscaling in hydrology witha specific focus on numerical methods used in hydrological models
and related computational challenges. They briefly discussed theseminal work of Horton (1933), Betson and Marius (1969), andDunne and Black (1970) in identifying the main mechanismsof runoff generation, which is closely related to the infiltration
process.
Zhao and Li (2015) reviewed the different approaches tomodel land surface processes across complex terrain. The main
focus was on the role of grid-scale spatial heterogeneity of land
surface variables and parameters (e.g., soil moisture content, net
primary productivity, leaf area index, fraction of vegetation cover)and the topographic impact on key atmospheric controls (e.g.,precipitation, air temperature, wind speed, air pressure). The
role of infiltration and its different parameterization was briefly
addressed in relation to the spatial variability of soil hydraulicproperties in complex terrains.
Clark et al. (2015) analyzed the state-of-the-art of infiltration
processes in LSMs. They concluded that the main challenges arein the appropriate treatment of the small-scale heterogeneity todescribe the large-scale fluxes of infiltration and surface runoff
and the need for an improved description of wetting front dynamics, which may lead to improved simulations of infiltration andsurface runoff during heavy rainfall. More recently, Mueller et al.
(2016) examined the potential of the LSMs SWAP, JULES, and
CH-TESSEL to produce surface runoff for intense rainfall events.
Based on the results, they recommend that future work should
consider a correction of the maximum infiltration rate in JULES
and investigate its numerical scheme to make it suitable for highvertical resolution. Recently, Morbidelli et al. (2018) reviewed therole of slope on infiltration. They pointed out the need to further
develop upscaling approaches up to catchment and subgrid scalesp. 4 of 53

and to establish a theoretical framework to describe infiltration
on hillslopes to better explain experimental observations that havebecome available in recent years.
In this review, we briefly recapitulate the main approaches andsolutions derived from soil physical theory to describe infiltration
processes at the point and field scales and the techniques used tonumerically solve the infiltration processes and finally to quantify
the impact of spatial variability on infiltration. Here, we focus on
infiltration processes in non-frozen soils. We present the infiltration approaches used in various LSMs and address how maximum
infiltration capacity is quantified, how soil moisture and spatial
variability of soil properties are parameterized, and how the areal
saturation fraction, Fsat, important for Hortonian runoff, is estimated. We also identify key soil parameters that affect the soil
infiltration and runoff and present upscaling approaches for soil
hydraulic parameters applicable to the grid scale of LSMs. Finally,
we will present an outlook and future perspectives on modelinginfiltration in LSMs. Table 1 shows the content and organizationof the review, referring to the specific sections.

6 Quantitative Expression of the

Infiltration Process at Point Scale
Basic Physical Models of Soil Water Flowand Infiltration

In the chronological development of the conceptual modeling of infiltration in porous media, there have been three main
approaches.
Darcy (1856) formulated the first quantitative description offlow through a saturated porous medium, known as Darcy’s law:
J =-K s

¶H
¶z

[1]

where J is the Darcian flux of water [L T−1] at time t [T], Ks is a
proportionality constant characterizing the medium and named
the saturated hydraulic conductivity [L T−1], and ¶H/¶z (dimensionless) is the hydraulic gradient calculated from the differencesin total hydraulic head (H) and the (vertical) distance z [L] within
the saturated porous medium (in saturated soils, the total head H
is the sum of pressure [h] and elevation heads [z]).
Buckingham (1907) proposed to extend Darcy’s law tounsaturated water flow, where the actual water content in the
porous medium, q, is lower than its maximum value at saturation,
q s. The main assumption is that the constant saturated hydraulicconductivity value, Ks, could be replaced by a function of soil watercontent, q, or matrix potential, h, as the characteristic of the unsaturated porous medium. That function was named the unsaturated
hydraulic conductivity function and given the symbol K(q) or K(h).
Following the notation of Eq. [1], the resulting unsaturated flow
equation resulting from Buckingham’s assumptions is
J =-K ( h )æ ¶h ö
¶H=-K ( h )ççç +1÷÷÷
è ¶z
ø
¶z

VZJ | Advancing Critical Zone Science[2]

Table 1. Contents.
General Introduction ........................................................................................................................................................1
Quantitative Expression of the Infiltration Process at the Point Scale ...............................................................5
Basic Physical Models of Soil Water Flow and Infiltration...............................................................................5
Empirical Infiltration Equations..............................................................................................................................6
Analytical and Semi-analytical Infiltration Description ..................................................................................6
Numerical Solutions of Profile Scale Infiltration Processes .............................................................................8
Numerical Methods................................................................................................................................................8
Governing Flow Equations ...................................................................................................................................8
Boundary Conditions ............................................................................................................................................8
Vertical and Temporal Discretization ...............................................................................................................9
Upscaling Approaches toward Larger Scales ........................................................................................................... 10
Upscaling Spatially Heterogenous Parameters Relevant for Infiltration ................................................. 10
Topography-based Aggregation ....................................................................................................................... 12
Homogenization................................................................................................................................................... 12
Similarity Scaling ................................................................................................................................................. 12
Bayesian Upscaling .............................................................................................................................................. 12
Machine Learning Based Upscaling ............................................................................................................... 13
Multiscale Parameter Regionalization ........................................................................................................... 13
Other Parameter Scaling Approaches ............................................................................................................ 13
Upscaling Infiltration Processes............................................................................................................................ 13
Similarity Scaling ..................................................................................................................................................14
Stochastic Upscaling ............................................................................................................................................15
Infiltration Processes in Land Surface Models ........................................................................................................16
Grid-Scale Infiltration Processes...........................................................................................................................16
The Maximum Infiltration Rate in Land Surface Models ........................................................................19
Spatial Heterogeneity in Soil Water Content Using Probability Densities ........................................ 22
Representing Spatial Heterogeneity of Surface Properties that Control Infiltration in a Single
Parameter ..................................................................................................................................................................... 23
Estimating the Areal Saturation Fraction ..................................................................................................... 24
Characterizing the Surface Saturated Hydraulic Conductivity ............................................................. 24
Spatial Distribution of the Saturated Hydraulic Conductivity ......................................................... 24
Estimators for the Saturated Hydraulic Conductivity ......................................................................... 25
Use of Soil Maps to Estimate the Saturated Hydraulic Conductivity ............................................ 25
Saturated Hydraulic Conductivity and Infiltration in Organic Soils ............................................. 25
Numerical Treatment of Infiltration in Land Surface Models ..................................................................... 26
Sensitivity of Infiltration–Runoff Process to Model Parameters ...................................................................... 27
Improving the Infiltration Process in Land Surface Models............................................................................... 28
Soil Structure .............................................................................................................................................................. 29
Hysteresis in the Soil Water Retention Curve and Thermal Effects on Hydraulic Properties ........... 29
Soil Water Repellency ............................................................................................................................................... 30
Compaction, Swelling, and Shrinkage ................................................................................................................ 30
Freeze and Thaw..........................................................................................................................................................31
Impermeable Soil Layers ...........................................................................................................................................31
Instability of Different Flow Regimes ................................................................................................................. 32
Solution of Numerical Issues .................................................................................................................................. 32
Summary and Conclusion ............................................................................................................................................. 33
Appendix A1 ..................................................................................................................................................................... 37
Appendix A2 ..................................................................................................................................................................... 37
Appendix A3 ..................................................................................................................................................................... 39
Appendix A4 ..................................................................................................................................................................... 42

where z [L] is the vertical coordinate being positive upward, andz = 0 represents a prescribed reference level. In this case, the totalhead H is the sum of the matrix potential in head units, h, and thegravitational head z.
Finally, Richards (1931) combined the flow equation of
Buckingham (1907) (Eq. [2]) and the principle of continuity
assuming an infinitely mobile air phase in the soil (zero resistanceto air flow). The resulting well-known and widely used one-dimensional expression for vertical water flow isæ ¶ h öù¶q ¶ é
ê K ( h )çç +1÷÷ú
=
ç
ê
è
ø÷úû
¶t ¶z ë¶z

[3]

Equation [3] requires the definition of an additional characteristic of the porous medium, which is the soil water retention
curve, where the actual water content is a function of the matrix
potential, q(h).
The q(h) and K(h) functions represent the hydraulic properties of the porous medium. An up-to-date presentation of themathematical expressions (e.g., Brooks and Corey, 1964; van
Genuchten, 1980) used to quantify these properties in soils canbe found in the review of Assouline and Or (2013).
The Brooks and Corey (1964) retention function is given byp. 5 of 53

æ h ö÷-l p
S e ( h ) = ççç ÷÷ ; h < hc
çè h ÷ø

[4]c

S e ( h ) = 1;h ³ hcwhere hc is the matric potential at air-entry value [L−1], l p is adimensionless pore size distribution index, h is the pressure head[L], and Se is the effective saturation (dimensionless). given by
Se =

q-q r
f-q r

[5]

where Se is the effective saturation (dimensionless), q is the actual
water content [L3 L−3], q r is the residual water content [L3 L−3],
and f is the porosity [L3 L−3], which can be related to q s as thesaturated water content.
The corresponding unsaturated hydraulic conductivity function is
3+2/l p

K = K s Se

[6]

where Ks is the saturated hydraulic conductivity [L T−1].
The water retention function proposed by van Genuchten(1980) is
q s -q rq( h ) = q r +m[7]n
1 + h hc vG()

Empirical Infiltration Equations
The necessity to describe quantitatively such a crucial hydrological process as infiltration combined with the complexity ofthe solution of the flow equation (Eq. [3]) led to the developmentof empirical expressions relating the infiltration rate, q, to time,
t. From Fig. 1 it can be seen that, for a constant water supply, r(t),the infiltration rate q(t) gradually decreases and tends towarda steady final infiltration rate, qf. Consequently, the different
forms of the suggested empirical equations describing infiltration
are all monotonically decreasing functions based on exponential
or power law decays, for which the parameters do not generally
have a physical meaning and are evaluated by fitting to experimental data.
In general, there are two main families of equations used todescribe the infiltration process—one is based on the original
formulation of Horton (1941) and the other on those presented
by Kostiakov (1932). Over time, both equations have been modified and extended by various researchers such as the Kostiakovmodel by Lewis (1937), Mezencev (1948), Smith (1972), Parhi etal. (2007), and Furman et al. (2006) to account for different initialand boundary conditions. The original formulation proposed by
Horton (1941) predicts an exponential decay of the infiltration
flux q over time t asq ( t ) = qf + ( qi - qf )exp (-a K t )[9]where m is a shape factor (dimensionless) often assumed to berelated to nvG by m = 1 − 1/nvG. This function is continuous in hand presents an inflection point, making it more appropriate forapplication in numerical solutions.
Applying the model of Mualem (1976) to Eq. [7] leads to thefollowing unsaturated hydraulic conductivity function:

where q is the infiltration rate [L T−1], qf is the final (constant)infiltration rate [L T−1], qi is the initial infiltration rate [L T−1],
and a K is the decay constant [T−1]. For t ® ¥, qf can be related
to the hydraulic conductivity of the wetted soil layer.
On the other hand, Kostiakov (1932) introduced a power law
equation with two fitting parameters, nK and bK, in the formm ù2éK ( h ) = K s S el ê1-(1- S e1/ m ) úêëúûq ( t ) = qf + nK b K t (b K -1)[8]

where l is a fitting parameter dependent on the soil type
(dimensionless).
The Richards equation (Eq. [3]) represents the actual physicalmodel that can be used to simulate and illustrate the infiltration
process in porous media. It addresses only the macroscale behavior
and is valid for a representative volume for which the prescribedhydraulic properties can be applied. The solution of the Richards
equation requires the definition of initial and boundary conditions.
When water is ponding at the soil surface, infiltration is governed by
the hydraulic head at the soil surface, a concentration type boundary condition (BC) known also as a Dirichlet BC. When the water
application rate is below the soil infiltration capacity, a flux or
Neumann BC can be applied at the soil surface. Because Eq. [3]is highly nonlinear due to the nonlinear character of the hydraulicconductivity function [K(h) or K(q)], analytical solutions can bederived only for specific initial and boundary conditions and withknowledge of the soil hydraulic properties. Consequently, solutions
to practical soil water flow problems generally require the use ofnumerical schemes designed to solve partial differential equations.

VZJ | Advancing Critical Zone Science

[10]

where nK and bK depend on the initial and boundary conditionsof the soil system.

Analytical and Semi-analytical
Infiltration Description
Compared with the empirical models, analytical and semianalytical equations based on the Richards equation were also
developed; Philip (1957) and Philip (1969) presented the first
analytical solution to the Richards equation. It considers infiltration as a sorption process with a perturbation generated by thepresence of gravity. Therefore, this method corresponds by definition to the first stages of infiltration into a relatively dry soil
profile where gravity plays only a minor role. It can therefore be
expected to be applicable for small and intermediate values of time
(Brutsaert, 2005; Hillel, 1998). Because the original solution, in
the form of a power series, diverges for large values of time, Philip(1957) proposed to use a truncated version of the original solution,
considering the first two terms of the series:

1
q ( t ) = S op t - 1/2 + A
2[11]

p. 6 of 53

where Sop is defined as the soil sorptivity [L T−2], t is time [T],and A is a fitting parameter (dimensionless). Note that Eq. [11]is equivalent to Eq. [10] for b K = 1/2, nK = Sop, and qf = A. For
long infiltration times (during ponded infiltration), the term A
approaches the value of Ks. Philip (1969) and Talsma and Parlange(1972) have shown that even for long infiltration times (of practical interest), the inequality ⅓Ks £ A £ ⅔Ks holds. The spatial
vertical water content distribution in the soil profile following
infiltration and wetting is quite complex and can be illustrated
by Fig. 3. Experimental results from Davidson et al. (1963) andRubin and Steinhardt (1964) have shown the sigmoid-like natureof the distribution of q(z) (here, the data of Davidson et al. [1963]
are depicted in Fig. 4).
Green and Ampt (1911) presented an approach that is basedon fundamental physics but makes strong assumptions about thesoil hydraulic properties and the shape of q(z). Their key assumption is the presence of a sharp infiltration front moving downward
with time instead of the sigmoid distribution of water content withdepth representing a more realistic wetting front, as shown in Fig.
3 and 4. The presented infiltration function integrated macroscale
physical entities such as pressure head differences with depth, soil
porosity, and saturated hydraulic conductivity over the wetter soil
layer. The Green–Ampt model, for the case of ponding infiltration
with a negligible ponded water depth at the surface, is represented
by an implicit expression for the infiltration rate, q(t), which needs
to be solved iteratively for the cumulative infiltration Q(t):é h ( q -q ) ù
q ( t ) = K s êê f s i +1úúëê Q ( t )ûúSwartzendruber (1987), Ross et al. (1996), and Barry et al. (1995,
2005). Selker and Assouline (2017) presented a simple explicit
solution for the position of the wetting front in time based on
approximating the term describing early infiltration behavior bymeans of the sum of gravitational flow and the exact solution forcapillary imbibition:q( t ) = K s

Af K s + nf K s hf 2 t
1 + Af ( K s t nf hf ) +( 2 nf hf ) K s t[13]where A f is a fitting parameter that can be approximatedby (2/3). For infiltration into a completely dry soil profile,nf = q s − q r, whereas for a profile with known initial water content q i, nf = q s − q i.
The result of Eq. [13] is within 1% of the exact implicit solution of vertical Green and Ampt infiltration (Eq. [12]). Hence, theproposed approximation adds essentially no error to the Green and
Ampt approach but greatly simplifies the computation of infiltration by allowing an explicit expression that is, in theory, easy to
implement in LSMs.
For near-constant rainfall rates, q(Q) is independent of theapplied rate r (Skaggs, 1982; Smith et al., 2002), and infiltration
capacity at any given time depends only on the cumulative infiltration volume, regardless of the previous rainfall history. The[12]

where Q(t) is the cumulative infiltration [L] at time t [T], hf isthe capillary head at the wetting front [L], Ks is the soil saturated
hydraulic conductivity [L T−1], and q i is the initial water content
of the soil profile. The matrix head hf is related to soil hydraulicproperties (Bouwer, 1964; Neuman, 1976).
Generalized and exact solutions for Q(t) were developed by Parlange et al. (1982, 1985), Haverkamp et al. (1990),

Fig. 3. Illustration of the two-dimensional spatial distribution of soil
moisture content within a wetted soil profile during infiltration (left)
and its corresponding soil moisture profile (right). A–A¢ is the verticaltransect depicted on the right-hand side of the left panel. (From Or etal. class notes with permission.)

VZJ | Advancing Critical Zone ScienceFig. 4. Infiltration into Hesperia sandy loam soil (initially air dry)showing the soil water content with depth and the propagation of thewetting front with time during infiltration. (Reprinted with permission from Davidson et al., 1963.)

p. 7 of 53

time invariance of q(Q) holds true also when a layered profile ora sealed soil surface is considered (Mualem and Assouline, 1996;
Smith 1990). This is an important characteristic of the infiltration
process and the basis of the so-called time compression analysis
for predicting the timing of ponding and runoff (Assouline et al.,2007; Brutsaert, 2005).

resolution of precipitation data is extremely important to capture
the generation of excess water (Hortonian excess). Finally, using
daily accumulated rainfall fluxes will lead to an underestimationof the generated excess water compared with the use of highly
resolved rainfall data (Mertens et al., 2002).

Numerical Solutions of Profile-Scale
Infiltration Processes

Most early applications of numerical methods for solving
variably saturated flow problems usually used the classical finite
differences method (e.g., Bresler, 1973; Bresler and Hanks, 1969;
Hanks et al., 1969; Rubin and Steinhardt, 1964). Integrated finite
difference, finite volumes, and especially finite element methods
became increasingly popular in the ‘70s and after (e.g., Huyakorn
et al., 1986; Neuman, 1973; Paniconi and Putti, 1994; Šimůnek
et al., 2008), accompanying the fast development of computerswith increasing computational speed and memory. While finite
difference methods are still used today in the majority of onedimensional models, finite volume methods and/or finite element
methods coupled with mass lumping of the mass balance term areusually used in two- and three-dimensional models (e.g., Healy,
2008; Pruess, 1991; Šimůnek et al., 2008). An overview of thesedevelopments with respect to infiltration prediction was given in
the review of Assouline (2013).
Furthermore, a number of LSMs use semi-implicit numerical solutions to Richards’ one-dimensional partial differential
equation such as in Best et al. (2011), OLAM-SOIL (Walko
and Avissar, 2008), ORCHIDEE (de Rosnay et al., 2002), andISBA-SURFEX (Boone et al., 2000; Decharme et al., 2011).
Conventional methods for solving this highly nonlinear equation inevitably lead to numerical and accuracy challenges that
impact their hydrological performance.

While the empirical, analytical, and semi-analytical models
discussed above provide descriptions of the infiltration process
for relatively simple conditions (e.g., a homogeneous soil profile,constant initial saturation with depth, constant ponding at thesoil surface), the quantitative analysis of real-world infiltration
problems (e.g., a layered soil profile, variable initial saturation,
time-variable rainfall, limited ponding) can be obtained only
using numerical solutions of the water flow equation [i.e., Richards’
equation (Eq. [3])]. The reason lies in the fact that the highly
nonlinear eliptic parabolic Richards equation cannot be solved
analytically, apart from for a very limited number of cases.
Since the 1960s soil water balance models have been developed that quantify and predict infiltration in soil by numericallysolving the Richards equation for different BCs and resulted ina wide range of software tools such as HYDRUS (Šimůnek et al.,2008), COUP (Jansson, 2012), DAISY (Hansen et al., 2012), andSWAP (van Dam et al., 2008), among many others. By using suchmodels based on the Richards equation along with sufficient vertical and temporal discretization, the infiltration rate and occurrenceof ponding should be a direct outcome of the numerical solution
in the spatiotemporal domain and do not have to be calculated“externally” using any of the introduced empirical or (semi-)analytical models above. Various simulations have been performed using
HYDRUS-1D to demonstrate water infiltration into a one-dimensional soil profile under natural BCs (Šimůnek and van Genuchten,
2008; Šimůnek et al., 2008, 2016). The HYDRUS-1D model was
selected for this purpose because it is one of the most widely usedand verified codes for unsaturated flow and solute transport modeling (Scanlon, 2004), and the full set of simulations can be foundat https://www.pc-progress.com/Downloads/Public_Lib_H1D/
Using_HYDRUS-1D_to_Simulate_Infiltration.pdf.
Depending on the intensity of precipitation, discretizationin both space and time need to be adapted to the rainfall flux,
whereby most current models rely on a predefined vertical discretization and only the time discretization will be changed byan adaptive time stepping routine. In general, the typical size ofthe vertical discretization is often smaller than 1 cm close to thesoil surface to adequately solve the Richards equation for rainfall
events with high intensities. In some models, the upper BC canbe automatically switched from a flux BC to a fixed pressure head
(Dirichlet) BC during evaporation conditions to avoid the development of extreme pressure head gradients close to the soil surface.
This approach stabilizes the numerical solution of the Richards
equation. Additionally, the appropriate choice of temporal

VZJ | Advancing Critical Zone Science

Numerical Methods

Governing Flow Equations
In general, the Richards equation can be formulated, andthus solved numerically, in three different ways: water content (q)based (often also denoted as the diffusivity form), pressure head(h) based, and in terms of a mixed formulation, when both q andh appear simultaneously in the governing equation as shown in Eq.
[3]. The most popularly used vadose zone flow models currently
utilize the mixed formulation of the Richards equation and solve
this equation using the mass-conservative method proposed by
Celia et al. (1990).

Boundary Conditions
Infiltration rates for a point source infiltration process canbe obtained by numerically solving the Richards equation for an
appropriate upper (soil surface) BC. In general, two types of BCs
can be used to simulate the infiltration process, i.e., the Dirichlet
or Neumann BC. The Dirichlet (pressure head based) BC fixes thepressure head, h, at the soil surface (z = 0) to a value h0 [L], which
can be either constant or variable with time:h( z , t ) = h0 ( t ) at z = 0[14]

p. 8 of 53

The value for h0 can be either negative (e.g., for tension disk
infiltration), zero, or positive (ponded infiltration). The majority of empirical, analytical, and semi-analytical models discussed
above represent conditions with zero (or slightly positive) pressure
head at the soil surface. The Neumann (flux-based) BC fixes thewater flux (infiltration), q, at the soil surface to a required water
flux, q 0 [L T−1], which can again be either constant or variablewith time:æ ¶h ö
-K ççç +1÷÷÷ = q0 ( t )
è ¶z
ø

[15]

where q0 is negative for infiltration and positive for evaporation
when the z axis is defined positive upward if z is defined negativewith depth. If water is allowed to build up on the soil surface after
the onset of soil surface ponding, a “surface reservoir” BC may be
applied (Mis, 1982; Šimůnek et al., 2008; van Dam et al., 2008):
æ ¶h ö
dh-K ççç +1÷÷÷ = q0 ( t )è ¶z
ø
dt

[16]

The flux q 0 in Eq. [16] is the net infiltration rate, i.e., thedifference between precipitation and evaporation. Based on Eq.
[16], the height h(0,t) of the surface water layer (ponding height)
increases due to precipitation and reduces because of infiltrationand evaporation of the ponding layer.
Kollet and Maxwell (2006) closed the problem of variably
saturated groundwater flow, infiltration, and surface water flow byapplying flux and pressure continuity conditions at the top boundary, leading to a free-surface overland flow BC:æ ¶ h ö ¶ h,0-K ççç +1÷÷÷ =-Ñ v h,0 + q0 ( t )
è ¶z
ø
¶t

[17]

where v is the depth-averaged velocity vector [L T−1], which canbe expressed in terms of Manning’s equation (e.g., in Chow et al.,1988) to establish a flow depth–discharge relationship.
The Dirichlet (Eq. [14]) and Neuman (Eq. [15]) BCs are
system-independent BCs for which prescribed quantities (i.e.,pressure heads or water fluxes) do not depend on the conditionsof the soil profile, its saturation status, or its infiltration capacity.
These BCs thus may not properly describe real-world conditions
in which infiltration or actual soil evaporation rates depend on
the conditions of the soil profile and its saturation status, which
may limit infiltration or evaporation. In many applications, neither
the flux across nor the pressure head at a boundary is known apriori but follows from interactions between the vadose zone andits surroundings (e.g., the atmosphere). External meteorologicalconditions thus control only the potential water flux across thesoil surface, while the actual flux also depends on the prevailing(transient) soil moisture conditions near the surface. This occurs,for example, when the precipitation rate exceeds the infiltration
capacity of the soil, resulting in the accumulation of excess wateron top of the soil surface and surface runoff, depending on soil
properties and on topographic conditions (for two- and threedimensional representation). Subsequently, the infiltration rate

VZJ | Advancing Critical Zone Science

is no longer controlled by the precipitation rate but instead by thesoil infiltration capacity.
Such conditions may be best described using system-dependent BCs, which take soil moisture conditions and the hydraulicconductivity of the soil near the soil surface into consideration.
For example, in HYDRUS-1D, such a system-dependent BC iscalled an “atmospheric” BC. For these conditions, the soil surface
BC may change from a prescribed flux to a prescribed head typecondition (and vice versa). The numerical solution of Eq. [3] is thenobtained by limiting the absolute value of the surface flux by thefollowing two conditions (Neuman et al., 1974):-K¶h- K £ ip
¶z

[18]

and
hA £ h £ hS

[19]

where ip is the maximum potential rate of infiltration or evaporation under the current atmospheric conditions [L T−1], andhA and hS are, respectively, minimum (i.e., for evaporation) andmaximum (i.e., for infiltration) pressure heads at the soil surface
allowed under the prevailing soil conditions [L]. When one of theendpoints of Eq. [19] is reached, a prescribed head BC will be used
to calculate the actual surface infiltration or evaporation flux.
The value for hA is determined from the equilibrium conditions between soil water and atmospheric water vapor (e.g., Feddes
et al., 1974). The value of hS is usually set equal to zero, which
represents conditions when any excess water on the soil surface
is immediately removed via runoff once ponding is reached. In
this case, the Neumann BC (Eq. [15]) is switched internally inthe model to the Dirichlet BC (Eq. [14]) with h = 0 once ponding
is reached and then back to the Neumann BC once rainfall stops(during redistribution) and the pressure head decreases below zero.
When hS is allowed to be positive, it then represents a layer of water,which can form on top of the soil surface during heavy rains before
the initiation of runoff. In such a case, the Neumann BC (Eq. [15])needs to be switched to the Dirichlet BC (Eq. [14]) similarly as
above, but the surface pressure head value h is calculated using thesurface reservoir BC. Once rainfall stops, infiltration (calculated
for a Dirichlet BC) continues until all water from the accumulated
water layer has infiltrated, when the BC is switched back to the
Neumann BC.

Vertical and Temporal Discretization
The numerical solution of the highly nonlinear Richards equation requires relatively fine spatial (on the order of centimeters)(Vogel and Ippisch, 2008) and temporal (on the order of minutes) discretization. Optimal spatial and temporal discretization
depends strongly on the intensity of precipitation–evaporation–
infiltration and the nonlinearity of the soil hydraulic properties,
as well as on numerical stability criteria, involving hydraulic orthermal diffusivity (see, e.g., Best et al., 2005). Simulations withhigh flux rates and strong nonlinearity of soil hydraulic properties
p. 9 of 53

require finer discretization in both space and time. Most current
vadose zone, including land surface, models rely on a predefinedvertical discretization, which is constant in time, while only thetime discretization changes during simulations by an adaptive time
stepping routine.
Spatial discretization should be made relatively small at locations where large hydraulic gradients are expected. Such a regionis usually located close to the soil surface, where highly variable
meteorological factors can cause rapid changes in soil water contents and corresponding pressure heads. Hence, it is generally
recommended to use relatively thin soil layers (small discretization)near the soil surface and then to gradually increase their thickness
with depth to reflect much slower changes in pressure heads atdeeper depths.

6 Upscaling Approaches

toward Larger Scales

Infiltration considered at larger scales is heavily affected by
the heterogeneity of the soils and the land surface. Assouline andMualem (2002, 2006) explicitly demonstrated the impact of heterogeneity, introduced by a combined effect of impervious areasand spatial variability in soil properties, on infiltration. Figure 5
shows the results in terms of the infiltration rate [q(t)] curves fora homogeneous and a heterogeneous field exposed to a constantrainfall rate.
Accounting for field spatial variability leads to shorter
ponding times and to a more gradual decrease in the infiltration flux with time (Smith and Hebbert, 1979). Consequently,
surface runoff will appear earlier in heterogeneous fields than inhomogeneous ones. This results from the fact that part of the heterogeneous field has much lower hydraulic conductivity than thatof the homogeneous one, and this generates the early runoff. As
LSMs typically incorporate infiltration at grid sizes from meters
to hundreds of kilometers, scaling approaches are necessary to
account for subgrid heterogeneity at the model resolution. Also,
parameters relevant for infiltration, as well as experiments forFig. 5. Mean infiltration in a heterogeneous (solid line) vs. a homogeneous (dashed line) field.

VZJ | Advancing Critical Zone Science

direct infiltration measurements, are typically observed and performed at the point scale (see above). Therefore, scaling approaches
are needed to translate the measurements to adequately addressinfiltration characteristics at larger scales. Two general strategies
are discussed in the following: (i) the scaling of infiltration-related
properties and (ii) the scaling of infiltration fluxes themselves. We
have chosen this categorization rather than a scale-based categorization because several of the methods proposed below, such as,e.g., similarity scaling, aggregation, and Bayesian upscaling, can beapplied at a range of scales. Figure 6 shows the different upscaling
approaches to obtain effective parameters at the scale of LSMs.
We distinguish four different categories: (i) the LSM upscaling
approaches that either assign uniform soil properties to a dominantsoil class or use probability density functions (pdfs) of parametersthat reflect subgrid variability, (ii) parameter upscaling methods,(iii) similarity upscaling methods, and (iv) stochastic upscaling
methods. The last three approaches have been mainly developedfor upscaling from the field to catchment scale, while the first
involves downscaling the parameters from grid to point scale
before upscaling the resulting infiltration from point to grid scale.
Within the LSM upscaling approaches, three main methods canbe distinguished: (i) uniform upscaling assuming the soil hydraulicparameters are constant, (ii) empirical upscaling that use pdfs todefine Fsat and the maximum infiltration rate, Imax, that are thenfurther used to calculate grid-scale infiltration, and (iii) physicalupscaling in which infiltration is calculated at the point scale andpdfs are used to upscale it to the grid scale.

Upscaling Spatially Heterogenous Parameters
Relevant for Infiltration
Land surface models need input parameters at the grid scale,i.e., estimates of effective parameters to quantify hydrologicaland energy balance fluxes and to generate soil infiltration-related
properties such as water storage at the grid cell level (see above).
To adequately represent nonlinear relationships between model
parameters and states and fluxes, it is generally agreed that spatial scaling ideally should take place after the model has been run.
However, this is also a matter of computational resources, which
are limited and generally do not afford simulations at the support
scale, which is the area or volume over which a measurement ismade or a state variable defined. Land surface models driven for
global or continental applications typically require upscaling of
information that is available at higher resolution to scales of a few
to tens of kilometers. Here, a simplification of the landscape heterogeneity by dominant class selection (e.g., USDA soil classes) or
simple parameter averaging typically does not account for nonlinear relationships of subgrid processes and may introduce importantbiases on specific LSM variables, including infiltration and runoff(Boone and Wetzel, 1999).
Another method proposed by Noilhan and Lacarrère (1995)
consists in computing grid-scale land surface parameters according
to observed soil textures (sand and clay) aggregated from the highresolution and using a continuous relationship derived from thep. 10 of 53

Fig. 6. Schematic overview of different upscaling methods of soil parameters (original land surface model [LSM] upscaling approach, parameter upscaling, similarity scaling, and stochastic upscaling) described in Infiltration Processes in Land Surface Models (A), Upscaling Spatially Heterogeneous
Parameters Relevant for Infiltration (B), Similarity Scaling (C), and Stochastic Upscaling (D). Differences in the infiltration model used are indicatedby different colors.

textural classification of Clapp and Hornberger (1978) or Cosby
et al. (1984). The scaled effective parameters are often dependent
on the spatiotemporal patterns of the unsaturated system at thesmaller scale and are not simply a function of the average parametervalues obtained from measurements.
General focus has been placed on upscaling soil hydraulicparameters because significant spatial variability of these properties was reported earlier by Nielsen et al. (1973), Warrick andNielsen (1980), Peck (1983), and Logsdon and Jaynes (1996). In
most of these cases, the field saturated hydraulic conductivity, Ks,
was lognormally distributed (Reynolds and Eldrick, 1985; Russo
et al., 1997; White and Sully, 1992). However, steady-state infiltration fluxes and soil surface water content are distributed either
normally or lognormally (Cosh et al., 2004; Kutilek et al., 1993;
Loague and Gander, 1990; Sisson and Wierenga, 1981; Vieira
et al., 1981). Several studies have dealt with modeling the effect

VZJ | Advancing Critical Zone Scienceof the spatial variability of soil hydraulic properties on infiltration (Assouline and Mualem, 2002; Dagan and Bresler, 1983;
Govindaraju et al., 2006; Milly and Eagleson, 1988; Sivapalan
and Wood, 1986; Smith and Hebbert, 1979; Warrick and Nielsen,
1980; Woolhiser et al., 1996).
While heterogeneity in soil hydrology properties, parameters,
and boundary conditions is ubiquitous, understanding of the sensitivity to different aggregation or upscaling methods is limited.
Zhu and Mohanty (2002a, 2002b, 2003, 2004, 2006), Zhu et al.
(2006, 2004), and Mohanty and Zhu (2007) investigated, in aseries of studies, the suitability of various soil hydraulic parameter
upscaling schemes by matching their prediction performances withensemble vadose zone fluxes (infiltration and evaporation) under
different hydroclimatic scenarios for horizontally and verticallyheterogeneous soil systems. Their synthetic experimental results
showed that soil texture, geological layering, groundwater depth,
p. 11 of 53

surface and profile soil moisture status, vertical flux direction(infiltration vs. evaporation), hydraulic parameter statistics (correlations and spatial structures), root distribution in the soil profile,and topographic features and arrangements conjointly determine the “upscaled” pixel-scale soil hydraulic parameters for theequivalent homogeneous medium that delivers the same amountof flux (infiltration or evaporation) as the natural heterogeneous
medium. Thus, different homogenization algorithms (rules) for
different hydrologic scenarios and land attribute complexities weresuggested for parameter upscaling. Several approaches have been
proposed in the literature to upscale soil hydraulic properties fromthe point to larger scales to obtain effective properties (Vereecken
et al., 2007), which are discussed in the following.

Topography-Based Aggregation
Expanding the power average operator of Yager (2001), Janaand Mohanty (2012a, 2012b, 2012c) coarsened the soil hydraulicparameters to the model grid scale in an attempt to study remotesensing-based soil moisture distribution at the watershed scale(Little Washita). Two types of aggregating methods were combinedin this topography-based aggregation technique. By combining thefeatures of both mode-like and mean-type aggregating methods,
the power average technique provided an ideal tool in scaling ofsoil hydraulic parameters for soil pedons. Power averaging operating across a number of spatial nodes uses a support function basedon linear distances in different Cartesian coordinates, allowing
data clustered around a particular value to combine nonlinearly
while being aggregated. Generally, soil pedons clustered around alocation tend to have similar properties; that correlation decreases
as the distance between two points increases. In other words, theaggregating method considers the mutual support the pedons
extend to each other when clustered.

Homogenization
Homogenization by spatially averaging the soil hydraulicparameters is a simple way of upscaling. Numerical studies suchas those presented by Zhu and Mohanty (2002b) and Mohanty
and Zhu (2007) and observational studies such as those of Janaand Mohanty (2012b) examined the impact of the areal heterogeneity of soil hydraulic parameters at the model grid scale on theensemble response of hydrologic fluxes. In particular, arithmetic,
geometric, and harmonic averages of the soil hydraulic parameters were tested with different parameter correlation structures.
Findings of these studies suggested that different averaging
schemes should be used for different soil hydraulic conductivity
parameters and functional forms.

Similarity Scaling
Another method for scaling is based on the theoreticalframework of geometric similarity for porous media introducedby Miller and Miller (1956). The basic idea is that porous media
that are geometrically similar in their microscopic arrangementof particles differ only in terms of their characteristic length scale

VZJ | Advancing Critical Zone Science

lSC . The underlying assumption is that the hydraulic behavior ofa porous medium can be transformed into the behavior of a reference medium by scaling (Nielsen et al., 1998). This scaling factor,
a SC , is defined as the ratio between the characteristic length lSC
of a geometrically similar soil and the characteristic length of thereference soil l r:a SC =

l SC
lr

[20]

This approach can also be used for spatial scaling of soil
hydraulic properties (Mohanty, 1999; Shouse and Mohanty,1998). Zhu and Mohanty (2006) used Miller–Miller scaling in
combination with a one-dimensional infiltration equation basedon the Haverkamp et al. (1990) model to derive effective hydraulicparameters at the grid scale of LSMs and global climate models(GCMs). The heterogeneity of soil hydraulic parameters was considered in the horizontal dimension, and the infiltration process
was described using the concept of parallel stream tubes without
lateral interaction. They found that the variability in Miller–
Miller scaling factors had a stronger influence on the effective gridcell infiltration than the saturated water content and the ponding
depth, which reflects microtopography. This approach was also
implemented by Montzka et al. (2017) to generate an optimized
soil hydraulic properties database at a coarse global grid resolutionof 0.25° from the soil texture information system SoilGrids1km
of Hengl et al. (2014). Further information about similarity andMiller–Miller scaling in terms of scaling infiltration processes isgiven below.

Bayesian Upscaling
Considered as a calibration method at the remote sensing
footprint scale, the Markov chain Monte Carlo based upscaling
algorithm introduced by Das et al. (2008) provides an alternative
to derive the upscaled effective soil hydraulic parameters from atime series of soil moisture observations and stochastic information of the fine-scale soil hydraulic parameter variability. The
Bayesian framework produces probability distributions of theeffective (pixel-scale) soil hydraulic parameters, where preexistingknowledge about the local-scale soil parameters (e.g., from the Soil
Survey Geographic [SSURGO] database) can be combined withdynamic hydrologic observations and model outputs.
Kim and Mohanty (2017) proposed a more general approachby accounting for the effects of mixed (weighted) physical controls(covariates) as well as the interactions between the controls on thesoil moisture distribution and subsurface flow including lateralflow between grid cells. They used a Bayesian averaging scheme
for effectively estimating the contributing ratios (weights) for thephysical controls and their interactions. This scheme provides
proper weights that show how the controls contribute to describing the spatial variability of soil moisture and thus effective soil
hydraulic properties. This approach underpins the concept ofhydrologic connectivity based on the probability of connected
local pathways of surface and subsurface flow, by which emergent
p. 12 of 53

catchment-scale behavior is depicted. The physical control based
connectivity index approach can be easily adopted in land surfacehydrologic and biogeochemical models, leading to an Earth Systemmodeling framework, where geophysical attributes such as geology,
ecotones, and topography are (or will be) the primary drivers forwater, carbon, and energy cycles.

Machine Learning Based Upscaling
Pedotransfer functions (PTFs) have been used as inexpensivealternatives for estimating soil infiltration and hydraulic properties
using soil textural and bulk density information at the local scale.
Expanding the concept of the PTF to estimate landscape-scale soil
hydraulic properties, Sharma et al. (2006) developed artificial neuralnetworks based pedo-topo-vegetation transfer functions with goodsuccess by adding topographic and vegetation characteristics to soil
textural and bulk density information. Advancing the idea within a
Bayesian framework, Jana et al. (2008) and Jana and Mohanty (2011)
developed multiscale Bayesian neural network (BNN) based PTFs.
Using a BNN, they upscaled or downscaled soil hydraulic parameters using soil texture and structure data at one scale to simulate thekey soil moisture contents related to soil water retention at another.
In that study, training inputs to the BNN consisted of the percentage of sand, silt, and clay, the bulk density of the soil, and digital
elevation map and vegetation data, while the targets were the soil
water content at 0 (saturation), 33 (field capacity), and 1500 (wilting
point) kPa. Using Markov chain Monte Carlo techniques, the BNN
provides a distribution of the output parameter instead of a single
deterministic value. Recently, Montzka et al. (2017) used PTFs toderive the effective spatial distribution of scaling factors and themean soil hydraulic properties that can be used in LSMs to quantify
the effect of spatial variability on infiltration fluxes.

Multiscale Parameter Regionalization
The multiscale parameter regionalization method proposedby Samaniego et al. (2010) is a two-step approach with initialregionalization and subsequent upscaling. The regionalization
establishes a priori relationships between model parameters anddistributed basin predictors at the fine scale, leading to linear or
nonlinear transfer functions. These functions are used as global
parameters to reduce overparameterization and ease transferability.
Soil texture and land cover can be used as predictors for infiltration (Samaniego et al., 2010). The upscaling toward coarser scales
is performed by a selection or combination of upscaling operatorssuch as arithmetic mean, maximum difference, geometric mean,
harmonic mean, and majority. With this approach, the subgrid
variability is used for facilitating transferability toward ungauged
and uncalibrated regions for improved model parameterization.
Samaniego et al. (2017) analyzed the applicability of the multiscaleparameter regionalization approach for several LSMs.

Other Parameter Scaling Approaches
Other parameter scaling methods are fractal approaches (e.g.,
Tyler and Wheatcraft, 1990) and the scaleway approach (e.g., Vogel

VZJ | Advancing Critical Zone Scienceand Roth, 1998). Although they have been used over small watersheds, they have strong limitations (Meng et al., 2006). Stochastic
upscaling or aggregation methods statistically allow betteraccounting for landscape heterogeneity conditions. Stochastic
upscaling makes use of geostatistical descriptors of spatial variables
and quantify a probability distribution for each state variable atlarger scales by stochastic perturbation rather than a deterministic
quantity (see below).
Govindaraju et al. (2006) suggested a semi-analytical model tocompute the space-averaged infiltration at the hillslope scale when
spatial variability in both soil properties and rainfall intensity areaccounted for. The soil spatial heterogeneity was characterized bya lognormal distribution of the saturated hydraulic conductivity,
while the rainfall spatial heterogeneity was simulated by a uniform distribution between two extreme rainfall intensities (lowand high). At each location, the soil Ks and the rainfall intensity
was assumed to remain constant during the rainfall event. The
main finding was that ponding time decreases with the increasein the coefficient of variation of Ks (see also Fig. 5).
Other researchers also developed and evaluated methods
that are potentially applicable for transferring infiltration
parameterizations across scales. Hailegeorgis et al. (2015) evaluated four regionalization methods for continuous streamflow
simulation. Their regional calibration method uses the maximum weighted average of a performance measure such as theNash–Sutcliffe efficiency to identify specific parameter sets
per subpixel by the DREAM algorithm (Vrugt et al., 2009).
Another method uses the regional median of each parameter forregionalization, where a limitation is that the correlation structure between parameters is lost. Scaling by the nearest neighborapproach includes the assumption that spatial proximity interms of Euclidian distance explains parameter similarity. The
physical similarity approach assumes that similarity of regionsin physical attributes such as land use, cumulative distributionfunctions of terrain slopes, and/or soil types can explain their
response in the hydrological variable.
The potential of most of these upscaling methods has not
yet been tested in currently used LSMs. As highly resolved soil
information becomes available that allows parameterization ofsubgrid variability, several of these methods may prove to be valuable in estimating effective soil properties controlling infiltration
processes at the pixel scale of LSMs. First results using, e.g., themultiscale parameterization method and similarity approach show
great potential in parameterizing key hydraulic properties at thegrid scale of a LSM.

Upscaling Infiltration Processes
To upscale infiltration processes themselves rather thanproviding a method for the averaging of infiltration-relatedparameters as discussed above, besides scaling the Richards equation for infiltration (Sadeghi et al., 2012; Warrick and Hussen,1993), the Miller–Miller similarity can be applied. The stochastic
upscaling of infiltration is based on the frequency distribution
p. 13 of 53

of infiltration-related parameters and allows the computation ofeffective infiltration fluxes using semi-analytical approximations.

Similarity Scaling
Besides providing a method for the averaging of infiltrationrelated parameters as presented above, similarity scaling can also
be applied to reduce the number of parameters in the analytical orsemi-analytical infiltration equations, which in turn also reduces
the data requirement. Within the similarity scaling approach for
infiltration, the infiltration data consists of a set of scaling factors, one for each location or grid cell, and an average infiltration
function. In the following, mainly two examples of similarity scaling for infiltration are presented. The first example, Sharma et al.
(1980), is based on the similarity scaling of cumulative infiltration
Q [L], which is estimated over time t (Philip and de Vries, 1957):
Q ( t ) = S op, i t 1/2 + ACi t[21]where Sop,i and ACi are the parameters fitted to each of the ith
infiltration measurements within a region. Sharma et al. (1980)
used two alternative approaches to realize the scaling of cumulative infiltration: either (i) by determination of two separate scaling
factors for the sorptivity and steady-state infiltration, or (ii) bydetermination of a single optimized scaling factor. For the first
approach, the two scaling factors a Sop and a AC (dimensionless),for sorptivity Sop [L T−0.5] and steady-state infiltration AC [L T−1],
respectively, were calculated according toæSö2
çç op, i ÷÷[22]÷÷a Sop,i = çççè S op ÷ø÷andæ AC ö÷1/2i ÷
a AC, i =ççççè AC ÷ø÷[23]

where áSopñ and áACñ represent the respective mean values of
those fitted infiltration parameters. The scaling of the infiltration data was subsequently performed by
Q * = aQ

[24]

and

t* =a 3 t[25]

where Q* and t* identify the scaled cumulative infiltration [L] andtime [T], respectively, and the scaling factor a is either a Sop ora AC. At this point, both scaling factors a opt and a AC can be used
to scale the cumulative infiltration data. However, both attempts,
using either a opt or a AC, fail to arrange the scaled infiltration datanear the reference cumulative infiltration curve áQñ defined byáSopñ and áACñ. For the second scaling approach, a single scaling
factor a opt was derived by least-squares fitting (using Eq. [24] and[25]). The sum of squares, SS, between scaled cumulative infiltration and reference cumulative infiltration, given by(

SS = å t * Q * t i * - Q t
i

)

2

i

VZJ | Advancing Critical Zone Science

[26]

was minimized by adjusting a opt based on the cumulative infiltration measurements available at each time step ti for each of themultiple measurements. Finally, with this second scaling approach,all the measurements reasonably coalesced about the reference
cumulative infiltration curve.
By using Eq. [26], Sharma et al. (1980) departed from thegeometric similarity proposed by Miller and Miller (1956), which
would require that a Sop and a AC are identical. A scaling methodology following the second approach (Eq. [26]) is empiricalin nature and usually referred to as functional normalization(Tillotson and Nielsen, 1984). Nonetheless, it is still relevant
because it provides a mathematical framework for the handlingof Philips’ infiltration in terms of spatial variability. Equations[24–26] can be used to account for the variability in infiltration
when only the variance of the scaling factors and the mean infiltration curve are known. Furthermore, this approach provides anopportunity to lump both infiltration parameters ACi and Sop,iinto a single scaling parameter a SC,i for each soil. Yet, to the bestof our knowledge, such an approach has not been implemented in
any land surface or hydrological model. However, for the spatial
scales of LSMs, it remains unclear how appropriate scaling parameters could be derived and implemented in practice.
The frequency distribution of a opt determined by Sharma
et al. (1980) had a mean of about 1, as expected, and a standard
deviation of 0.58. The spatial distribution of a opt did not show
a distinct pattern, which was explained by Loague and Gander(1990) later on, as they determined a spatial autocorrelation lengthof <20 m for infiltration in the catchment under consideration.
This, in turn, would support a purely stochastic approach to thevariability of infiltration. However, this conclusion is based on theassumption of the occurrence of vertical fluxes only.
Hopmans (1989) simplified the scaling of Philips’ infiltrationsuggested by Sharma et al. (1980) by substituting Eq. [22] and [23]into Eq. [21]. Accordingly, the scaled instantaneous infiltration q*[L T−1] yields- 1/2
2
q* ( t ) = a1/2
+a SC
AC
SC S op t[27]

where q* is the scaled infiltration rate [L T−1], and áSopñ [L T-0.5]
and áACñ [L T−1] represent the parameters of the reference infiltration curve, now determined by nonlinear regression to all
infiltration measurements available. Hopmans (1989) also performed functional normalization according to Tillotson andNielsen (1984) with the modified empirical three-parameter
infiltration equation of Kostiakov (see Eq. [10]). Both scaling
approaches were successful and Hopmans (1989) stressed thatthe attractiveness of this lies in the potential to lump the spatial variability of infiltration into a single parameter, which, forexample, allows correlation with other environmental variables.
Furthermore, it potentially removes any intercorrelations between
the infiltration parameters.
Haverkamp et al. (1998) provided inspectional analyses ofMiller–Miller similarity scaling for infiltration according to thetheories provided by Green and Ampt (1911), Philip (1957), andp. 14 of 53

Talsma and Parlange (1972). However, the analyses presented by
Haverkamp et al. (1998) remain at a theoretical level and were nottested in any practical application.

Stochastic Upscaling
A stochastic approach to upscale infiltration can be implemented via Monte Carlo simulations for any type of infiltration
equation and any given distribution of relevant parameters. Thishas been demonstrated, e.g., by Smith and Hebbert (1979) and
Sharma and Seely (1979). A sufficient number of random samples needs to be drawn from a given distribution of infiltration
parameters and, after computation of every Monte Carlo model
realization, the average and/or standard deviation of predictedinfiltration can be obtained. However, due to the assumption thatthe infiltration process is ergodic, large sample sizes are usuallyrequired, resulting in a rather large computational demand, which
in reality prohibits any practical application in large-scale studies
involving hydrological models or LSMs. Notwithstanding, manysemi-analytical solutions to infiltration problems are still validated
against Monte Carlo results (e.g., Corradini et al., 2002; Craig etal., 2010; Govindaraju et al., 2001).
Starting with the infiltration theory according to Philip(1957), while making some considerable approximations, Sivapalan
and Wood (1986) provided semi-analytical solutions to calculate
average infiltration across a grid box for two cases. For the firstcase, a spatially constant rainfall rate and a lognormal distribution of saturated hydraulic conductivity Ks were assumed. In thesecond case, a spatially variable rainfall rate and spatially uniform
distribution of Ks were considered. In both cases, run-on was notaccounted for and the soil profiles were assumed to be verticallyhomogeneous. The approximations derived by Sivapalan andWood (1986) were tested against Monte Carlo simulations, but an
in-depth analysis of the uncertainties was not presented. Moreover,due to the mathematical complexity of the approximations, practical application was not feasible.
Smith and Goodrich (2000) simulated the ensemble behavior
of infiltration according to Parlange et al. (1982) by making useof Latin hypercube sampling. They assumed that the lognormal
probability density functions of Ks and the capillary parameter G
were divided into n equal sub-areas, and infiltration was computed
using the average of each sub-area. The areal average infiltration
was then given as the equally weighted sum of all subdivisions.
This method does not involve any analytical solution and has
the advantage of a significantly reduced computational demandcompared with the Monte Carlo technique. Govindaraju et al.
(2001) assumed spatial autocorrelation of lognormally distributed
saturated hydraulic conductivity and provided semi-analyticalsolutions to the Green–Ampt equation based on series expansion.
As a follow-up, the semi-analytical solutions of Govindaraju et al.
(2001) were combined by Corradini et al. (2002) with a kinematicwave expression to account for run-on effects.
Craig et al. (2010) provided approximations for the estimation of the average infiltration rate áqñ [L T−1] according to the

VZJ | Advancing Critical Zone ScienceGreen–Ampt approach as a function of the precipitation rate Jw[L T−1] and mean saturated hydraulic conductivity, áKsñ, and standard deviation, sKs, of a standard lognormal distribution, f k(Ks), of
Ks. They estimated a dimensionless error term e for Jw, Ks, and X,
where X represents a dimensionless time parameter computed fromtime t, Jw, and a wf, which is defined as the product of the absolutevalue of the wetting front suction head and the difference between
saturated water content and initial water content at the beginningof the infiltration event. Finally, the averaged infiltration rate, q,
is approximated byé ln ( J w X )- K s ùJwú
erfc êê
ú2s ks 2
êë
úû

q=

æ
és
ln ( J w X )- K s úù
s2 ÷ö
1expççç K s + ks ÷÷erfc êê ks ú÷2X
2 øs ks 2
èç
ëê 2
ûú

+

+ Jw òX (t )

0[28]

e ( X ( t ), ks ) f k ( ks )d ks

According to Craig et al. (2010), a comparison between thesemi-analytical upscaling approach given in Eq. [28] with Monte
Carlo simulations revealed errors of <3%, which appears to be
acceptable for any practical purpose. They also provided approximations of infiltration for a given normal distribution of aggregate
parameter a wf. From Monte Carlo simulations, they concluded
that the heterogeneity in Ks is of much greater relevance than theheterogeneity in a wf.
Craig et al. (2010) upscaled the Green–Ampt equation to abasin or (LSM) grid cell. Thereby, the upscaled Green–Ampt equation is based on the information of the distributions of Ks and/or
the initial soil water deficit of the basin or grid cell. Additionally,
preferential flow, which is defined as the uneven and rapid
movement of water through the soil, characterized by regions ofenhanced flux, e.g., wormholes, root holes, and cracks, can also be
incorporated by considering a bimodal distribution of Ks as
q=ò

¥

0

¥

ò0 éë X ( t , a GA ), K s ùû f ka ( K s , a GA )d a GA d K s
GA

[29]

with
X=

1
1 +1 ( J w a GA ) t[30]

with a GA = |hf |(q s − q i), where hf is the pressure head at the wetting front [L], q s is the saturated water content, and q i is the watercontent at the wetting front [L3 L−3]. The ponding time can becalculated as

tp =a GA K s
Jw ( Jw -Ks )[31]

This approach also allows calculating of the saturated land
surface fraction in a basin or LSM grid cell. To our knowledge,
none of the proposed methods related to stochastic upscaling havebeen tested or used in full LSMs or hydrological models.
Finally, Choi et al. (2007) developed a three-dimensional volume-averaged soil moisture transport model based on the Richards
p. 15 of 53

equation to account for local-scale variation of topographic attributes such as elevation, slope, and curvature on subsurface soil
moisture fluxes and thus also the infiltration flux. In this approach,the coordinate system xi Î {x, y, z} is replaced by xi* Î {x, y, z*},where z* follows the terrain elevation, defined as[32]z * = Zg - z

where Zg is the ground surface elevation [L]. The local terrain surface slopes, Sxl*, are defined asæ ¶ Zg ¶ Zg
ö-1÷÷÷S x l * = ççç÷øèç ¶ x ¶ y[33]Second-order Taylor expansion of the hydraulic conductivity
and diffusivity with respect to soil moisture were formulated andsubstituted in a diffusion, pressure-based Richards equation, andthe soil water content and terrain slope at any point were approximated by a first-order perturbation approach. Scale-dependentfunctions are used to characterize subgrid variability, incorporating statistical properties that relate the dependence of soil moisturevariability to terrain attributes. The covariance between soil watercontent and terrain slopes is defined asq ¢ S x l ¢ = rsq s S xl

[34]

where r is the correlation coefficient between the terrain slope andthe soil water content. The correlation coefficient was expressedby Choi et al. (2007) as
2

æ q ö
æ q ö
r = g1 + g 2 ççç ÷÷÷ + g 3 ççç ÷÷÷ 		çè q s ÷øçè q s ÷ø[35]

where g1, g 2 , and g 3 are fitting parameters estimated from observations. Note that other approaches are available to estimate thevariance of the soil water content (e.g., Choi et al., 2007; Qu et al.,2015; Vereecken et al., 2007).
Choi et al. (2007) combined this modeling approach with a
one-dimensional diffusion wave model for surface overland flow,
called it the conjunctive surface–subsurface flow model (CSSF),
and proposed it for use in mesoscale climate simulations. The
CSSF was tested through implementation in the common landmodel (CLM) (Oleson et al., 2008) using a set of offline simulations for catchment-scale basins around the Ohio Valley region.
They showed that CSSF simulated a strong contribution of theeffects of surface or overland flow depths of water on infiltration
estimates that led to increased baseflow generation. In addition,a better representation of the surface–subsurface flow processes
improved the representation of soil moisture spatial variability andmay improve the partitioning of energy at the land surface.

66Infiltration Processes

in Land Surface Models

In the following, we provide an overview of how widely used
LSMs represent infiltration processes at the global scale. Therefore,
we selected 12 LSMs, namely: (i) OLAM-SOIL (Walko et al.,VZJ | Advancing Critical Zone Science2000), (ii) ParFlow-CLM (Kollet and Maxwell, 2006, 2008a),(iii) ORCHIDEE rev4783 (d’Orgeval et al., 2008; Ducharne
et al., 2017), (iv) Catchment Land Surface Model (CLSM) (DeLannoy et al., 2014; Ducharne et al., 2000; Koster et al., 2000), (v)ISBA-SURFEX (Boone et al., 2000; Decharme et al., 2011, 2016;
Decharme and Douville, 2006), (vi) Noah-MP version 3.0 (Niuet al., 2011; Schaake et al., 1996), (vii) JULES version 4.6 (Best etal., 2011; Clark et al., 2011), (viii) CLM4.5 (Oleson et al., 2013),(ix) CABLE (Decker, 2015; Kowalczyk et al., 2013, 2006; Wanget al., 2011), (x) SSiB (Li et al., 2010; Sun and Xue, 2001; Xue etal., 1991, 1996; Zhan et al., 2003), (xi) CH-TESSEL (Balsamo etal., 2009; Entekhabi and Eagleson, 1989), and (xii) JSBACH version 3.0 (Hagemann and Stacke, 2015; Roeckner et al., 2003). All
models are listed in Table 2 and the labeling is used throughout
the appendix for consistency.
First, we review how these LSMs address the challenge oftreating infiltration processes in a large grid cell and associatedrunoff generation. Second, we review how infiltration processes
are treated numerically in different LSMs, considering approaches
ranging from empirical to analytical equations or direct solutionof the Richards equation with appropriate numerical schemes,boundary conditions, and vertical discretization. These approaches
have been introduced above.

Grid-Scale Infiltration Processes
Here we focus on parameterization of infiltration and runoff
generation processes and their definition at the grid scale. These
include the maximum infiltration rate, the characterization of
the soil water content distribution at the grid cell level, the characterization of grid cell heterogeneity with respect to infiltration
controls, and the saturated surface fraction. Although most LSMs
use a Richards type formulation to describe the soil water flow at
the grid scale, it is used primarily to redistribute soil water verticallyin the profile rather than for calculation of infiltration rates at eachtime step. This is justified by the computationally demanding solution of the pressure head at the land surface required for obtaining
accurate and stable solutions of the Richards equation. This calculation is particularly demanding during high-intensity rainfall
events, where accurate solutions require very small spatial and temporal discretization (see above) leading to very high computational
demand. Various parametric approaches have been devised in LSMs
to handle the process of infiltration efficiently within the constraintsof data availability (e.g., spatially available Ks values) to avoid such a
computational burden. However, the implementation of a Richards
based solution for infiltration allows (i) direct calculation of excesswater flux (e.g., rainfall) that cannot infiltrate, (ii) a physically based
determination of the time to ponding, and (iii) direct accounting
for the effect of variable soil properties on infiltration at the gridcell (when such information is available). Table 2 gives an overview
of 12 different LSMs and their numerical treatment of infiltration
processes. Four of the 12 LSMs derive infiltration rates directly fromsolving the Richards equation in its mixed form (ISBA-SURFEX,
OLAM-SOIL, CLSM, and ParFlow-CLM) and thus appear below
p. 16 of 53

VZJ | Advancing Critical Zone Sciencep. 17 of 53

ParFlow-CLM

ORCHIDEE
rev4783

Catchment land 6, 7, 8
surface model(CLSM)

ISBA-SURFEX 9, 10, 11, 12 <1 d (down to 5min)

Noah-MPversion 3.0

JULES
version 4.6

3

4

5

6

7

15, 1613, 144, 52, 31-mm top layer, 11 layers
total

10 layers for first 3 m,
additional 5 layers fornext 100 m‡

2–5 cmfour layers: 0.1, 0.3, 0.6,
and 1.0 m<1 d (typically 1 h flexible, set by user;
but set by user)
standard top layer =
10 cm

<1 d

mixed form

flexible but never >1 d;
usually 1–6 hdiffusivity form

offline mode 30 min–1 h; diffusivity formin coupled mode down
to seconds

0.01 m (by default, but user max. 3 h but can be smaller mixed form
can modify this)depending on available
data

hourlydiffusivity form

mixed form

variable

30 mindiffusivity form

minutes

Vertical resolution of the Temporal resolution of
first compartment
precipitation input<1 d (typically set 1 cm for Richards’
to <1 h)equation in off-line
simulations, 2–5 cm
for catchment-scale
simulations<1 d (30 min)variable<1 d (?1 min–6
h)

2

1

OLAM-SOIL

1

Key
Model timereferences† stepping

LSM

No.

no

no

yes¶

yes§

no

yes

yes

Infiltration
rate from
Form of Richards’
solution ofequation used for 1D Richards’
water movement
equation

Table 2. Overview of land surface models (LSMs) and their numerical treatment of the infiltration process.

no

no

no

no

no

yes

yes

Approaches to pointscale infiltration

Approaches to grid-scale
infiltration

runoff occurs if
precipitation exceeds
Ks

runoff occurs if
precipitation exceeds
Ks

exponential pdf on Ks

uniform upscaling of
infiltration

runoff occurs if
precipitation exceeds
Ks (Horton) or soil issaturated (Dunne)
when TOPMODEL
runoff occurs if
approach is used, there
precipitation exceeds
is an additional deep
Imax = bKs (Horton)water store below the soil
or soil is saturated
column, in which saturated (Dunne)
conductivity is assumed toreduce exponentially withdepth

no impact due to freedrainage LB

runoff occurs if
default is free drainage(hydraulic conductivity inprecipitation exceeds
lowest model layer) except
Ks (Horton) or soil issaturated (Dunne)
when ISBA is coupled toa groundwater model inwhich capillary rise canoccur from groundwatercontinued on next page

TOPMODEL approachfor Fsat

exponential pdf on
precipitation, andTOPMODEL approachfor Fsat

Fsat based on reflected
power pdf on soil
moisture capacity(Dümenil and Todini,1992) or TOPMODEL
approachwater table depth influences runoff occurs if water
Richards’ equation
the soil moisture profile
table depth is above
combined withand therefore infiltration
surface or precipitation TOPMODEL
rates
exceeds infiltration
formulations tocapacity
derive parameters forcatchment-scale vertical
moisture transfer

by controlling the soil
moisture profile andtherefore infiltration
capacity

Richards’ equation is solved
globally (saturated
and variably saturated
conditions); different
types of BCs can be used
(Dirichlet, Neumann,
gravity drainage)

uniform upscaling of
LB of soil model is
Richards’ equation isimpermeable but >100 m
coupled with surface
infiltration
deep; water table governedwater (ponding) upperby Richards’ equation
boundary condition;
infiltration occurs
when infiltration
capacity will be
exceeded

Is ponding
allowedbased onpressure
Impact of lower boundary
head?
(LB) on infiltration

VZJ | Advancing Critical Zone Sciencep. 18 of 53

CLM4.5

CABLE

SSiB

CH-TESSELJSBACHversion 3.0

8

9

10

11

12

<1 d (30 min)29, 3027, 2822, 23, 24,
25, 26<1 d (hourly)

1 h or less

<1 ddiffusivity form

diffusivity form

6.5 cmhourly, daily

diffusivity form

no

no

no

no

(form is not provided) no

increasing layer thickness flexible but never longer
diffusivity form
with depth, upper layer
than 6 h offline;= 7 cmtypically 3-hourly or
hourly offline and£atmospheric time stepwhen coupledincreasing layer thickness 30 min or hourlywith depth, upper layer= 1.75 cm

depending on forcing
data or coupling modeltime stepincreasing layer thickness 3–6 hwith depth, upper layer= 1.75 cm

Vertical resolution of the Temporal resolution of
first compartment
precipitation input

18, 19, 20, 21 fixed, depending increasing layer thickness
on forcing datawith depth; first layer
or coupling
2.2 cm; soil thickness
model time step 4.6 m

17Key
Model timereferences† stepping

Infiltration
rate from
Form of Richards’
solution ofequation used for 1D Richards’
water movement
equation

no

no

no

no

yes

runoff occurs if
precipitation exceeds
Ks (Horton) or soil issaturated (Dunne)

Approaches to pointscale infiltration

TOPMODEL approachfor Fsat

Approaches to grid-scale
infiltration

no impact of LB

no impact of LB

runoff occurs if soil
moisture exceeds
capacity

runoff occurs if soil
moisture exceeds
capacity

by interaction of baseflow
runoff occurs if
and the overlying soil layers precipitation exceeds
Ks

Fsat based on reflected
power pdf on soil
moisture capacity(Dümenil and Todini,
1992)

Fsat based on reflected
power pdf on soil
moisture capacity(Dümenil and Todini,
1992)

exponential pdf on
precipitation

infiltration depends on
air space in first three soil gamma pdf on soil
vertically averaged relative
layers is filled from themoisture
saturation over the grid
top with infiltration
cell, which includes thewater at the beginninglowest soil layer.
of time step

no impact of LB

Is ponding
allowedbased onpressure
Impact of lower boundary
head?
(LB) on infiltration† 1, Walko et al. (2000); 2, Kollet and Maxwell (2006); 3, Kollet and Maxwell (2008a); 4, d’Orgeval et al. (2008); 5, Ducharne et al. (2017); 6, Koster et al. (2000); 7, Ducharne et al. (2000); 8, De Lannoy et al. (2014); 9, Boone et al. (2000); 10, Decharmeand Douville (2006); 11, Decharme et al. (2011); 12, Decharme et al. (2016); 13, Niu et al. (2011); 14, Schaake et al. (1996); 15, Best et al. (2011); 16, Clark et al. (2011); 17, Oleson et al. (2013); 18, Kowalczyk et al. (2006); 19, Kowalczyk et al. (2013);
20, Wang et al. (2011); 21, Decker (2015); 22, Xue et al. (1991); 23, Xue et al. (1996); 24, Zhan et al. (2003); 25, Sun and Xue (2001); 26, Li et al. (2010); 27, Balsamo et al. (2009); 28, Entekhabi and Eagleson (1989); 29, Roeckner et al. (2003); 30,
Hagemann and Stacke (2015).
‡ ParFlow-CLM vertical discretization is in general flexible and can be set by the user, the given example is used as the standard setting in TerrSysMP.
§ Parameters for catchment-scale vertical moisture transfer between surface and root zone, and root zone and catchment deficit, are derived from detailed Richards’ equation simulations conducted off-line combined with TOPMODEL formulations.
¶ Calculation of infiltration rate from solution of Richards’ equation is done only for small time steps. For large time steps and intense rain events, a Green–Ampt-based approach can permit infiltrated water to attain several layers during a time step.

LSM

No.

Table 2. continued.

only when dealing with specific issues of grid-scale parameterizationof infiltration and runoff generation. Here it has to be noted thatmost reviewed LSMs (except OLAM-SOIL and Parflow-CLM) aredesigned to be used globally, with a typical grid cell size of 0.5° or
more. These models first define uniform soil parameters at the gridcell scale (Table 3), used to calculate either infiltration or Hortonianand Dunne runoff at the point scale. In the latter case, grid-scale
infiltration is defined as the incoming water that does not run off.

The Maximum Infiltration Ratein Land Surface Models
In LSMs, the maximum infiltration rate is used to partition
the water flux reaching the land surface into infiltrable flux andexcess water that generates surface runoff. We focus on infiltrationof vegetated or bare land surfaces and exclude glaciers, lakes, andurban areas. An interesting example was provided by Entekhabi
and Eagleson (1989), who used a Darcy-based approach to estimate
the maximum infiltration rate or infiltrability for a heterogeneous
grid cell in a LSM. They defined surface runoff, qs, as the sum of
Hortonian infiltration excess (Px − f * for Px > f * and s < 1, i.e., thefirst term on the right-hand side in Eq. [36]) and Dunne saturation
excess (the second term on the right-hand side of Eq. [36]), where
saturation excess is equal to Px when s ³ 1:
1é 1 ¥
qs = ê ò ò ( Px - f * ) f Px ( Px )d Px f s ( s )d s
k êë 0 f *¥ ¥ù+ ò ò ( Px ) f Px ( Px )d Px f s ( s )d s ú
1
0
úû

[36]

Table 3. Overview of approaches used for the spatial distribution of saturated hydraulic conductivity, Ks, in the land surface models (LSMs).

LSM

Pedotransfer function to estimate
for grid-scale Ks

Underlying soil maps to provide coverage of Ks

Use of soil
structural
information

Ks
Variation
constant in of Ks withgrid cell
depth

SoilGrids (Hengl et al., 2014)

yes

yes

yes

OLAM-SOIL

Weynants et al. (2009) or de Boer(2016)

ParFlow-CLM

Schaap and Leij (1998) for first 10
FAO (1988) and Gleeson et al. (2011)
soil layers; Gleeson et al. (2011) for
deeper layers

no

yes

yes

ORCHIDEE rev4783

lookup tables for van Genuchten soil Zobler (1986) at 1° is default soil texture map; other soil
parameters for each USDA classmaps can be used: Reynolds et al. (2000) at 1/12° or
SoilGrids1km (Hengl et al., 2014)

no

no†

yes

Catchment Land Surface Campbell (1974) or Wösten et al.
Model (CLSM)(2001)

Harmonized World Soil Database version 1.21 (FAO,
2009) and the State Soil Geographic (STATSGO2)
database (De Lannoy et al. 2014)

no

yes

yes

ISBA-SURFEX

from Clapp and Hornberger (1978)or Cosby et al. (1984) (optional);
for organic soils (Decharme et al.,
2016)

Harmonized World Soil Database (FAO, 2009) at 1-km
resolution

no

yes

yes

Noah-MP version 3.0

from Clapp and Hornberger (1978)

Harmonized World Soil Database (FAO, 2009) or Milovac noet al. (2014)

yes

no

JULES version 4.6

from Clark and Gedney (2008);
for organic soils, Chadburn et al.
(2015b)default configurations of JULES uses soil ancillaries based no
on sand, silt, clay fractions from the Harmonized World
Soil Database using the function of Cosby et al. (1984)

yes

yes

CLM4.5from Clapp and Hornberger (1978) based on IGBP soil dataset for mineral soils (IGBP, 2000), no
and Cosby et al. (1984); for organic with organic soil described from ISRIC-WISE (Batjes,
soils, Lawrence and Slater (2008)
2006) and the Northern Circumpolar Soil Carbon
Database (Hugelius et al., 2013)‡

yes

yes

CABLE

same as horizontal

CABLE rev 2.0 uses Zobler soil class information (Zobler, no
1986, 1999); Decker (2015) uses sand, clay, silt fractions
from the Harmonized World Soil Database (FAO,
2009), calculating soil properties from the textures as in
CABLE rev 2.0 from Zobler (1986)

yes

no

SSiB

from Clapp and Hornberger (1978)
and Sellers et al. (1986)

Harmonized World Soil Database v.1.1 (FAO, 2009)

no

yes

no

CH-TESSEL

based on tabulated values for 7
textural classes including organicsoils (Balsamo et al., 2009)

FAO (2003) soil texture map

no

yes

yes

JSBACH version 3.0

from Beringer et al. (2001)

based on improved FAO soil type dataset (Hagemann and no
Stacke, 2015)

yes

no

† Based on exponential pdf (see Table 3 and Appendix A, A2.6).
‡ The International Geosphere–Biosphere Program soil dataset (Global Soil Data Task, 2014) of 4931 soil mapping units and their sand and clay content for eachsoil layer were used to create a mineral soil texture dataset (Bonan et al., 2002). Soil organic matter data are merged from two sources. The majority of the globe isfrom ISRIC-WISE (Batjes, 2006). The high latitudes come from the 0.25° version of the Northern Circumpolar Soil Carbon Database (Hugelius et al., 2013). Bothdatasets report C down to 1-m depth.

VZJ | Advancing Critical Zone Sciencep. 19 of 53

VZJ | Advancing Critical Zone Science

no
ORCHIDEE
rev4783

amended TOPMODEL
NA
approach: grid cell soil
moisture variability tied
to the catchment watertable depth distribution;
moisture distribution
shifted by root zone waterin excess (or deficit) of
equilibrium conditions continued on next page
amended TOPMODEL approach: local
max. infiltration rate not imposed per se in model concept;
water table depth related to local
it is a complex outcome of off-line Richards’ equation
topographic index, mean topographic
simulations, TOPMODEL formulations, and the model’s
index, and mean water table depth;
dynamically varying spatial moisture fields (see Appendix
saturated fraction under equilibrium
A2, A2.4)
equals the land fraction with a water table
above ground surfacenonoCatchment Land amended TOPMODEL
Surface Modelapproach: grid cell soil
(CLSM)
moisture variability tied tothe catchment water table
depth distribution; moisture
distribution shifted byroot zone water in excess(or deficit) of equilibrium
conditionsNA
not described for thehorizontal variability
no Fsat, i.e., no saturation-excess runoff, but max. infiltration rate not imposed per se in model concept; itponding can occur because of infiltrationis the outcome Ks distribution and infiltration from top tobottom via a modified Green–Ampt model (see Appendix
excess runoff (see Appendix A3, A3.3)
A2, A2.3)
exponential
pdf

no
no
no
no
ParFlow-CLM

noponding is permitted, and the hydraulic head homogeneous within grid box
at the bottom of the pond is the Dirichlet
upper boundary condition for Richards’
equation in the soil layers; infiltration isthe direct solution of Richards’ equationnonono
OLAM-SOIL

no

NA
NA

NA
NA†

Grid cell soil moisturevariability equation used
Distributionof Ks
Saturated fraction of grid cell
Distribution of
precipitation
Distribution of soil moistureand h is the matric potential [L], Dz is the thickness of the first soil
layer [L], s is the effective relative saturation (q/q s), and Ks is thesaturated hydraulic conductivity [L T−1]. For gravitational flow, v= 0 and hence f * » Ks. It is important to state that Ks is assumed
to be uniform for a grid cell.
Other, more empirically based approaches have been implemented in LSMs to quantify the maximum infiltration rate todetermine the excess water for surface runoff. Land surface models
differ in the way they define this maximum infiltration rate andthe input parameter needed to estimate the infiltration flux. An
overview of the different equations used in LSMs to estimate Imax
is given in Table 4. Appendix A2 provides a more detailed presentation of the concepts and equations. In general, all LSMs thatuse Imax require, in one way or another, information about theactual soil water status of the land surface. In Noah-MP and CLM,
this is embedded in the saturated water fraction, Fsat (see below).
H-TESSEL and CH-TESSEL require knowledge about the actualand maximum water content of the first 50 cm of the soil profile
to calculate Hortonian infiltration. ORCHIDEE, JULES, andISBA-SURFEX require information on the soil hydraulic conductivity to estimate Imax. CLSM does not use the concept of a
maximum infiltration rate but assumes that the amount of waterthat can infiltrate in a certain time at the catchment scale is a function of the model’s dynamically varying spatial moisture fields.
Based on the detailed description provided in Appendix A2
on how the maximum infiltration rate is estimated in LSMs, andthe information provided in Table 2 on the numerical simulationof soil water flow at the grid cell scale, we conclude that currently
used LSMs parameterize infiltration processes by (i) estimating thesaturated area in a grid cell to calculate Dunne saturation excess—all precipitation that falls on the saturated surface fraction of agrid cell becomes immediately available for runoff and does not
infiltrate into the soil profile, (ii) calculating the unsaturated land
surface area available for Hortonian infiltration excess, and (iii)
using different approaches to simulate the maximum amount ofTable 4. Overview of definition of parameters and properties affecting infiltration in the different land surface models (LSMs).

where
dh
1
v=
d s s =1 D z

Description or formulation of grid-scale Imax

[37]

LSM

f * = K s vs + K s (1- v )b parameter

where Px [M L2 T−1] is the flux of water incident at the soil surface,f * is the infiltrability of the first soil layer, s is effective relative
saturation, k is a scaling factor needed to redistribute the GCM
grid-scale precipitation across the scale of precipitation events, andf represents the spatial pdf of the respective variables (rainfall andrelative saturation). The first term on the right-hand side refers to
the amount of point precipitation intensity that exceeds the infiltration rate of the soil, f *, and therefore represents the maximum
infiltration rate at that moment in time. The second term on theright-hand side refers to the rainfall that falls on saturated surfaces
and cannot infiltrate.
Based on the Buckingham–Darcy equation, they derived f * asp. 20 of 53

VZJ | Advancing Critical Zone Sciencep. 21 of 53

Gamma pdf (Entekhabi andEagleson, 1989)

noreflected power pdf

reflected power pdf

CABLE

SSiB

CH-TESSEL
H-TESSELJSBACHversion 3.0† NA, not applicable.

TOPMODEL approach

CLM4.5nono

exponential pdf

no

no

no

modified TOPMODEL
approach based on local
topographic index

JULES
version 4.6

noexponential pdf

reflected power pdf orTOPMODEL approach

ISBA-SURFEX

Distribution of
precipitation

Noah-MP version TOPMODEL approach
3.0

Distribution of soil moisture

LSM

Table 4. continued.

no

no

no

no

no

no

no

no

æ
ö
ç K s D t ÷÷ c SSiBlog çç
÷bSSiBçè Pdrop aSSiB ø÷÷ aSSiB
1- 0.5 f d ( z wt -z bot )+ FfrzPd + Dx (1- e -kdtDt )Dx (1- e -kdtDt )b +1 ö
æ
1/(b +1)
ïìïæççé P + P ù ïüï ÷÷
W ÷ö
M úï
÷÷
÷÷
+ max çç 0,Wsat ïíççç1- êê Týççïïçè Wsat ÷ø( b +1)Wsat úúû ïï ÷÷÷
ê
ë
çè
îï
þï ø

I max = (Wsat -W )NANA

I max = (1- Fsat ) K s C ice

I max = (1- Fsat )b s K s

I max = Pd

see Appendix A, A2.5

Description or formulation of grid-scale Imax

Arno scheme following Dümenil and Todini for each time step, infiltration is limited to the difference(1992)
between water holding capacity of the root zone and watercontent of the root zone at the beginning of the time stepb /(b +1)æW ö÷÷÷Fsat = 1-ççç1çè Wsat ø÷

Fsat =

æ K ö÷
çp ÷
Fsat = 1- erf çç÷ççè l s ÷ø÷

Fsat = as exp-cs l crit

Fsat = (1- Ffrz ) Fmax e

two options available:
1. Fsat based on a TOPMODEL typeapproach
2. Arno scheme following Dümenil andTodini (1992)

Distributionof Ks
Saturated fraction of grid celll sac ac -1 -l sse
G( ac )NA

æ s -smin ÷ö÷b = max ççç ov
çè sov +smax ÷÷ø

homogeneous within gridbox NA

reflected power density
function

homogeneous within gridbox constants based on
different plant
functional type (Xueet al., 1996)

fs ( s)=Entekhabi and Eagleson(1989)

homogeneous within gridbox NA

é T ( z local ) ùú +l
l local = ln êêmúëê T ( zÑ ) ûú

NA

NA

NA

modified TOPMODEL
approach based on local
topographic index

can be set to 0.2 or 0.5
depending on model
applicationb parameterreflected power density
function

Grid cell soil moisturevariability equation used

water that can enter the soil during a certain time step, any excess
being Hortonian infiltration excess.
With respect to Hortonian infiltration excess, the majority of the LSMs use an empirical parameterization of Imax . The
vertical redistribution of the amount of infiltrated water in mostLSMs is based on a diffusion form of the Richards equation. A
primary consideration in using the simplified infiltration parameterization presented above is the large computational burdenof solving the more accurate pressure head based Richards
equation. Such a direct approach, however, has the advantageof obviating the need for defining the maximum infiltration
capacity, and that the amount of Hortonian infiltration excess
and Dunne saturation excess are immediate outcomes of solving
the Richards equation (albeit subgrid information and runoff–
run-on processes must be represented). CLSM overcomes thecomputational burden of the mixed form of the Richards equation by running detailed Richards equation simulations solely
off-line prior to LSM runs to derive parameters of catchmentscale vertical moisture transfer using an amended TOPMODEL
approach. Also, Walko et al. (2000) developed a global LSM(part of OLAM-SOIL) capable of calculating infiltration processes using a pressure head based Richards equation and highly
resolved spatial information of soil properties.

Spatial Heterogeneity in Soil Water Content
Using Probability Densities
Characterization of subgrid soil moisture variability in LSMs
is especially important to estimate the generation of Dunne saturation excess. We focus here on the different approaches used in
LSMs to quantify heterogeneity in soil water content with respect
to quantifying infiltration–runoff processes at the grid scale level.
These approaches are based on the underlying assumption thatspatial variability in soil water content and derived properties suchas soil water storage or soil water deficit can be described by pdfs
without considering their spatial patterns or specific locations.
Different types of pdfs have been proposed and used in LSMs
such as standard reflected power distribution functions, Gamma
functions, and exponential functions.
Several researchers used standard reflection power distribution functions to describe the spatial variability of soil watercontent and soil water content dependent variables that affectinfiltration. These functions are a special case of the more generalb distribution function written asf ( x;a G , bG ) =G(a g +bG )G( a G ) G( b G )cumulative distribution function to characterize the spatial variation of the storage capacity, c, written as
b
æ
c ö÷
÷÷
F ( c ) = 1-ççç1çè c max ø÷

0 £ c £ c max ; b > 0[39]with the corresponding density functionb-1

f (c )=dF ( c )
b æç
c ö÷
=
çç1÷÷
dc
c max çè c max ø÷[40]

where c is the water storage capacity at a certain location definedas the depth of water that can be stored [L], cmax is the maximum
water storage capacity [L], and b controls the spatial variability ofthe storage capacity across the basin, whereby b = 0 indicates a constant value of the storage capacity and b = 1 means a capacity that
follows a uniform distribution between 0 and cmax. Moore (1985)
did not specify how to calculate the storage capacity from basic soil
properties. Based on Eq. [40], we can derive the maximum water
depth Smax [L] across the basin as
S max = ò

c max

0

c
é1- F ( c )ù d c = max
ë
û
b +1

[41]

The concept proposed by Moore (1985) has been usedin many LSMs using a bucket type of soil water model tocharacterize the spatial variability of infiltration and storage
capacity such as in VIC (Liang et al., 1994), HD (Hagemann
and Gates, 2003), or Richards equation based LSMs (e.g., JULES,
H-TESSEL/CH-TESSEL). Clark and Gedney (2008) compared
the PDM approach of Moore (1985) and a modified TOPMODEL
approach in generating surface runoff using the MOSES (nowJULES) LSM by comparing model output against the observed
stream flows in three catchments. TOPMODEL performed best
as it allowed a better response to subsurface flow contributing topeak flows but also capturing slower changes in recession times.
PDM only improved the calculation of the surface runoff without
improving subsurface flow. In Appendix A1 we briefly describe the
PDM scheme of Moore (1985).
Rather than using a ref lected power density function,
Entekhabi and Eagleson (1989) proposed a two-parameter g pdf
to characterize the spatial heterogeneity of soil water content at
the level of the grid cell:f s ( s; a c , l s ) =l s a c a c -1s
exp-ls ,G( a c )ac ,ls , s > 0[42]where s is the surface layer point soil water saturation defined byx

a G -1bG -1(1- x )

[38]

0 £ x £1; a G , bG > 0
where G(…) is the Gamma function and a G and b G are the parameters of the Gamma function. Moore (1985) was among the first to
propose the probability-distributed model (PDM) as a concept tocharacterize the spatial variability of soil water content related variables at the catchment scale. He used a standard reflected power

VZJ | Advancing Critical Zone Science

s=

q
qs

[43]

where q is the actual and q s is the saturated volumetric soil watercontents [L3 L−3], respectively. The two-parameter Gamma distribution, G, in Eq. [42] is related to s via
l=

ac
s

[44]

p. 22 of 53

where áf (s)ñ is the grid mean relative saturation of the surface soil
layer [L3 L−3]. This concept has been implemented, for example,
in the CABLE model (Decker, 2015).
For the CLSM, the basic land element is the (irregularly
shaped) hydrological catchment rather than the GCM grid cell.
A number of these catchments lie within a given cell, and grid-cellfluxes are computed through areal weighting of the component
catchment fluxes. Within an individual catchment element,
CLSM simulates a dynamic water table depth of which the distribution is related to catchment topography characteristics using
the TOPMODEL formulation (Eq. [A25]) (Beven and Kirkby,1979). Each catchment is characterized by its topographic index
distribution, which in effect is used to diagnose the spatial variability of soil moisture within the catchment from the catchment
element’s three bulk water prognostic variables: catchment deficit representing the catchment air-filled pore volume assuminghydrostatic equilibrium conditions for a given water table depth
distribution, and surface layer excess and root zone excess representing the water in excess (or in deficit) of the correspondinghydrostatic equilibrium conditions. The information on soil
moisture variability is used to define three distinct hydrologicalregimes: (i) the wilting area, (ii) the subsaturated-but-transpiringarea, and (iii) the saturated fraction. The definition of these three
regimes is key to CSLM because different evaporation and runoff
physics are applied to each.

Representing Spatial Heterogeneityof Surface Properties that Control Infiltration
in a Single Parameter
Several LSMs use a single parameter in combination with pdfs
to characterize the spatial heterogeneity of infiltration processes.
Many of these pdfs, presented above, can be derived from Eq. [38],which was used above. to describe the spatial variability of soil watercontent and related variables. The estimation of this exponent(denoted as B or b) in such pdfs remains an open question because it
cannot be immediately derived from available soil properties. Please
note that this exponent has sometimes received a different notation
in some LSMs. Liang and Xie (2001) used the b notation to describe
the variability of soil and grid cell properties affecting Dunne excess
saturation in the VIC model and the B parameter to describe thevariability of soil and grid cell properties affecting Hortonian infiltration excess in the VIC model. Three main approaches have beenreported in the literature to determine this exponent, from here onreferred to as the b parameter for simplicity.
In the first approach, the b parameter is obtained by model
calibration to available hydrological time series. Huang et al.
(2003) present a brief discussion of the literature dealing withthe estimation of the b parameter in the VIC model. They
showed that the calibration approach suffers from the problemof equifinality, as so many other LSM calibration exercises do.
Furthermore, they did not manage to establish meaningful relationships between model parameters and physical characteristicsof the catchment or regions.

VZJ | Advancing Critical Zone ScienceIn the second approach, the b parameter is derived from available soil and/or topographic information. For example, Dümenil
and Todini (1992) suggested calculating the b parameter from thesubgrid standard deviation of topography byé s -sminùb = max ê h;0.01ú
ê s +s
ú
max
ë h
û

[45]

where sh refers to the standard deviation of the topography withina model grid cell. Balsamo et al. (2009) stated that s h varies
between 0.01and 0.5 and smin and smax can be set to values proposed by van den Hurk and Viterbo (2003). On the other hand,
Habets et al. (1999), in this case for the ISBA model, assumed theb parameter to be constant.
The third option fits a pdf to the observed cumulative distributions of soil properties as a function of the occupied space in agrid cell or basin. Sivapalan and Woods (1995) fitted the so-called
Xinanjiang distribution to data of soil profile depth to estimate
the b parameter in
b
æ
Z ö÷
÷÷
Fs ( z * ) = 1-ççç1è Zm ÷ø[46]

where Fs is the cumulative distribution function of the scaled
infiltration capacity (scaled by its maximum value), Z is the soil
profile depth to the bedrock [L], and Zm is its maximum valueof the bedrock depth within the grid cell [L]. Here, it is assumed
that the available porosity is constant in space and time. To provetheir concept, Sivapalan and Woods (1995) used a regionalized
cumulative soil depth distribution for six landforms observed in
the Serpentine catchment and obtained the best fit with a valueof b = 4.03 and Zm = 10 m.
Huang et al. (2003) used a method based on a self-organizing
neural network map combined with a K-means clustering method
to develop transfer functions that are able to transfer b values ofdata-rich to data-poor areas. The data-rich area was used to fit therelation to values of soil water capacity data derived from soil map
information to the surface areas of catchments by
1/b
w = wmax éê1-(1- A ) ùú[47]
ëû
where w and wmax are the point and maximum point soil watercapacity, A is the area for which the soil water capacity is lessthan or equal to w, and b is defined as a soil water capacity
shape parameter. In their study, Huang et al. (2003) used the
STATSGO (Digital General Soil Map of the United States)
database to represent the data-rich areas. Here, it has to be
noted that appropriate soil information at resolutions <1 km
is presently available, whereby the global SoilGrids1km 250-m
database of Hengl et al. (2014) is only one example. These newdevelopments will allow estimating the value of the b parameterat high resolution globally, without the need for applying transfer approaches as proposed by Huang et al. (2003). Additionally,such highly resolved global map information combined with
PTFs would allow the b parameter to be directly related tothe soil properties to be determined, as well as the subgrid

p. 23 of 53

variability depending on the size of the grid cells used in LSMs.
This is discussed below.

Estimating the Areal Saturation Fraction
Estimation of the areal saturation fraction, Fsat, as used forexample in the equations listed in Appendix A2 (Eq. [A24], [A35],[A36], [A40], [A48], [A50], [A52], [A56]), is key in determining
the contribution of Dunne saturation excess runoff and Hortonian
infiltration excess runoff to the overall runoff generation of agrid cell in many LSMs. Here we discuss briefly the different
approaches used in various LSMs. A correct representation of Fsat
typically requires a good knowledge about subsurface properties
like the depth to bedrock, the location of the groundwater table,soil porosity, and the actual state of the soil water content in the
profile. One approach that is frequently used in LSMs is basedon the assumption that the saturated fraction of a grid cell canbe determined from topographical characteristics and the soil
moisture status of a grid cell. This approach is closely related tothe concept introduced in TOPMODEL and adopted in manystudies of infiltration runoff generation. It has been used in LSMs
such as CLM, JULES, Noah-MP, CSLM, and ISBA-SURFEX
among others. Some LSMs require only information on the soil
moisture distribution such as CABLE. ORCHIDEE does not use
the concept of Fsat but introduces the notion of a ponded fraction,
which serves as the opposite of Fsat. Instead of increasing runoff, it
enhances infiltration since it allows the latter to develop over several time steps. JSBACH uses the Arno scheme to calculate surface
runoff and infiltration and uses soil water capacity to determine
Fsat. The equations used to calculate Fsat are listed in Table 4, anda more extensive description of the use of Fsat in various LSMs,
including the concepts used in TOPMODEL, are presented inAppendix A2.

Characterizing the Surface Saturated
Hydraulic Conductivity
The Ks is a key soil hydraulic property that controls soil water
fluxes and thus the infiltration of water in soils. However, Ks is ascale-dependent parameter that strongly depends on the support
size of the measurement (Ghanbarian et al., 2015), and it exhibits anisotropic behavior (Pachepsky and Hill, 2017). Despite its
importance, point-scale Ks is a parameter that is still very difficult
to estimate from PTFs that are typically based on simple basic soil
information such a soil texture and bulk density. However, thisinformation only partly explains the variability observed in pointscale Ks values. In the widely used ROSETTA software (Schaap etal., 2001), which contains soil hydraulic properties and basic soil
information for various soils and is used to estimate soil hydraulicproperties, the logarithmically transferred point-scale Ks value isestimated with a coefficient of determination (R2) of 53.5% whenusing bulk density, sand, silt, and clay percentages. Including soil
moisture content at field capacity and the wilting point, the R 2
increased to 64.7%. The HYPRES database of soil hydraulic properties for European soils developed by Wösten et al. (1999), which

VZJ | Advancing Critical Zone Science

contains information on 1136 horizons for hydraulic conductivity,
provided a PTF to estimate point-scale Ks with an R2 of 19%. Thisfunction uses silt and clay content, bulk density, organic matter
percentage, and whether the soil sample used to derive the properties originated from the topsoil or subsoil. A similar R 2 value wasobtained by Vereecken et al. (1990) using clay, sand, and carbon
content, as well as bulk density, for 127 undisturbed point-scale
samples. To further increase the prediction of point-scale Ks, structural properties in addition to currently used soil properties need
to be considered (Vereecken et al., 2010).
Currently, PTFs for point-scale Ks are based on relatively
small sample volumes. In addition, point-scale Ks shows a highspatial variability as it is not only controlled by textural properties
but also by soil structural properties and by the management ofsoils (Strudley et al., 2008; Van Looy et al., 2017). Some of these
aspects are elaborated in more detail below. Here we review thedifferent approaches that currently are being used in LSMs to estimate Ks and to represent its spatial distribution at the grid scale.
Spatial Distribution of the Saturated Hydraulic Conductivity.
Table 3 provides an overview of the approaches used to quantify
grid-scale Ks in 12 different LSMs. It must be stressed that theseestimations of grid-scale K s are based on PTFs (using the texturalclass as input for prediction), although the support scale for thedevelopment of these PTFs is in all cases the point scale.
In all models studied in the analysis except for ORCHIDEE
and ISBA-SURFEX, K s is considered horizontally uniform
across the grid cell. ORCHIDEE considers an exponential distribution of the infiltrability across the grid cell and with depth.
Infiltrability is defined as the arithmetic mean of K s in the
deepest fully saturated layer and the hydraulic conductivity inthe topmost unsaturated layer. ISBA-SURFEX also considers
an exponential distribution of Ks within a grid cell but only for
Hortonian runoff (Decharme and Douville, 2006). As alreadymentioned above, heterogeneous soil surfaces will generate
runoff earlier than homogeneous soil surfaces, which is partly
caused by the heterogeneity in Ks.
With respect to the vertical heterogeneity, JSBACH,
Noah-MP, and CABLE assume a constant Ks value with depth.
CLSM, JULES, ISBA, ORCHIDEE, and OLAM-SOIL allow for
both constant and depth-variable Ks values. Vertical heterogeneityof Ks in CLM is derived directly from vertical heterogeneity in soil
texture. In ORCHIDEE, Ks decreases exponentially with depth
using the Ks value defined at the 30-cm depth as a reference value.
The extent of the exponential decrease depends on an extinction
factor. Furthermore, the Ks profile may be modified dependingon root density. For example, CLSM differentiates between surface and root zone layer hydraulic properties (including Ks) in the
Richards equation simulations that are used in the derivation oftime-scale parameters of catchment-scale vertical moisture transfer.
For subsurface runoff, CLSM uses a TOPMODEL based approachwith a depth-dependent Ks. ISBA can consider up to 14 soil layerswith a varying Ks value.
p. 24 of 53

Estimators for the Saturated Hydraulic Conductivity. In all
LSMs, the estimation of grid-scale Ks is typically based on soil
textural information but using classical PTFs developed for the
point scale (details on these PTFs are provided in Appendix A4).
These estimations are often completed with specific estimation
functions for organic-rich soils (e.g., CLM 4.5, ISBA, JULES, seedetails below and Table 3).
Rahmati et al. (2018) compared the Ks values for the point
scale provided by Clapp and Hornberger (1978) with Ks estimatesderived from field-scale infiltration experiments using the global
soil infiltration database SWIG. The Ks derived from these experiments did not show a dependence on soil textural properties. Onlyfor sand and loamy sand were the values of Clapp and Hornberger(1978) and the SWIG data on the same order of magnitude. For
all other textural classes, SWIG values were larger by a factor of
two to a factor of >1000. Weynants et al. (2009) pointed out that
improved estimates of point-scale Ks may be obtained by including
structure-related information such as pedologic data and information about land use and crop management, which strongly affectsurface soil properties as a result of tillage, root growth, soil trampling by cattle, etc. (see below).
A notable advance is offered by OLAM-SOIL (Walko andAvissar, 2008; Walko et al., 2000), which uses the PTFs developed by Weynants et al. (2009) or de Boer (2016) to estimate Ks.
The PTFs of Weynants et al. (2009) are based on a dataset of 182
point-scale samples of Belgian soils (Vereecken et al., 1989, 1990).
They estimated the texture-dependent Ks value by extrapolating
the unsaturated soil hydraulic conductivity data to saturation andthen calculated the ratio between the textural-dependent Ks andthe measured Ks value, which was also available. This ratio is ameasure of the effect of related structural properties on the measured Ks and was shown to depend on the sand fraction in the soil.
Use of Soil Maps to Estimate the Saturated Hydraulic
Conductivity. The previous section shows that the basic soil information that is used to derive the soil properties needed to estimate
point-scale Ks differs among the models. Importantly, each LSM,
in particular when used for weather forecasting or GCM studies,
can use different soil maps to inform the model on soil hydraulic parameters (so-called ancillary data). For example, CABLE
uses the Zobler soil class information (Zobler, 1986). The version
used by Decker (2015) uses the soil texture from the HarmonizedWorld Soil Database (FAO, 2009), and so does JULES. OLAMSOIL uses the SoilGrids data published by Hengl et al. (2014,2017). The SoilGrids databases provide information on basic soil
properties at seven depths in the soil profile up to 2 m. They also
provide estimates of the depth to bedrock and the distributionof soil classes. The gridded predictions were based on ?150,000
soil profiles and remotely sensed data. The input data used to estimate the soil hydraulic properties by CLSM or ISBA-SURFEX are
obtained from the Harmonized World Soil Database version 1.21
(HWSD1.21) and the State Soil Geographic (STATSGO2) project(De Lannoy et al., 2014). PARFLOW-CLM uses the UNESCOVZJ | Advancing Critical Zone Sciencesoil map of the world (FAO, 1988) to provide a global coverageof Ks values.
The default soil texture map used in ORCHIDEE is the map
presented by Zobler (1986) at 1°, and it is used in the currently
ongoing CMIP6 simulations. CMIP6 refers to the Coupled Model
Intercomparison Project Phase 6 that aims to address (i) how the
Earth system is responding to forcing, (ii) the origins and consequences of systematic model biases, and (iii) the assessment of
future climate change in the context of internal climate variability, predictability, and scenario uncertainty (Eyring et al., 2016).
ORCHIDEE can also use two other soil texture maps. The firstone is the map developed by Reynolds et al. (2000) at 1/12° andthe second is the surface SoilGrids maps at 1 km (Hengl et al.,2014). All these maps are used to define the dominant texture at
the model resolution using the USDA texture classification.
CLM 4.5 uses the soil dataset produced by the International
Geosphere–Biosphere Program (Global Soil Data Task, 2014) thatcomprises 4931 soil mapping units and contains the textural composition (sand and clay content). These data were used by Bonanet al. (2002) to create a global soil map of mineral soil texturaldata. The global data on soil carbon content are obtained from
ISRIC-WISE (Batjes, 2006), whereas the soil carbon content for
the higher latitudes are obtained from the Northern Circumpolar
Soil Carbon Database (Hugelius et al., 2013).
Saturated Hydraulic Conductivity and Infiltration in Organic
Soils. Four models listed in Table 3 adjust the calculation of Ks
for soils rich in organic matter because these soils typically show
higher Ks than mineral soils. Organic soils mostly refer to soils in
the northern latitudes and tropical areas and include permafrost
and peat soils. The first few centimeters of these soils typicallyhave a soil organic C content of up to 100% (Lawrence and Slater,2008). CLM 4.5 uses the approach proposed by Lawrence andSlater (2008) to calculate Ks:
K s, i = (1- f sc, i ) K s,min, i + f sc, i K s,sc, i[48]

where fsc,i is the soil C density of a soil normalized by the soil C
density of peat in the ith layer, Ks,sc,i is the saturated hydraulic conductivity [L T−1] of the organic C fraction in the ith layer, Ks,min,iis the saturated hydraulic conductivity of the mineral fraction [LT−1], whereby Ks,min is obtained from the point-scale PTFs ofClapp and Hornberger (1978) and Cosby et al. (1984) and given by
K s,min, i = 0.0070556´10- 0.884 +0.0153( %sand )i[49]

Values for Ks,sc were given by Letts et al. (2000).
In JULES and ISBA-SURFEX, a geometric average of thepoint-scale Ks is used to account for the effect of soil C following the approach of Chadburn et al. (2015b) and Decharme et al.
(2016), respectively, with
(1- f sc )K s = K s,min
K s f,scsc[50]

For CLM 4.5, soil organic matter data are obtained fromthe Harmonized World Soil Database, except for the northernp. 25 of 53

high latitudes, which come from the Northern Circumpolar
Soil Carbon Database (Hugelius et al., 2013). In ISBA, the soilorganic C content is determined from two soil horizons (0–30
and 30–100 cm) of the HWSD1.21.
In JSBACH, three additional classes are added to the 11 soil
textural classes, namely peat, moss, and lichen, each with an average value of Ks according to Beringer et al. (2001).
In CLSM, soils are stratified into four levels of organic C
content (De Lannoy et al., 2014), with the highest one (>8.72%
C, i.e., >15% organic matter) representing peat. Peat parameters
were taken from the point-scale information provided by Wösten
et al. (2001) that represents highly decomposed peat. In a recently
developed peatland module for CLSM, parameters and model
structure were further revised, accounting for the high macropore
fraction of undecomposed peat by allowing direct infiltration to
the water table when the topsoil layer becomes saturated, i.e., effectively turning off the Hortonian runoff mechanism over peatlands(Bechtold et al., 2019).

In the original JULES version listed in Table 2, soil water
fluxes are calculated by numerically solving the diffusivity formof the Richards equation for the soil water content increment
Dq, where K(q) and h(q) are given by either Brooks–Corey (Eq.
[4–6]) or van Genuchten (Eq. [7–8]) parametric models. The
choice of parametric model is set by the user. Note that currently
the JULES van Genuchten scheme uses parameters (hc and nvg)
that are directly derived from the Brooks–Corey soil ancillaryparameters as used for the Cosby et al. (1984) K(q) and h(q) equations, rather than using the van Genuchten-specific PTFs given
by Wösten et al. (1999), for example. The numerical scheme uses
an implicit “forward time step weighting” for numerical stability,in which the water fluxes are first calculated as a first-order finite
difference scheme using Eq. [3], but then the moisture increments
are recalculated with K and h given by the water contents after
the time step, effectively solving (see Best et al. [2011] for details)

Numerical Treatment of Infiltrationin Land Surface Models

The vertical discretization in JULES is flexible and set by
the user. For the standard operational configurations and current
Earth system model, only four soil layers are used, with thicknessesof 0.1, 0.25, 0.65, and 2 m. For greater numerical accuracy, a finer
discretization has been applied (e.g., Chadburn et al., 2015b).
Model time stepping in JULES is typically <1 h but is also set by
the user. Recently, Haverd et al. (2016) and Cuntz and Haverd(2018) implemented a new soil (Haverd and Cuntz, 2010) andsnow model with physically accurate freeze–thaw processes within
CABLE, which solved the Richards equation in the mixed form.
In JSBACH (Hagemann and Stacke, 2015; Roeckner et al.,2003), the diffusive form of the Richards equation is used withhourly model time stepping. The soil profile is discretized in thefirst compartment with 6.5-cm layer thickness, and vertical infiltration is calculated based on the ARNO scheme according to
Dümenil and Todini (1992), whereby water will infiltrate while
precipitation is below the difference between the local storage
capacity and the initial water content within the root zone. In
JSBACH, no ponding is allowed and the lower boundary has no
impact on the infiltration rate.
In SSiB (Sun and Xue, 2001; Xue et al., 1991; Zhan et al.,2003), the water flow is solved by the diffusivity form of the
Richards equation, whereby the soil is discretized nonuniformly,
with smaller layer thicknesses close to the surface (upper layer =1.75 cm). Time stepping is <1 h, and runoff will occur if precipitation exceeds Ks. No ponding of water is allowed at the soil surface.
Noah-MP (Niu et al., 2011; Schaake et al., 1996), JULES(Best et al., 2011), and CABLE (Decker, 2015; Kowalczyk etal., 2013, 2006) also solve the Richards equation in its diffusivity form, and none of these models calculate the infiltration rate
directly from the Richards equation. Vertical discretization varies
between 2.2 cm (CABLE) and 10 cm (Noah-MP and JULES) for
the upper layer, and time stepping is < 1 d for all models (see Table
2). No ponding is allowed for any of these three models, and for

The precipitation simulated in LSMs either infiltrates into thesoil or becomes runoff when the soil is saturated (Dunne overland
flow) or if the precipitation rate exceeds the infiltration capacity ofthe soil (Hortonian runoff). As discussed above, various modeling approaches have been developed to describe infiltration andthe transition to surface runoff. In the following, we provide an
overview of the numerical treatments used by the LSMs listed inTable 2.
The version of ORCHIDEE described here (Ducharne et al.,2017) solves water infiltration and redistribution based on the diffusivity form of the Richards equation, using an 11-layer vertical
discretization (de Rosnay et al., 2002). Since d’Orgeval et al. (2008),
point (profile) scale infiltration is no longer calculated based on the
Richards equation; instead, it invokes a piston-like wetting front
inspired by the Green and Ampt (1911) formulation to simplify
the boundary conditions under time-varying rainfall rates and soil
moisture profiles. In ORCHIDEE, the speed of the wetting front
propagation is simplified in two ways, which both tend to reduce
infiltration compared with the classical Green and Ampt formulation: (i) the suction at the front is neglected, and (ii) the hydraulicconductivity at the wetting front is not the soil Ks but an arithmetic
mean of Ks in the lowest fully saturated layer and K(q) in the topmost unsaturated layer (the resulting value is termed Ki). An iterative
procedure is used to account for saturation of one soil layer after
the other, and infiltration-excess (Hortonian) runoff results when
the rainfall rate cannot infiltrate during a single ORCHIDEE time
step of 30 min. In this framework, infiltrability depends on the profiles of Ks and K(q), which itself depends on Ks and q based on the
Mualem–van Genuchten model (Eq. [7] and [8]). The value of Ks
at the 30-cm depth is taken from Carsel and Parrish (1988) for the
12 USDA soil textural classes, but Ks exponentially decreases withdepth following Beven and Kirkby (1979).

VZJ | Advancing Critical Zone Science¶ h( q +dq ) ù
¶q
¶ éê
=
+ 1ú
K ( q +dq )ú¶ t ¶ z êë
¶z
û

[51]

p. 26 of 53

Noah-MP and JULES the lower boundary does not affect infiltration, whereas for CABLE the infiltration will be modified bychanges in the soil moisture profile due to groundwater influence.
Another set of LSMs solve the Richards equation in its
mixed form such as Parflow-CLM (Kollet and Maxwell, 2006,
2008b), OLAM-SOIL (Walko et al., 2000), ISBA (Boone et al.,
2000; Decharme and Douville, 2006), CLSM (Ducharne et al.,
2000; Koster et al., 2000), and CABLE enhanced by Haverd etal. (2016) and Cuntz and Haverd (2018). Because these models
solve the Richards equation in the mixed form, the infiltration
rate can be directly computed from solving the Richards equation.
Nevertheless, ponding is only allowed in OLAM-SOIL, ParflowCLM, and the enhanced CABLE.
In CLSM, infiltration rates are not solely an outcome of
the Richards equation simulations. They are one component ofan overall catchment-scale model concept, with other processes
affecting infiltration rates. CLSM uses a non-traditional framework that strongly emphasizes the subgrid horizontal variabilityof the land surface hydrological processes. The results of the onedimensional Richards equation simulations are combined with
TOPMODEL formulations that control the varying water table
depth and moisture fields at the catchment scale. From this combination, time-scale parameters of catchment-scale vertical moisturetransfer are derived. It is to be emphasized that in this approach
the spatial water table depth distribution is an important factor
influencing catchment-scale infiltration rates.

66Sensitivity of Infiltration–Runoff

Process to Model Parameters

Most sensitivity studies that have been performed with LSM
models with respect to infiltration–runoff processes have focusedmore on the analysis of runoff and river discharge (e.g., Huang etal., 2017; Materia et al., 2010; Zhang et al., 2016) and less on thesensitivity of the infiltration process to model parameters. Here we
review the results from sensitivity studies on LSMs that provide
information on key parameters controlling infiltration processes
and thus ultimately the whole water and energy balance.
One of the first studies to analyze the sensitivity of LSMs toinfiltration processes was conducted by Dirmeyer and Zeng (1999).
They analyzed the sensitivity of infiltration to the treatment of
convective precipitation and the choices made with respect tothe vertical resolution of the soil profile and soil properties. They
found that the choice of the thickness of the surface soil layer
impacts the simulation of infiltration, with thinner surface layers
causing infiltration excess to be more likely because the thinnersurface layer has a much smaller capacity. Basic information about
the impact of the vertical discretization of the hydrological components is also provided above. Unfortunately, in most LSMs thediscretization is predefined and often fairly coarse. In addition,
Dirmeyer and Zeng (1999) found that a “realistic” distributionof convective rainfall in space and time at the grid cell scale isneeded to adequately represent the infiltration and thus surface

VZJ | Advancing Critical Zone Sciencerunoff. In addition, evaporation of intercepted canopy water will
be overestimated if “unrealistic” distributions of convective rainfall
are assumed. They also analyzed in detail the impact of having adepth-dependent soil porosity (they used three layers for the SSiB
model), with a higher porosity for shallow soil layers and lower
porosity (more compaction) for deeper ones. By doing so, theymodeled larger infiltration amounts and reduced gravitational
drainage. Finally, thinner soil layers (2 instead of 5 cm) were foundto generate more infiltration excess, i.e., higher surface runoff,
during high insensitivity rainfall events and soil melt.
Soet et al. (2000) analyzed different conceptualizations of
the land surface scheme and parameter values for three sites with
contrasting soils and climate using the ECMWF TESSEL LSM
developed by Viterbo and Beljaars (1995). A sensitivity analysis,
set up to explore the impact of using standard parameter values
instead of site-specific ones, found that implementing site-specific soil hydraulic properties had a significant effect on runoff
and infiltration at all three sites. On the other hand, the use ofstandard soil parameters led to a systematic underestimation of
evapotranspiration and biases in surface runoff that differed insign for the three different locations.
The sensitivity of the infiltration shape parameter b in the VIC
model (see above) as well as the exponent in the Brooks–Corey equation (Eq. [4–6]) were found to be key for correct representation of
the hydrological system and the partitioning of rainfall between infiltration and runoff under dry soil conditions (Demaria et al., 2007).
Ducharne et al. (1998) found a similar sensitivity to the b parameterfor the bucket model version of ORCHIDEE. On the other hand,
the impact of these parameters on surface runoff generation andstream flow simulations in wet regions was not significant. Shi etal. (2014) analyzed the sensitivity of the catchment outlet discharge
rate to soil properties in the Pennsylvania State University model
Flux PIHM and found an important impact of the van Genuchten
parameters hc and nvg a and n (see Eq. [7] and [8]) on both discharge
rate and soil water content. Runoff simulations of 10 state-of-the-art
hydrological and land surface models including H-TESSEL, JULES,
and ORCHIDEE were compared by Beck et al. (2017), and theyargued for the need to better calibrate, parameterize, and regionalize the parameters of these macroscale models. Most models were
found to generate snowmelt runoff that occurred too early, due toeither the underestimation of precipitation or incorrect descriptionsof input snowfall, snow physics, and meltwater infiltration into thesoil (Bierkens, 2015). Getirana et al. (2014) calibrated river routing
parameters and stated that one of the most important aspects togetting the runoff timing (by 1 or 2 mo!) and runoff volumes right
was the specification of the soil water threshold when runoff occurs(using the Habets and Saulnier [2001] option for runoff in ISBA).
Getirana et al. (2017) extended the study to a large group of LSMs
including those presented here (e.g., CLSM, ISBA, H-TESSEL,JULES, ORCHIDEE) and found that this was also the case for mostof the other LSMs in the ALMIP2 (AMMA Land Surface Model
Intercomparison Project 2; https://www.umr-cnrm.fr/amma-moana/amma_surf/almip2/index.html) ensemble analysis.
p. 27 of 53

In another study, Gudmundsson et al. (2012) compared ninelarge-scale LSMs including H-TESSEL/CH-TESSEL, JULES,
and ORCHIDEE to predict observed runoff percentiles of 426
small catchments throughout Europe and found that the differences in performance among the models became more pronounced
for low runoff percentiles. They concluded that this might be
explained by the uncertainty associated with the representationof hydrological processes, such as the depletion of soil moisturestorage by root water uptake. It is likely that differences in thetreatment of infiltration and calculation of hydraulic properties
will also have played a role. The performance of three LSMs wasanalyzed by Sahoo et al. (2008), including HySSiB, Noah, and
CLM, and they found substantial differences in the prediction of
surface and subsurface runoff for the Little River experimental
watershed, Georgia, United States, which was caused by differencesin the partitioning of the precipitation into infiltration, surface
runoff, and evaporation. An extensive analysis was presented by
Zhou et al. (2012), who compared a set of 14 LSMs (including
VISA, CABLE, ISBA, CLMTOP, and Noah) and six Budyko-typemodels against the observed mean annual runoff from 150 largebasins. They showed that the LSM biases in the prediction of thesimulated mean annual runoff were caused by errors in forcing
data and model parameterizations but also by structural model
errors. The largest biases between the LSM estimates and observed
runoff were found in regions with low mean annual runoff, whichcorresponds with the findings of Gudmundsson et al. (2012).
Hogue et al. (2006) evaluated the model performance andparameter sensitivity for varying levels of LSM complexity acrossfour different biomes using five LSMs including Noah-MP. They
found a large variability among porosity, Ks, and the b parameterused in these models. Based on the impact of these parameterson the simulation results, they advocate either a rigorous calibration or the development and integration of improved vadose zone
water flow models. Especially, calibration of these parameters at
different experimental sites led to differences with respect to thestandard values. A study in the same direction was performed by
Cuntz et al. (2016), who analyzed the role of hard-coded model
parameters (i.e., providing the user with no option to change
values) on the hydrological fluxes in Noah-MP. They found thatthe total runoff was sensitive to both plant and soil parameters(e.g., soil porosity) and that therefore these parameters should beconsidered for calibration. They also stated that surface runoff
is affected by subsurface runoff, which is dependent on available
soil water in the soil profile.
Yang and Niu (2003) compared three different schemes of
topography-based runoff production for the LSM VISA (which
is based on the LSM of Bonan, 1998) and analyzed their sensitivities to key parameters using two catchments. They found that the
decay factor, f, which controls the timing and partitioning of subsurface runoff by rescaling Ks with depth, is a highly importantparameter controlling water table depth and the saturated fraction of the grid cell. Shellito et al. (2016) compared calibrated soil
hydraulic parameters in Noah using in situ soil moisture network

VZJ | Advancing Critical Zone Sciencedata and surface soil moisture from SMOS satellite observations
obtained from seven sites in the United States. Most calibrated
simulations led to higher surface runoff than simulations basedon hydraulic parameters estimated from textural information
using a PTF. The calibrated soil hydraulic parameters included
pore size distribution index, saturated soil water content, Ks ,and the saturated matric potential (or air-entry value, hc) in the
Brooks–Corey equations (Eq. [4–6]). Finally, Yang et al. (2005)
concluded that the characterization of the vertical soil hydraulic heterogeneity is highly important to correctly describe soil
water and soil temperature at the land surface and thus indirectly
infiltration and surface runoff. Based on numerical simulations
and experimental data, they concluded that it was not possible
to replace vertical soil heterogeneity by a homogeneous soil witheffective parameters.
In conclusion, there is relatively little information provided
in the literature on how well the infiltration process and the generation of Dunne or Hortonian overland flow are modeled using
different LSMs and which model parameters mostly impact the
infiltration process. In addition, this literature review indicatesthat it is difficult to identify sources of errors in handling infiltration estimation due to the complexity and the different ways inwhich the infiltration process is being described. One approachto address comparisons of different approaches was proposed by
Clark et al. (2015), who advocated the development of models that
include different parameterizations of the infiltration process sothat parameters and parameterizations can be evaluated in a controlled manner.

66Improving the Infiltration Processin Land Surface Models

The balance between parameterization of a complex heterogeneous soil structure and exogenic processes that affect infiltration
with the operational performance to compute infiltration–runoff
processes that are embedded within LSMs requires an understanding of the trade-off of adding more complex physics to describe
the infiltration–runoff process and the reality of the technicalaspects of computing land surface processes and the determination of related parameters. Also, infiltration and runoff are justtwo of many processes impacting the land–atmosphere interaction.
We therefore aim here to provide an overview to contextualize thecomplexity of the derivation of soil hydraulic parameters rather
than to point out the shortcomings of LSMs in terms of modelinginfiltration and runoff.
In general, there are many soil characteristics broadly related
to soil composition and structure (including macro- and biopores)
and also to exogenic processes, including water repellency, wetting
and drying, swelling and shrinkage, air entrapment, freeze–thaw,
thermal gradients, impermeable layers, and anthropogenic perturbations (e.g., tillage, harvesting) that impact infiltration andrunoff at the point scale (see Young and Crawford [2004] andHannes et al. [2016], for example).
p. 28 of 53

Most of these are presently not considered in most hydrological and land surface models, or not in enough detail, and there islittle ongoing work in this sense such as the implementation ofsoil structure in the OLAM-SOIL model (see Table 3) by the useof a dual-porosity model. Unfortunately, the main challenge in
implementing soil structure into the LSMs lies in the lack of PTFs
considering, e.g., soil structure explicitly and also the temporalchange in soil structure.
In the following, we provide an overview of these processes
and features and discuss the impact on infiltration and runoff
generation briefly.

Soil Structure
The physical soil structure is formed by the combinationof the size, shape, and arrangement of voids and solids, which
ultimately affect water infiltration and runoff, mainly throughthe soil hydraulic properties (water retention and hydraulic conductivity curves) that are part of most LSMs and hydrologicalmodels. In general, the amount of water that infiltrates into asoil is dependent on the available void space (represented by themodel soil layer porosity), which is the cross-sectional area offlow. Greater soil aggregation and pore connectivity increase
bypass or preferential flow, therefore increasing the hydraulicconductivity and movement of water to deeper soil layers (e.g.,
Franzluebbers, 2002; Nissen and Wander, 2003). However, theprocess has not been implemented in most LSMs (Le Vine etal., 2016), apart from the efforts described by Rahman and
Rosolem (2017). The formation of aggregates and the stabilityof the intra-aggregate void spaces is dependent on the rearrangement, flocculation, and cementation of soil and is mediated bysoil organic C, soil biota, ionic bridging, soil clay content, andcarbonates. Additionally, macroorganisms facilitate soil porosity,
infiltration, and aggregate stability by ingestion of soil (Brown etal., 2000). Factors affecting soil aggregation are summarized inFig. 7. Macropores, defined as large continuous openings formedby macroorganisms (e.g., earthworm burrows, old root channels)

Fig. 7. Factors affecting soil aggregation (modified from Bronick andLal, 2005).

VZJ | Advancing Critical Zone Science

also have an important influence on infiltration and subsurface
storm flows as reviewed by Beven and Germann (1982) and forsnowmelt infiltration by Mohammed et al. (2018).
Growing vegetation modifies the soil matrix, affecting soil
hydraulic conductivity and soil water storage. Roots alter the poresize distribution and the connectivity between pores as they push
into the soil matrix, and they also release complex organic compounds into the soil (Bengough, 2012). The continuous networkof branched roots that permeate the soil, with new roots frequently
forming while old ones decay, causes hydrological processes tochange. Root length distribution is a key property that controls
connectivity and preferential flow pathways within the rooting
zone and thus impacts infiltration (Lange et al., 2009).

Hysteresis in the Soil Water Retention Curveand Thermal Effects on Hydraulic Properties
Many LSMs assume that the difference between the soil water
retention behavior between wetting and drying phases in unsaturated soils (hysteresis) can be ignored and that the soil can beconsidered as having one unique soil water retention curve, which
is used to solve the Richards equation (Eq. [3]).
However, hysteresis can play a crucial role in the accurate
description of the flow processes within a soil profile (Glass et al.,
1989; Hanks et al., 1969; Ibrahim and Brutsaert, 1968; Scott et al.,1983). Hysteresis is a process that describes the nonidentical natureof equilibrium soil water content in relationship to matrix potential during the wetting or drying phases. The relationship between
actual soil water content and matric potential can be obtained indesorption, i.e., drying of wet soils, or sorption, i.e., gradual wetting of dry soil. The resulting desorption–sorption curves are
generally not identical because equilibrium soil water content isgreater at a given suction during drying than during wetting. The
relationship of actual water content and matrix potential has been
extensively studied by Haines (1930), Everett (1955), Poulovassilis(1962), Topp (1971), Mualem (1974), Mualem and Dagan (1975),
Parlange (1976), Hogarth et al. (1988), Nimmo (1992), Bachmannand van der Ploeg (2002), Huang et al. (2005), and Mualem andBeriozkin (2009).
Thermal gradients also induce significant changes in theestimated water fluxes because temperature affects soil hydraulic properties (Ben Neriah et al., 2014; Gardner, 1955; Grant andBachmann, 2002; Grant and Salehzadeh, 1996; Hopmans andDane, 1986; Nimmo and Miller, 1986; Parlange et al., 1998; Philipand de Vries, 1957; She and Sleep, 1998). For example, increasingwater temperature decreases water viscosity, causing an increasein hydraulic conductivity (Levy et al., 1989) and thermal swelling of solid particles that change soil pore characteristics and thesolid–liquid interface between soil particles (Gao and Shao, 2015).
The advancement of innovative modeling that includes thehysteretic nature of the soil water retention curve was reviewed
and further developed by Nuth and Laloui (2008), but it has tobe mentioned that for large-scale LSMs the inclusion of hysteretic complexity requires greater computing capability and thep. 29 of 53

knowledge of input parameters from observations and databases
or availability of appropriate PTFs.

Soil Water Repellency
Soil water repellency (SWR) or hydrophobicity reduces theaffinity of soils to infiltrating water such that they resist wetting
for periods ranging from a few seconds to hours or even weeks (e.g.,
Doerr and Thomas, 2000; King, 1981). Additionally, soil water
repellency is spatially and temporally very variable (Regalado andRitter, 2008; Ritsema and Dekker, 1998; Täumer et al., 2005). Soilwater repellency is mostly caused by the coating of the soil particles
by hydrophobic substances, whereby different organic compounds
derived from living or decomposing plants or microorganisms canbe responsible for SWR. Soils below particular vegetation types(such as needle leaf trees), soils with higher soil C content, coarsetextured soils, as well as areas with frequent wildfire are more prone
to SWR than others. A review of factors affecting SWR was given
by Doerr et al. (2000). As mentioned, SWR will reduce the soil infiltration capacity (e.g., Imeson et al., 1992; van Dam et al., 1990) andtherefore will increase overland flow (e.g., Crockford et al., 1991;
McGhie and Posner, 1981; Witter et al., 1991). Topsoil SWR may
cause Hortonian overland flow (runoff) even during precipitation
events with rates much smaller than the saturated hydraulic conductivity. In LSMs where Ks is derived from PTFs and where water
repellency is not included, this may lead to overestimation of infiltration for SWR-prone soils. In some areas, water-repellent layers
underlie highly permeable hydrophilic surface layers, and here, theinfiltrating water may pond above the water-repellent layer; subsequently, the infiltration water can be stored above this hydrophobic
layer and be used for evapotranspiration.
This process can also cause saturated excess overland flow if thepermeable layer becomes fully saturated, can cause lateral water flow
either through structural gaps or along the slope of the hydrophilic

layer, or the water can move downward through the hydrophilic
layer along preferential flow paths (Doerr et al., 2000). A schematic
illustration of the possible hydrological responses caused by top- andsubsoil hydrophobicity is provided in Fig. 8. According to our knowledge, SWR has not yet been implemented in any LSM.

Compaction, Swelling, and Shrinkage
Compaction is the process of reducing the volume of voids in asoil, mainly those filled with air, by packing the soil particles closer
together. It can result from natural processes such as soil overburden or from anthropogenic causes, such as the use of cultivation
machinery or cattle grazing. Compaction is often characterized byan increase in soil bulk density. This has often been considered as
an appropriate independent variable to quantify the decrease in thesoil saturated hydraulic conductivity (Ahuja et al., 1989; Assouline
and Or, 2008; Laliberte et al., 1966; Or et al., 2000) or changes inthe soil hydraulic functions (Ahuja et al., 1998; Assouline, 2006a,
2006b; Stange and Horn, 2005) following compaction. Theseestimates can thus be applied to evaluate the impact of soil compaction on infiltration. Soil swelling during wetting and shrinking
during drying induces dynamic changes in porosity with changing
water content of the soil and changes in the hydraulic properties,
which consequently affect infiltration (Giraldez and Sposito, 1985;
Philip, 1970; Raats and Klute, 1969; Smiles, 1974; Sposito, 1975).
In general, the macroporosity, and to a lesser extent the microporosity, of swelling and shrinking soils is affected by their shrinkage
and swelling behavior (Alaoui et al., 2011), whereby exactly thesevoids in the pore system are highly important for rapid water
infiltration into the soil and the separation between infiltration
and runoff. Electrolyte concentration of the applied water also
can have a significant impact on soil hydraulic properties and on
the infiltration process. The way the soil structure responds toelectrolyte concentration depends on pedogenic processes and theFig. 8. Schematic illustration of possiblehydrological responses of the soil under wettable and water-repellent soils, layer bounds,
spatial variability of soil hydraulic properties, and macropore flow induced by soil
fauna and vegetation. (From Or et al. classnotes with permission.)

VZJ | Advancing Critical Zone Sciencep. 30 of 53

nature of the parent material. For example, a high proportion of
Na+ ions relative to other cations weaken the bonds between soil
particles, decreasing hydraulic conductivity (Frenkel et al., 1978;
McNeal and Coleman, 1966; Quirk, 1994; Rengasamy and Olsson,1991). This process and its impacts on soil physical and chemical
properties has been described in several studies (Assouline andNarkis, 2011; Assouline et al., 2016; Bresler et al., 1982; Greeneand Hairsine, 2004; Jury et al., 1991; Kim and Miller, 1996; Quirk
and Schofield, 1955; Russo, 2005).
Additionally, expansive soils, including peat, can adsorb large
quantities of water during rainfall and therefore reduce surface
runoff. According to the USDA soil classification, clayey soils withclay content >30% (often Vertisols) cover around 320 million ha
globally and are sensitive to swelling and shrinkage (Dinka andLascono, 2012). Several studies have looked at the dynamics ofshrinking and swelling and associated crack changes for the purpose of improving hydrological models (e.g., Arnold et al., 2005;
Bronswijk, 1991; Kishné et al., 2010) but not for LSMs. An extensive review of these models was provided by Adem and Vanapalli(2015). Unfortunately, the shrink–swell properties of Vertisols varyalso as a function of soil properties, climate, topography, vegetation, cropping management, and management practices (Davidson
and Page, 1956; Lin et al., 1998; Thomas et al., 2000; Vaught etal., 2006), which complicates representation in hydrological andland surface models.

Freeze and Thaw
Many soils at higher elevation or latitudes freeze and thaw
seasonally, impacting the soil physical properties and therefore
substantially affecting water movement in the landscape. The
main effect of freeze–thaw cycles (FTCs) on soil properties lies intheir impact on the soil structure, which, as shown above, regulates infiltration and runoff to a large extent (e.g., Chamberlain
and Gow, 1979; Fouli et al., 2013; Qi et al., 2006). Freezing andthawing processes induce uneven stress within the soil, but theconclusions in the literature about the effects on soil structure andwater flow are not unanimous. There are indications that FTCs
decrease soil stability (Edwards, 1991; Kværnø and Øygarden,
2006), whereas Lehrsch (1998), Lehrsch et al. (1991), and Park etal. (2011) observed increasing stability after a few FTCs, while anincreased number of FTCs caused a decrease in soil stability, leading to changes in soil hydraulic parameters over time. On the other
hand, there seems to be more consensus that the effect of FTCs
on clayey soils is much larger than on coarse-textured soils (Bisal
and Nielsen, 1967; Kværnø and Øygarden, 2006). Unger (1991)
additionally stated that FTCs decrease soil bulk density.
It is also known that the hydraulic conductivity of frozen
soil decreases rapidly as the temperatures fall (Williams andBurt, 1974), and some models do take this into account (e.g.,
CLM, Noah-MP, SSiB, SURFEX, CABLE, OLAM-SOIL,
ORCHIDEE). Additionally, even water in the liquid phase isimpacted by temperature changes as the viscosity of the pore water
increases significantly with decreasing soil temperatures (Hillel,

VZJ | Advancing Critical Zone Science1998), leading to lower fluidity and water percolation even before
freezing. Finally, if the freezing front is near the soil surface, ponding is likely to occur at the soil surface after a precipitation event,resulting in runoff because the amount of liquid water-filled pathways has been reduced. Most LSMs take this effect into account.
Also, as the freezing front moves down the soil profile, soil water
will migrate toward the freezing front, leaving a drier soil behind,resulting in a larger matric potential gradient pulling the water
toward the freezing front (Jame, 1977). Even though some LSMs
account for the direct impact of freezing on Ks, temporal changes
in the hydraulic parameters due to structural changes induced by
FTCs have not been implemented yet. This might be problematic
for regions where FTCs might become more frequent under futureclimate change, as stated by Eigenbrod (1996).

Impermeable Soil Layers
Impermeable layers, or more precisely soil horizons withextremely low saturated hydraulic conductivity, frequently occur
in natural or managed soils. Often these layers are denoted as
hardpans, hard layers, or compacted horizons located at either thesurface or subsurface (Busscher, 2011). These layers can be causedby traffic, tillage practices, trampling of livestock, or soil-forming
properties that result in layers with high density or cemented soil
particles (Hamza and Anderson, 2005; Silva et al., 2000). For
example, the extent of compacted soil is estimated worldwide at
68 million ha of land from vehicular traffic alone (Flowers andLal, 1998). Some of these compacted or extremely dense soil layers
are relatively thin and are therefore often neglected in soil maps
at coarser scales. Additionally, changes in the saturated hydraulicconductivity due to soil compaction is often not accounted for in
PTFs if bulk density is not used for the prediction of the hydraulicparameters (Van Looy et al., 2017). Nevertheless, these layers are ofutmost importance because they control the infiltration of water
into the soil and its redistribution to greater depth. In general, thepresence of impermeable soil layers will lead to the same hydrological response as shown for the hydrophobic layers depicted in Fig. 8,
generating more overland flow (runoff) and subsurface flow events.
Impermeable layers complicate the naturally occurring soil
vertical heterogeneity, where generally successive distinct layers ofsoil with different hydraulic properties occur. Several studies haveproposed solutions for infiltration into layered soil systems (Childs
and Bybordi, 1969; Colman and Bodman, 1945; Hanks and Bowers,
1962; Miller and Gardner, 1962; Philip, 1967; Raats, 1973; Warrick
and Yeh, 1990; Zaslavsky, 1964). Chu and Marino (2005) presented
a solution for determining ponding conditions and simulating
infiltration into a layered soil profile based on the Green and Ampt
approach for unsteady rainfall. Beven (1984) and Selker et al. (1999)
also extended the Green and Ampt model for infiltration into soil
profiles where pore size varied with depth. A review of the applications of the Green and Ampt model to vertically heterogeneousconditions was provided by Kale and Sahoo (2011).
A special case of a layered soil profile occurs when a seal layeror crust develops on the soil surface, resulting from the destructive
p. 31 of 53

action of raindrop impacts on the soil, which alters the soil structure and soil hydraulic properties, especially Ks. This process andits impacts on physical and chemical properties was described by
Quirk and Schofield (1955), Bresler et al. (1982), Jury et al. (1991),
Kim and Miller (1996), Greene and Hairsine (2004), Russo (2005),
and Assouline et al. (2015). A review of the approaches proposed
to model infiltration into sealed (or crusted) soils can be foundin Mualem and Assouline (1992), Mualem and Assouline (1996),
and Assouline (2004). The direct effect of the presence of the
impeding seal layer at the soil surface is to reduce ponding time
and the infiltration rate during rainfall (Römkens et al., 1986a,
1986b). Hillel and Gardner (1969, 1970) first addressed the problem of infiltration in the case of sealed soils. They presumed that asealed soil can be modeled as a uniform soil profile capped with asaturated thin layer of low permeability with constant prescribed
physical properties such as the saturated hydraulic conductivity.
Their simplified solution was based on the Green and Ampt model,
assuming a constant water content (or suction) at the interface
between the seal and the soil beneath. It was further applied indifferent studies (Ahuja, 1974, 1983; Moore, 1981a; Parlange et al.,1984). Variations and extensions of this basic approach included
the simulation of infiltration with time-dependent seal hydraulic conductivity functions (Ahuja, 1983; Brakensiek and Rawls,
1983; Chu et al., 1986; Farrell and Larson, 1972; Moore, 1981b;
Vandervaere et al., 1998; Whisler et al., 1979). An additional
conceptual model, based on the model of Corradini et al. (1997),
was suggested by Smith et al. (1999). Römkens and Prasad (1992)
applied the solution of Prasad and Römkens (1982) based on the
spectral series approach to solve the infiltration equation in soils
topped by a constant or transient crust.
Note that in all these studies the hydraulic properties of theseal layer were arbitrarily chosen. Mualem and Assouline (1989) as
well as Baumhardt et al. (1990) have addressed the problem of infiltration into sealed and sealing soils by attributing to the seal layer
hydraulic functions that evolved from those of the undisturbed
soil and that were related to the specific rainfall kinetic energy andintensity involved in the seal formation. The impact of soil surface
sealing on infiltration is illustrated in Fig. 9, which depicts that

Fig. 9. The impact of soil surface sealing on infiltration.

VZJ | Advancing Critical Zone Sciencesoil surface sealing reduces the ponding time and the infiltration
rates into the soil profile, including the final quasi-steady rate. As
a result, much more runoff is formed by a given rainfall event when
the soil surface sealing is accounted for.

Instability of Different Flow Regimes
Wetting front instability occurring under certain flow regimes
can also significantly affect the infiltration process (DiCarlo, 2004;
Jury et al., 2003; Or, 2008; Parlange and Hill, 1976; Philip, 1975;
Raats, 1973). Wetting front instability refers to a splitting up of theinfiltration front into several fingers along which water is rapidlytransported downward. Because a part of the soil pore volume is
bypassed by the infiltration through fingers, wetting front instability leads to considerably deeper infiltration than in the case of stable
wetting fronts. Raats (1973) explained that an increase of soil waterpressure with depth above the wetting front in general leads to instabilities of the wetting front. Entrapment of air, the presence of layerswith higher water-entry values, water repellency, but also the reversalof pressure gradients during redistribution just after infiltration at
the soil surface has ceased can cause such an increase in pressure
above the wetting front that leads to unstable wetting fronts (Wanget al., 1998, 2003a, 2003b). Another process, which originates at thepore scale and which can explain the persistence of individual fingers due to pressure increase or pressure overshoot above the wetting
front of a single finger, is the dynamic pressure–water content relation that results from a rapid filling of larger pores and a subsequent
redistribution (DiCarlo, 2013).

Solution of Numerical Issues
Rainfall of different intensity also affects infiltration depths
and runoff ratios (Frauenfeld and Truman, 2004). For example,
varying intensity rainfall simulations yield larger runoff ratios andpeak runoff rates than uniform rainfall simulations (Dunkerley,2012). Using Horton equations, predicted runoff rates were significantly improved during intra-event time variation of fluctuating
rainfall simulations (Dunkerley, 2017).
Conventional solution methods to the highly nonlinear
Richards one-dimensional partial differential equation used in
LSMs inevitably lead to numerical and accuracy issues, which
impact their hydrological performance. Land surface models
may consider the implementation of alternative, one-dimensional,
unsaturated zone flow solution methods (such as those provided by
Ogden et al., 2015). The Ogden Soil Moisture Velocity Equation(SMVE) approach uses the hodograph method to transform
Richards’ equation into a differential equation for a velocity anduses a discretization of the resulting equation in the form of “bins”
containing values of the water content. The scheme is computationally efficient, although the explicit time steps are limited by stability
considerations because there are no convergence limits as imposedby implicit schemes. Ogden et al. (2015) considered the transportof three regimes of soil water in detail, namely infiltration, wetting
fronts disconnected from the surface, and groundwater recharge.
The SMVE method offers accuracy comparable to, or in some cases
p. 32 of 53

exceeding, that of the numerical solution of the Richards partial
differential equation method but without the numerical complexity and in a form that is robust, continuous, and suitable for use inmodels of coupled climate and hydrology at a range of scales.

66Summary and Conclusion
Infiltration processes are at the core of LSMs, representing
the complex and highly dynamic coupling between precipitationand land surface properties where soil, vegetation, initial soil conditions, and topography interact. Although the formulation of
infiltration representation for the soil profile scale is well established and tested by the soil physical community, there are still
issues that need resolving with regard to the parameterization of
infiltration in LSMs. In particular, the extension of the concepts
to the catchment and global grid cell scales remains challenging
and is in some cases tentative, and with various different solutions
that are currently in use. In this perspective, we reviewed and analyzed the different approaches used in current LSMs to predict soil
infiltration processes. Specific attention was given to the underlying physical principles and concepts used to predict infiltration
at the point and grid scales and the approaches used to describe
spatial heterogeneity and upscaling of key parameters controllingthe infiltration process in LSMs. We identified several topics andprocesses that warrant further attention in advancing the prediction of infiltration processes.
First, there is the prediction of Ks, a key parameter in describing infiltration. Currently, Ks estimates in LSMs are derived from
PTFs that are typically based on the textural composition of soils
but do not consider the impact of soil structure on the infiltration process in general. Recently, Rahmati et al. (2018) publisheda global database of infiltration measurements that clearly showsthat Ks derived from field experiments cannot be predicted fromsoil texture alone. Therefore, research needs to be directed toward
the development of PTFs that consider the effect of structuralproperties on K, e.g., using land use and tillage treatment as proxies(Jorda et al., 2015). This might be even more important as a recent
study of Hirmas et al. (2018) indicated that drier climates induce
the formation of greater soil macroporosity than do more humid
ones and that such climate-induced changes occur across shorter
timescales than have previously been considered. Translation of
the effects of these different, largely exogenic, processes to timevarying hydraulic properties (currently hydraulic properties in
LSMs are kept constant in time) is one of the greatest challenges
in current land surface modeling. For example, the increase inhigh-frequency rainfall events under future climate conditions
will make crust formation for certain soil types more likely, which
will cause a decrease in infiltration and increase in surface runoff.
Ignoring these aspects will add further uncertainties to predictions of future land–atmosphere interactions. This issue needs tobe addressed urgently and in a coherent fashion whereby othersoil properties (e.g., thermal properties) and vegetation parametersthat depend on, or affect, soil properties (such as rooting depth)

VZJ | Advancing Critical Zone Scienceare changed concurrently. Additionally, numerical simulations areneeded to quantify the effect of Ks estimates considering the roleof soil structure on the energy, water, and matter cycles.
Second, due to the availability of spatially highly resolved
soil map information at the global scale with a spatial resolutionof 250 m or even less, quantification of the subgrid variability is
now within reach. The use of this information in combination
with PTFs allows direct estimation of b, a lumped parameter usedin several LSMs to describe the spatial variability of infiltration
capacity. In addition, this highly resolved spatial information canbe used to derive effective soil hydraulic parameters such as theMualem–van Genuchten parameters, which are used in the solution of the Richards equation. The increasing availability of highly
resolved spatial data poses questions about how to effectively andefficiently represent subgrid soil and landscape information in
LSMs. The strengths and weaknesses as well as the validity andapplicability of the methods presented here with respect to modeling land surface processes at the continental and global scales still
has not been addressed.
The requirement to correctly represent Hortonian infiltration
but also redistribution processes of water in the subsurface (e.g.,
due to root water uptake or capillary rise) is best fulfilled by usinga Richards-equation-based approach. Stochastic analyses of waterflow in spatially heterogeneous soil fields (Mantoglou and Gelhar,
1987a, 1987b; Vereecken et al., 2007) have shown that the upscaled
Richards equation at the field or larger scale has a form similar
to the local-scale equation. However, the spatial variability of soil
hydraulic properties introduces a hysteretic behavior of the larger
scale system, as the effective hydraulic conductivity is a function of
the hydraulic gradient and of its history reflecting nonequilibrium
conditions. How far this nonequilibrium behavior is relevant for
grid-scale infiltration processes needs to be further studied. The
definition of the effective parameters in the upscaled Richards
equation, however, requires detailed knowledge of the spatial
statistics of local-scale hydraulic parameters. The availability ofhighly resolved soil maps in combination with PTFs opens up newopportunities to define subgrid variability of hydraulic parameters
and thus to quantify effective hydraulic parameters at the scaleof LSMs. Also for heterogeneous porous media, the solution of
the Richards equation for an infiltration problem remains stable(Egorov et al., 2003). To represent the impact of soil structure,
macropores, cracks, or other well-connected structures on water
infiltration in soils, several modifications ranging from changing
the typically used unimodal pore size distribution to a dual or
multimodal pore size distribution to introducing an extra flow
equation that represents the infiltration in the macropore pore network and that is coupled with the flow equation in the soil matrix
have been proposed (see reviews of Jarvis [2007] and Šimůnek et
al. [2003]). These well-connected and highly conductive structures
could also be represented in three-dimensional Richards models(e.g., Vogel et al., 2006). But even for such media, the solution of
the Richards equation leads eventually to relatively stable infiltration profiles that could be represented fairly well by an upscaledp. 33 of 53

Richards equation with effective parameters (e.g., Schlüter et al.,2012). It must be noted, though, that the local water fluxes above
the infiltration front can be very heterogeneous (but the wetting
front is relatively homogeneous). However, the Richards equation
cannot reproduce unstable infiltration fronts that are observedat the local scale as a consequence of pore-scale dynamic effects.
Phenomena like finger development in gravity-dominated flow,which can have an important impact on the vertical distributionof the infiltrated water and how it varies with the infiltration rate
at the soil surface are therefore not represented by the Richards
equation. Several approaches to account for these dynamic andnonequilibrium processes by adding additional terms to the continuum Richards equation have been proposed (Cueto-Felgueroso
and Juanes, 2009; DiCarlo, 2013; Eliassi and Glass, 2002).
Although these approaches describe experimentally observed
non-uniform infiltration fronts, it still requires further investigation how the upscaled unstable infiltration can be described bya continuum model and what its consequences are for the water
distribution during an infiltration event at the LSM grid scale.
Also, correct representation of Hortonian infiltration requiresconsideration of the vertical heterogeneity of soil hydraulicparameters, a vertically variable discretization with the finest discretization near the surface, and the use of a pressure-head-based
Richards equation. Recently, LSMs such as Parflow/CLM, ISBA,
CLSM, and OLAM-SOIL have developed approaches that allow
these requirements to be fulfilled. Introducing similar approaches
in other LSMs automatically avoids the need to define a maximum
infiltration capacity of soils, leading to a more physically consistentdescription of infiltration.
We observed a disparity between the approaches used at
the field or small catchment scale and the approaches applied at
the grid scale in LSMs. This, of course, is a result of having two
different scientific communities working predominantly at different spatial scales. The soil physics community mainly focuseson the field scale and typically uses semi-analytical solutions
or full implementations of the Richards equation to explicitly
solve for infiltration flux. This in general requires fine vertical discretization (approximately millimeters near the surface
boundary) and short time steps (seconds) to calculate infiltration fluxes. Modeling of infiltration in LSMs is performed ata much larger scale, which usually does not allow a fine spatial
and temporal discretization to keep the models computationally efficient. The majority of the LSM community has therefore
taken the approach of parameterizing the infiltration process
at the land surface and using the Richards equation, mainly inthe diffusive form, to redistribute infiltrated water in the soil
profile. The common basis for both approaches is the Richards
equation, even though for different reasons. It is, however, thegoal of this review to foster the cooperation and the exchangeof ideas between the two communities. As a first step, the workof Montzka et al. (2017) provides a global concept of subgrid
variability of soil hydraulic properties along with the methods ofsimilarity scaling. As the need increases to account for subgrid

VZJ | Advancing Critical Zone Sciencevariability in LSMs, these above-mentioned methods provide
options for incorporating this uncertainty.
Furthermore, there is a large diversity among the analyzed
LSMs in estimating key properties such as soil moisture capacity
and in the treatment of the heterogeneity of soil moisture at the gridscale. In the case of soil moisture heterogeneity, three mathematicalformulations have been used: (i) reflection power distribution functions, (ii) g distributions, and (iii) exponential distributions with avariable number of parameters (two or three). In some cases, LSMs
also use different approaches to derive the saturated fraction (Fsat)of a grid cell, which is used to partition between Dunne saturation excess and Hortonian infiltration excess. Besides differencesin concepts used to formulate the saturated fraction, there is a large
divergence in the way the saturated fraction of the land surface
within a pixel, Fsat, is being parameterized.
Differences include whether or not the groundwater depth is
explicitly simulated, and if so, how, and the treatment of the storage capacity of the soil between the land surface and the bedrockor groundwater table depth.
In addition, our analysis showed that basic soil informationthat is used to obtain spatial coverage of key soil hydraulic properties differs strongly among LSMs but also among the various
versions of one single LSM. The impact of using different spatial
soil maps combined with the wide range of approaches used toestimate Ks and other soil hydraulic properties is not yet known.
Further research is needed in this direction to quantify the impactof this input variability.
Also, many LSMs use a prescribed parameterization of the
maximum infiltration capacity to partition precipitation between
infiltrable water and runoff (exceptions being Parflow/CLM,
ISBA-SURFEX, ORCHIDEE, CLSM, and OLAM-SOIL). These
approaches have been heavily tuned by each LSM to ensure thatthey fit with runoff observations. The lack of a general framework
for this central hydrological process leaves a serious gap in the present LSM parameterizations and hinders simple and transparent
updating of soil information when it becomes available.
Moreover, re-infiltration, called run-on, is ignored in mostLSMs; runoff production that occurs at sites where the infiltration
capacity is exceeded may re-infiltrate in the grid cell due to soil
and land surface heterogeneity so that not all of the runoff thatis generated at a grid cell needs to be routed out of the cell or toa receiving water body. A classic example is the run-on in vegetation patches or bands (strips) in semiarid regions (Assouline et al.,2015). Roots can increase the local infiltration capacity so thatrunoff from sealed, unvegetated areas can infiltrate in vegetatedareas (Nimmo et al., 2009). In addition, run-on leads to a scaling
behavior of rainfall–runoff relations, with generally less runoff
produced at a larger scale than what would be derived from smallerscale rainfall–runoff relations. A crucial property that defines
the rainfall–runoff relations at larger scales is the connectivity ofregions that generate runoff (Herbst et al., 2006).
Finally, several processes that control infiltration and thus
impact the soil water balance, and ultimately the energy balancep. 34 of 53

and related land–atmosphere interactions, at the grid cell scalerequire more attention. This includes the role of vegetation in
the infiltration process, the role of the runoff–run-on process at
the grid cell scale, and the dynamics of soil structural properties.
Correct representation of the runoff–run-on process will needspatially distributed information about parameters controlling
Hortonian infiltration excess generation and the formulation andparameterization of redistribution mechanisms within the gridcell. The role of vegetation is related to the effect it exerts on the
structural status of the vadose zone, leading to soil properties thatare changing over time. In addition, changes in land use and management may affect the structural status of the vadose zone andthus affect water infiltration into soils. One way forward would
be to develop PTFs that consider time-dependent soil properties.
This, however, requires a basic understanding of how vegetationand management practices change soil hydraulic properties of soils,
suggesting the need for greater integration of soil physics, plantscience, and land management. There is increasing evidence thatspatial variability in water infiltration may also be attributed to
dynamics of vegetation-driven spatial heterogeneity (Archer et al.,2012, 2013; Puigdefábregas, 2005) leading to increased infiltration capacity of soils. These processes that may lead to a decreasein Hortonian infiltration have not yet been introduced in LSMs.
Currently, activities have been initiated between ISMC
SoilMIP, and the Global Energ y and Water Exchanges(GEWEX) project (https://soil-modeling.org/activities/events/
the-gewex-soilwat-initiative-first-planning-workshop-for-scopeand-interactions-advancing-integration-of-soil-and-subsurfaceprocesses-in-climate-models) to advance the implementation ofhigh-quality soil information and the description of soil processes
in LSMs. These improved LSMs, in turn, will feed into Earthsystem models for global prediction and closure of water, energy,
and carbon budgets. This review is a part of this initiative and oneof the first outcomes of this joint activity based on a workshop
held in Leipzig, Germany, in 2016. Further activities are presently
ongoing, such as the analysis of the effect of incorporating soil
structure on the soil water balance of the terrestrial system, andnew ones are being initiated and developed.

List of Symbols
a
ar
as
aSSiB
A
ACi
Af
Agf
A1
b
bBC
br
bSSiB
B

fraction of surface runoff
fitting parameter in CLSM
fitting parameter for each grid box to estimate Fsat (JULES)
constant in Eq. [A54]fraction of an area for which the soil moisture capacity is £I
steady-state infiltration at the ith location [L T−1]
fitting parameter
grid box fraction of soil water capacity
fitting parameter [L T−1]
parameter that reflects the grid cell heterogeneity
parameter in the Brooks–Corey equation (Clapp and Hornberger, 1978)
fitting parameter in CLSM
constant in Eq. [A54]potential infiltration rate shape parameter—a measure of the spatial
variability of the potential infiltration rate, defined as the maximum
infiltration rate of each point when the surface is ponded (Liang and Xie,
2001)

VZJ | Advancing Critical Zone Science

B¢
BTW
c
cmax
cs
cSSiB
C
Cice
Cv
do
dx
ád2ñ
Db
Dbmax
Dx
Dsm
Dt
Dq
f
f
f*
fc
fd
fh20sfc
fH
fm
fmm
fsno
fs(s)
fsc,i
funf,i
f0
F(c)
Ffrz
Fmax
Fs(z*)
Fsat
hsat
h
hA
hc
hf
ho(t)
hS
H
i
i*(t)
ic
ip
I
I*
Ic,t
áImñ
Imax
Ir
Iroff
Iunf,iJ
Jw

parameter related to the Brooks–Corey exponent of the water retention
characteristic
exponent of the tension water capacity distribution curve
water storage capacity at a certain location [L]maximum water storage capacity [L]fitting parameter for each grid box to estimate Fsat (JULES)
constant in Eq. [A54]fraction of an area for which the potential infiltration rate is less than orequal to f (Liang and Xie, 2001)
ice impedance factor
tunable parameter
maximum local water deficit [L]length of the grid cell (km)depth of the root zone [L]spatially averaged soil moisture storage deficit
maximum soil water storage [L]soil water storage [L]local soil moisture deficit [L]mean water storage deficit [L]the amount of water stored between surface and groundwater [L]
pdf of the respective variables
point potential infiltration rate [L T−1]
infiltrability of the first soil layer [L T−1]
final infiltration capacity [L T−1]
decay factor [L−1]fraction of the area that is inundated
Hortonian infiltration capacity at time t [L T−1]
maximum potential infiltration rate [L T−1]
average potential infiltration rate [L T−1]
fraction of the grid cell covered by snow
gamma function of the mean surface layer point soil water saturation
fraction of soil carbon in the ith layer
pdf representing the spatial variability of the local maximum infiltration
rate in a grid cell
initial infiltration capacity [L T−1]
standard reflected power cumulative distribution function of the spatial
variation of the storage capacity cfractional impermeable area as a function of soil ice content of the surface
soil layer cumulative distribution of mean infiltrability (ORCHIDEE)
fraction of pixels in a grid cell whose topographic index is ³ the grid cellmean topographic index
cumulative distribution function of the scaled infiltration capacity (Eq.
[46])
saturated fraction of the catchment or grid cell in a LSM
pressure head at air-entry value [L]matric potential [L]minimum (i.e., for evaporation) pressure head at the soil surface allowed
under the prevailing soil conditions [L]matric potential at air entry [L]capillary pressure at the wetting front [L],matric potential at the soil surface at time t [L]maximum (i.e., for infiltration) pressure heads at the soil surface allowed
under the prevailing soil conditions [L]total water head [L]point soil moisture capacity [L]scaled infiltration rate [L T−1]
infiltration capacity [M L2 T−1]
maximum potential rate of infiltration or evaporation under the current
atmospheric conditions [L T−1]
spatially averaged actual infiltration [L T−1]
scaled cumulative infiltration [L]cumulative infiltration capacity at time t [L]mean maximum infiltration rate over the grid cellmaximum infiltration rate [M L2 T−1 or L T−1] or maximum soil moisturecapacity [L]infiltration rate for the whole soil column [L T−1]
normalized runoff in SSiB [L T−1]
local maximum infiltration rates for unfrozen soil [L T−1]
Darcian water flux [L T−1]
precipitation rate [L T−1]

p. 35 of 53

Kdt
Kiint

constant parameter set equal to 3.0
average of hydraulic conductivity at the wetting front and the deepest
saturated node [L T−1]
intscaled Kiint [L T−1]
Ki *hydraulic conductivity of the soil surface layer [L T−1]
Koreference K value (2 ´10−6 m s−1) [L T−1]
KrefKs
saturated hydraulic conductivity [L T−1]
Ks,min,i saturated hydraulic conductivity of the mineral phase of the ith layer [LT−1]
saturated hydraulic conductivity of the organic soil carbon of the ith layerKs,sc,i[L T−1]ltortuosity index in the Mualem modelmfitting parameter in the van Genuchten modelmpore size distribution index in van Genuchten modelsurface layer excess [L]Msefitting parameter in the Kostiakov equation
nvg
Pdrop
precipitation reaching surface after canopy interceptrion in SSiB [L T−1]
snowmelt [L, M L2 T−1, or L T−1]
PM
rainfall rate [L T−1]
PRPT
throughfall precipitation [L, M L2 T−1, or L T−1]
flux incident at the soil surface [M L2 T−1]
Px
infiltration capacity of the soil [L T−1]
qcap
q(t)
infiltration flux [L T−1]
q*(t)scaled infiltration rate [L T−1]
Dunne saturation excess runoff [L T−1]
qD
qevapo.soil water loss due to evaporation [L T−1]
final (constant) infiltration rate [L T−1]
qf
qH
Hortonian excess infiltration runoff [L T−1]
initial infiltration rate [L T−1]qiqin,soil infiltration rate of water into the soil [L T−1]
qin,surface surface moisture flux remaining after surface runoff has been removed [LT−1]qliq,omoisture input into the grid cell [L T−1] and is the sum of liquid
precipitation reaching the surface and snowmelt
recharge rate [L T−1]
qresurface runoff (ML2T−1) or [L T−1]
qsqudownslope subsurface flow rate per unit contour width [L T−1]
Q(t)cumulative infiltration at time t [L]input of water (rainfall, snowfall, dew) [L T−1]
Qwat
r(t)
water supply rate [L T−1]roamount of water that needs to infiltrate (throughfall, snowmelt, andponded water) [L T−1]
R(i)mean infiltration excess runoff [L T−1]swater content expressed as saturation
grid mean surface layer soil water saturation
ásñS
grid cell mean storage
effective saturationSe
surface gradientSg
Slmean slope in the grid cells
maximum threshold slope
Slmax
Smax
maximum water depth across the basin [L]minimum storage below which there is no surface saturation [L]Sosoil sorptivity at the ith location[L T−2]
Sop,iSW
equilibrium root zone wetness index
minimum value of the soil wetness index distribution
SWo
SWI
soil water index
local terrain surface slope
Sxl*ttime [T]t*scaled time [T]local surface topographic slope
tanb
time needed to infiltrate a certain amount of water [T]tftime to ponding [T]tp
Ttransmissivity [L T−1]
transmissivity of the soil profile [L T−1]
To
T(zÑ) transmissivity in Eq. [A38]v
depth averaged velocity vector [L T−1]
vK

parameter describing the decrease of saturated hydraulic conductivity withdepth

VZJ | Advancing Critical Zone Science

wact
wmin
wmax
W
Wsat
x

X
zbot
z(i)
zl
zpdm
zwt
zÑ
Z
Zg
Zm
a AC
aG
aGA
aK
a opt
as
a sc
a Sop
aT
awf
b
bG
bK
bs
bT
g
gf
g1,2,3
dt
e
q
qi
q ice
ql
qr
q rz
qs
q wilt
k
kp
l
lcrit
lm
lp
lr
ls
l sc
r
sKs

subgrid water content that corresponds to the fractional saturation
minimal local subgrid soil water capacity
maximum local soil water capacity
vertically integrated soil water content (over the first 50 cm of the soil
profile) [L]vertically integrated saturated soil water content (over the first 50 cm of
the soil profile) [L]mean catchment value of ln(a/tanb)
dimensionless time
middle point of the bottom layer (1.5 m) [L]
thickness of the ith soil layer [L]local water table depth [L]depth over which soil moisture is considered [L]water table depth [L]mean water table depth [L]soil profile depth to bedrock [L]ground surface elevation [L]maximum value of profile depth to bedrock [L]Miller–Miller scaling factor for final infiltration capacity
parameter in the Gamma distribution
parameter proportional to the matric potential at the wetting front
decay constant in Kostiakov equation
single Miller–Miller scaling factor for infiltration
shape parameter in the spatial distribution function of soil wetness
Miller–Miller similarity scaling factor
Miller–Miller scaling factor for sorptivity
upstream area that contributes flow through a unit contour positioned at
the point or specific catchment area
product of the absolute value of wetting front suction head and thedifference between saturated water content and initial water content at the
beginning of the infiltration event
parameter depending on grid resolution
parameter in the Gamma distribution
fitting parameter in the Kostiakov equation
correction parameter for Ks
slope parameterparameter related to topography
fraction of surface runoff
fitting parametersmodel time step [L]dimensionless error
volumetric soil water content
initial water content
the content of ice
the content of liquid water
residual moisture content
mean diagnosed root zone moisture content
saturated soil water content
water content at wilting point
scaling factor needed to redistribute the GCM grid scale precipitation overthe scale of precipitation events
related to the standard deviation of the subgrid slope
wetness index
critical topographic index
mean topographic or wetness index
dimensionless pore size distribution index
characteristic length of Miller–Miller scaling reference soil [L]subgrid heterogeneity of soil
characteristic length of a Miller-Miller similar porous medium [L]correlation coefficient between the terrain slope and the soil water content
standard deviation of saturated hydraulic conductivity [L T−1]

p. 36 of 53

smax
smin
sor
sslopesSxl
sqtæç K int
F ( K iint ) = 1- expççç- iintççè K imaximum standard deviation of the orography
minimum standard deviation of the orography
standard deviation of the orography
standard deviation of the subgrid slope variability
standard deviation of local surface slopes
standard deviation of soil water content
time scale of transfer of surface layer moisture into the root zone [T]66Appendix A1

The Probability Distributed Model of Moore (1985)
In the PDM-based scheme (Moore, 1985) for calculating
Dunne runoff, Fsat was described asb /(b +1)
æS - S0 ö÷÷÷Fsat = 1-ççç1èç S max - S0 ÷ø[A1]where S is the grid cell mean water storage [L], S 0 is the minimum water storage below which there is no water saturation at thesurface [L], Smax is the maximum possible grid cell water storage[L], and b is a shape parameter proposed by Moore (1985) thatreflects the heterogeneity in a lumped manner (see also Eq. [40]).
Parameters b and S 0 were obtained from model calibration using
catchment data, and Smax was obtained from available data andcalculated from
S max = q sat z pdm[A2]

where zpdm [L] is the soil depth throughout which the soil watercontent is considered for PDM modeling. Clark and Gedney(2008) assumed zpdm to be 1 m.

66Appendix A2

Description of Maximum Infiltration Rate
in Different Land Surface Models
Note that only those models using the maximum infiltration
rate in their concept are described here and that the last number ofthe label (e.g., A2.3) refers to the number given in Table 2 to facilitate reading. Additionally, Imax depends on Fsat for those models
relying on the Fsat approach.

A2.3 ORCHIDEE
ORCHIDEE (Ducharne et al., 2017) includes a subgrid distribution of infiltration, which reduces the effective infiltration
rate into each successive layer of the wetting front. In practice,the mean infiltrability of a layer over the grid cell is spatially distributed using an exponential pdf, then compared locally to theamount of water that needs to infiltrate (called r0 and comprised ofthroughfall, potentially increased by snowmelt and ponded water).
As a result, infiltration-excess runoff is produced over the fraction of the grid cell where r 0 [L T−1] is larger than the local Kiint[L T−1] defined by the exponential distribution of mean saturated
hydraulic conductivity áKiintñ, applying the following cumulative
distribution function:

VZJ | Advancing Critical Zone Science

ö÷
÷÷
÷÷
÷ø

[A3]

The Kiint value is calculated as the average of the actual hydraulicconductivity at the wetting front and the deepest saturated node.
A spatial integration is conducted for each soil layer that becomessaturated when the wetting front propagates, giving the mean
infiltration excess runoff, Re, i, produced from the saturation ofeach ith soil layer:éæ r ö÷ù
Re, i = r0 - K iint êê1- expççç- 0int ÷÷úúçè K i ø÷ú
ëê
û

[A4]

By reducing the effective conductivity,éæ r ö÷ù
K iint * = K iint êê1- expççç- 0int ÷÷úúçè K i ø÷úëêûcompared with the uniform case (Kiint* = Kiint), this subgrid distribution increases surface runoff, given by the sum of Ri from all
the layers saturated during the time step. The model also considers
the mean slope of the grid cell, with a reinfiltration of excess water
possible only at low slopes. This subgrid distribution can be seen as
opposite to the parameterization of Warrilow et al. (1986) because
the actual hydraulic conductivity, K, rather than the precipitation
rate is spatially distributed within the grid cells.

A2.4 CLSM
The Catchment Land Surface Model (CLSM) (Koster et al.,2000), the land model component of the NASA Goddard Earth
Observing System (GEOS-5) coupled Earth system model, does
not impose a priori a maximum infiltration rate in its formulation. The amount of water that can infiltrate within a certain
time at the catchment scale is a function of the model’s dynamically varying spatial moisture fields. Infiltration in CLSM isconsidered here in two steps: (i) precipitation throughfall into
the near-surface soil layer (2 or 5 cm), and (ii) the subsequenttransfer of this soil water into the root zone. It is important torecognize that CLSM is designed to emphasize a descriptionof horizontal moisture variability that is linked to the simulation of a spatially variable dynamic water table depth. This isdiscussed in more detail above. In effect, the land surface areain CLSM is divided into distinct (and dynamically changing)hydrological regimes. Regarding the throughfall into the nearsurface layer, all rainwater runs off the surface in the “saturated
fraction” regime, effectively as Dunne runoff and without
infiltration. In recent versions of CLSM, the other two regimes(the “subsaturated-but-transpiring” and the “wilting” regimes)
allow all precipitation water in a given time step to infiltrate
and thereby increase the surface soil moisture, although if thelayer becomes fully saturated, the excess does run off the surface,effectively as Hortonian runoff.
The transfer of surface layer moisture into the root zone iscontrolled by a time scale, t [T], computed asp. 37 of 53

t=

ar

3

( q rz + br M se )[A5]

where ar and br are fitted parameters, q rz is the mean diagnosed
root zone moisture content [L3 L−3], and Mse is the surface layer
excess [L] (see above). With this timescale defined, the water transferred from the surface layer to the root zone, DMse [L], isD M se =- M se D t
t[A6]

The empirical equation for the timescale t was fitted to resultsfrom high-resolution (1-cm) solutions of the vertical one-dimensional Richards equation, conducted off-line prior to runningclimate or land surface simulations. The simulations behind these
off-line solutions used a comprehensive set of values for the CLSM’swater prognostic variables (see Ducharne et al., 2000) appropriately
downscaled to 1-cm vertical resolution, and a comprehensive set ofsoil classes parameterized by the Campbell (1974) equations, withcorresponding hydraulic parameters based on lookup tables or using
the PTFs of Wösten et al. (2001) (De Lannoy et al., 2014).
In CLSM, the catchment-scale infiltration rate decreases by
two mechanisms in which the actual groundwater level is crucial:(i) the unsaturated area into which rainfall can infiltrate at thecatchment scale decreases for rising water levels, i.e., higher areal
fractions of the saturated regime, and (ii) the hydraulic gradient
between surface and root zone in the unsaturated area decreases
when the root zone fills up due to infiltration and rising waterlevels. The combination of both mechanisms results in a dynamic
prediction of catchment-scale maximum infiltration rates into thesurface layer that range from high values (larger than Ks) underdeep water level conditions to values that drop below Ks undershallow water level conditions.

A2.5 ISBA-SURFEXIn ISBA-SURFEX (Decharme and Douville, 2006), the local
maximum infiltration rates for unfrozen soil are given byéb h æ qö ù
I unf , i = K s, i êê BC c ççç -1÷÷÷ +1úú
[A7]
÷ø ú
êë D z çè q s
û
where hc is the matric potential at air entry [L], q s is the soil
porosity or saturated water content [L3 L−3], bBC is the pore size
distribution index from the Brooks–Corey equation (dimensionless) (see Eq. [4]), Dz is the top layer soil thickness of 0.1 m, andKs,i is the saturated hydraulic conductivity at the ith location[L T−1]. The variable Iunf,i comes from the equation presented in
the study by Abramopoulos et al. (1988) for calculating infiltration and evapotranspiration in global climate models, whereby the
maximum infiltration was defined aséæ q -q s ö÷ ùú
ê dh
[A8]
I max = K s, i ê iççç i÷÷ +1ú
ê d q i q i =q s çè z i ÷ø úëûwhere i refers to the top soil layer (i = 1), q i is the volumetric soil
water content of the top soil layer [L3 L−3], and zi is the thickness

VZJ | Advancing Critical Zone Science

of this layer [L]. The mean maximum infiltration rate, Imax, is used
to calculate the surface runoff generated by Hortonian overland
flow asqHH =é 1 ¥ù [A9]m ê ò ò ( Px , i - I unf , i ) f Px , i ( Px , i )d Px , i f unf , i ( I unf , i )d I unf , i úëê 0 f *ûú
where Px is the incident flux reaching the soil surface and

æ -I
ö
1
f unf , i ( I unf , i ) = expççç max, i ÷÷÷çè I max ÷ø
I

[A10]

where Imax is the mean maximum infiltration rate across the gridcell [L T−1]. Here funf,i(Iunf,i) represents the spatial variability of
the local maximum infiltration rate in a grid cell.
The parameter m in Eq. [58] is given by
m = 1- exp(-b Px )[A11]where b is a parameter depending on the grid resolution according tob = 0.2 + 0.5exp(- 0.01dx )[A12]where dx is the length of the grid cell in kilometers.
Combining all equations, we obtain
qHH =

Px

1 + I unf ( m Px )[A13]

Overall, the infiltration rate is calculated from the difference between the throughfall rate and surface runoff, whereby
the throughfall rate has three components, namely interception,
snowmelt, and dripping from the interception reservoir.

A2.6 Noah-MP
Wang et al. (2016) defined the maximum infiltration rate,
Imax [L T−1], for a grid cell in Noah-MP asPx Dx éë1- exp(-kdt d t )ùû
[A14]
Px + éë1- exp(-kdt d t )ùû
where d t is the model time step [T], Dx is incorrectly termed by
Wang et al. (2016) as soil water diffusivity, when in fact it is simply
the soil water storage in length units [L], calculated as
I max =

4

Dx = åD z i ( q s -q i )

[A15]

i =1

and
kdt = kdt ref

Ks
K ref

[A16]

where Px is the precipitation rate [L T−1], Ks is the saturated
hydraulic conductivity [L T−1], q s and q i are the volumetric saturated water content and the actual water content [L3 L−3] at the
ith time step, respectively, kdt is a parameter in Eq. [A14], kdtref
is a constant parameter set equal to 3.0 (dimensionless), and Kref
is a reference K value (2 ´ 10−6 m s−1).
In Noah-MP (Niu et al., 2005, 2011), Imax is used to calculate
Hortonian excess infiltration as
p. 38 of 53

qs = Fsat Q wat + (1- Fsat )max éêë 0,( Q wat - I max )ùúû[A17]

where Fsat is calculated from TOPMODEL (see below), Q wat[LT−1] is the input of water (rainfall, snowfall, dew), and Imax isthe maximum soil infiltration capacity calculated according to
Entekhabi and Eagleson (1989), where it is defined as infiltrability.

A2.7 JULESIn the JULES model (Best et al., 2011), the maximum surface
infiltration rate is defined as Imax = b sKs, where b s is the enhancement factor (dimensionless), generally set equal to 0.5 for bare soil,
whereas larger values are used for vegetated grid cells (4 for trees, 2
for grasses and shrubs), to account for infiltration enhancing factors
such as root macropores, and Ks is the saturated hydraulic conductivity [L T−1]. This approach may lead to an underestimation of theinfiltration rate because sorptivity forces are not taken into accountand essentially all the water reaching the soil surface will infiltrate.
Hortonian runoff is calculated as the difference between throughfall
plus snowmelt and infiltration. Runoff can increase under certain
configurations to avoid supersaturation of the upper soil layer—
thus explicitly representing Dunne runoff at the point scale. If the“large-scale hydrology” scheme is used, Dunne runoff is calculated inaddition to Hortonian runoff, based on the surface saturated fraction, Fsat (see Appendix A3), according toqdunne = Fsat q ( t )[A18]where q(t) is the infiltration rate [L T−1], which is multiplied by
1 − Fsat to account for this additional runoff term. Note, thatponding is not simulated.

A2.8 CLM
The Community Land Model (CLM) version 4.5 (Oleson
et al., 2013) calculates the maximum soil infiltration capacity as
I max = (1- Fsat ) K sC ice[A19]

where Ks is the saturated hydraulic conductivity [L T−1] and Cice isan ice impedance factor (dimensionless). Hortonian excess infiltration runoff, qH [L T−1], is generated as

qH = max éë qin, soil -(1- f h20sfc ) I max ,0 ùû[A20]where f h20sfc is the fraction of the area where ponded water exists
and thus equals Fsat. Note that Hortonian excess infiltration runoff
is only generated and ponded water only occurs on the diagnosed
unsaturated fraction (1 − Fsat) of the soil column. The variable
qin,soil refers to the infiltration rate of water into the soil, defined as

qin,soil = (1- f h20 s fc ) qin,surface- qH -(1- f sno - f h20 s f c ) qevapo,soil

[A21]

and

qin,surface = (1- Fsat ) q liq,o[A22]

where qliq,o is the water input into the grid cell [L T−1], which isthe sum of liquid precipitation reaching the surface and snowmelt,

VZJ | Advancing Critical Zone Scienceqevapo,soil [L T−1] is the water loss due to evaporation, qin,surface[L T−1] is the surface water flux after surface runoff has been
removed, and fsno is the fraction of the grid cell covered by snow.

A2.11 H-TESSEL and CH-TESSELIn H-TESSEL (Balsamo et al., 2009) or in the more recent
CH-TESSEL, both based on the same hydrological principles(Boussetta et al., 2013), the maximum infiltration rate is calculated
to determine Hortonian overland flow as
I max = (Wsat -W )æìæüb+1 ö÷÷
ïçç
ö÷1/(b+1) é PT + PM ù ïWïï ÷ [A23]ïçúï÷÷+ max çç 0,Wsat íçç1- êê
÷
úý
çç
÷
ï
ï ÷÷
ç
1
W
b
+
W
(
)
è
ø
ê
sat
sat ûú ï
ï
ë
çè
ï
ï ø÷
î
þ

where P T is the throughfall precipitation [L], P M is the snowmelt [L] leading to Px as the total water reaching the surface(Px = P T + P M [L]), and W and Wsat [L] are the vertically integrated soil water contents (equivalent to q and q s , albeit withdifferent units) over the first 50 cm of the soil profile. The b
parameter reflects the grid cell heterogeneity (see above). In
H-TESSEL and CH-TESSEL, qs = Px − Imax , where qs is surface runoff (defined as R in Balsamo et al., 2009). Notice thatthe units of Imax are length, although it is defined as a rate inBalsamo et al. (2009).

66Appendix A3

Description of the Saturated Water Fraction

Note that only those models using the maximum infiltration rate in their concept will be described here and that the last
number of the label (e.g., A3.3) refers to the number given in
Table 2 to facilitate reading. Because the TOPMODEL concept
to calculate Fsat is embedded in several LSMs presented in thisappendix, we briefly present the basic equations and ways to calculate basic properties.

A3.0 TOPMODEL
In the TOPMODEL approach (Niu et al., 2005), Fsat can be
defined as

Fsat = òl³( l m + f d zÑ )pdf ( l )d l

[A24]

where fd is the decay factor [L−1], which is a measure of the declineof Ks with increasing depth, zÑ is the mean water table depth (positive for values below the ground surface), and the local water table
depth zl is given byz l = zÑ -1
( l m -l )
fd

[A25]

where l m is the mean topographic (wetness) index and l is the
local topographic index, defined asæ a öl = ln ççç T ÷÷÷è tan b T ø÷

[A26]

p. 39 of 53

where aT is the specific catchment area and tanbT is the local surface topographic slope.
To calculate the average water table depth, Chen and Kumar(2001) proposed an iterative procedure where D q is the amount ofwater stored between surface and groundwater:
Dq = ò

zÑ0

é q s ( z )-q ( z )ù d z
ë
û

[A27]

where the water contents q s(Z) and q(Z) denote the catchment orbasin average [L3 L−3]. Further, D q is calculated fromDq =N LAYERSå
i =1

é q s -q ( z ) ù D z i
i ûú
ëê

[A28]

The water content q(z) is calculated using, e.g., the Brooks–
Corey equation (Eq. [4–6]) and the assumption that the soil watercontent profile is an equilibrium with the groundwater byé hc -( zÑ - z ) ù - 1/ Búq ( z ) = q s êê
ú
h
c
ë
û

¢

[A29]
¢

æ hc -( zÑ - z ) ö÷- 1/ B			q ( z ) = q s ( z ) ççç
÷÷
hc
èç
ø÷

[A30]

where B¢ is related to the Brooks–Corey exponent of the water
retention characteristic and hc is the pressure head at air entry (cm).
This approach was used by Niu et al. (2005) to develop a runoff
scheme for global climate models, and it has also been implemented in CLM (see below)

A3.3 ORCHIDEE
The ORCHIDEE model does not use Fsat to generatesurface runoff, which is the complement of upscaled local infiltration rates (see above). In contrast, it uses the concept of theponded fraction to reduce the surface and enhance infiltration.
A fraction of surface runoff, g f, is allowed to pond in flat areas,
and it is kept to be infiltrated at the following time step withthroughfall and snowmelt to account for the effect of ponding on
infiltration (d’Orgeval et al., 2008). This fraction g is constantover time but varies spatially, based on the mean slope Sl in thegrid cells and a threshold slope Sl max (with a default value of0.5%), such that the ponding fraction decreases from 1 when Sl= 0 to 0 when Sl ³ Sl max:g f = 1- min (1,Sl Sl max )[A31]This leads to grid-scale surface runoff reduced by g f.

A3.5 ISBA-SURFEX
In the ISBA-SURFEX model, the saturation excess runoff
can also be computed using the TOPMODEL assumption insteadof using the Arno scheme, so only Dunne runoff is affected by
Fsat (Decharme and Douville, 2006). It is described according to acorrected approach of the original TOPMODEL framework proposed by Saulnier and Datin (2004) and written as
qD = Px Fsat

[A32]

where Fsat is defined as the saturation fraction of a grid cell, being
inversely proportional to the mean water storage deficit, Dt, of thegrid cell, whereby Dt can be written as
Dt = ( q s -q ) d2

[A33]

The value of Dt is bounded between 0 and do, the maximum local
water deficit, defined by
do = ( q s -q wilt ) d2

[A34]

where q s and q are the saturated and mean volumetric water content averaged across the depth d2 [L] at grid scale [L3 L−3] that canbe, optionally, the depth of the entire root zone (Decharme andDouville, 2006) or the depth of the layer in which the cumulated
root profile reached 90% (Decharme et al., 2013), and q wilt is thewater content at the wilting point [L3 L−3]. Here, Dt will be dowhen Fsat = 0 or vice versa.

A3.6 Noah, Noah-MP
We briefly describe the way infiltration is handled in Noah andthen present Noah-MP and the approaches used to calculate Fsat.
In the Noah model (Schaake et al., 1996), the spatiallyaveraged actual infiltration rate, Ir, depends on the cumulative
infiltration capacity (Ic,t) at a certain moment in time t; Ic,t isexpressed as

I c, t = Db éë1- exp(-K d t D t )ùû

[A35]

where Dt is the model time step, Kdt is a constant, and Db represents the spatially averaged soil water storage for the whole soilcolumn. The value of Db in each soil layer is computed byæ q ( i ) +q i ( i )-q wilt ( i ) ö÷÷÷
Db ( i ) = Dbmax ( i )ççç1- l
q s ( i )-q wilt ( i ) ø÷
èç

[A36]

where Dbmax(i), q l(i), q i(i), q wilt(i), and q s(i) are the maximum soil
water storage, the content of liquid water and ice, the wilting point,and soil porosity or saturated water content, respectively, in the ithsoil layer; Dbmax(i) is defined by

A3.4 Catchment Land Surface Model

Dbmax ( i ) =-z ( i ) éë q s ( i ) -q wilt ( i )ùû

In the Catchment Land Surface Model (CLSM), the saturated land fraction under equilibrium conditions is effectively
computed as the fraction of the area for which the water table
depth lies “above” the ground surface, based on the TOPMODEL
framework. The strategy for calculating the saturated land fraction
was described above.

where z(i) is the thickness of the ith soil layer.
Then, the infiltration rate for the whole soil column is given by

VZJ | Advancing Critical Zone Science

Ir =

[A37]

PR I c, t ( PR + I c, t )

[A38]

Dt

where P R is the precipitation input into the whole column.
p. 40 of 53

Soil freezing has a significant effect on the permeability ofsoils because ice impedes the infiltration rate. Noah computes theimpermeable area factor FC to consider the influence of freezing
on soil infiltration following Koren et al. (1999).
Noah-MP is an improved version of Noah, whereby Noah-MP(Niu et al., 2011) offers four options for computing infiltration andsurface–subsurface runoff. Option 1 is the TOPMODEL‐based
runoff scheme with the simple groundwater scheme (Niu and Yang,2007) implemented by Cai et al. (2014), where the fraction of thegridbox that is saturated, Fsat, is given by

Fsat = (1- Ffrz ) Fmax exp éë - 0.5 f d ( z wt - z bot )ùû + Ffrz

[A39]

where zbot is the soil profile thickness (default = 2.0 m), zwt (L)is the water table depth, Fmax is set to a global mean value of0.38, and Ffrz is the fraction of impermeable area as a functionof soil ice content of the surface soil layer. For urban areas, Ffrz
is set to 0.95. The runoff decay factor, fd , in Eq. [A39] equals
6 m−1.
Option 2 is a simple TOPMODEL‐based runoff schemewith an equilibrium water table (Niu et al., 2005). As in Option
1, this scheme parameterizes both surface and subsurface runoff
as functions of the water table depth but with a sealed bottomof the soil column (zero‐f lux lower boundary condition) inaccordance with one of the TOPMODEL assumptions, i.e., theexponential decay of saturated hydraulic conductivity; Fsat is
calculated as
Fsat = (1- Ffrz ) Fmax exp(- 0.5 f d z wt ) + Ffrz

[A40]

where fd is set to 2.0 m−1.
Option 3 is an infiltration‐excess‐based surface runoff schemewith a gravitational free‐drainage subsurface runoff scheme as usedin the original Noah (Schaake et al., 1996). Surface runoff (R) is
computed as the excess of precipitation (Pd) not infiltrated into
the soil (R = Pd − Imax). The maximum infiltration rate, Imax, is
computed as

I max = Pd

Dx éë1- exp(-kdt D t )ùû
Pd + Dx éë1- exp(-kdt D t )ùû

[A41]

N

[A42]

i =1

and
kdt = kdt ref

Ks
K ref

[A43]

where Dt is the time step [T], Ks [L T−1] is the saturated hydraulicconductivity, which depends on soil texture and is prescribed ina lookup table, and N is the number of layers. The parameters
kdtref = 3.0 and Kref = 2 ´ 10−6 m s−1 were determined in theframework of the PILPS 2(c) experiments for the Red–Arkansas
River basins in the Southern Great Plains region of the UnitedStates (Wood et al., 1998).

VZJ | Advancing Critical Zone Scienceæ q tot ö÷4÷÷ + Ffrz
Fsat = (1- Ffrz )çççèç q ÷ø

[A44]

s

with
q tot =

1NåD z i q iz bot i=1

[A45]

A3.7 JULES
There are two options in JULES to account for the spatial
heterogeneity of soil water content and thus determine Fsat. The
first option is based on a modified form of TOPMODEL (Beven
and Kirkby, 1979), described by Gedney and Cox (2003), where
the assumption of an exponential decay of Ks with depth (leading
to Eq. [A24–A25]) is relaxed. Grid box mean water table depth, zÑ,is calculated using the approach described for A3.0 TOPMODEL(Eq. [A27–A29]). This is used to estimate a critical topographic
index l crit using the relation
é T (0) ùú +ll crit = ln êêmúëê T ( zÑ ) ûú

[A46]

where transmissivity T(zÑ) is given by
T ( zÑ ) = ò¥zK s ( z )d z

[A47]

The value of Fsat is then calculated as the fraction of the grid box
for which l > l crit, assuming l follows a g distribution with mean(l m) and standard deviation read in from observational datasets.
To speed up the computation, the integral over the probability distribution of l is calculated for a range of mean water table depths
during model initialization and approximated by an exponentialfunction as
Fsat = as exp (-c s l crit )with
Dx = åD z i ( q s -q i )Finally, Option 4 is the BATS runoff scheme, which parameterizes surface runoff as a fourth power function of the top 2-m
soil wetness (expressed as degree of saturation) and subsurface
runoff as gravitational free drainage (Yang and Dickinson, 1996),
described by[A48]

where as and cs are fitted to approximate the full integral (i.e., Eq.
[A24]). Recently, this scheme has been modified to account for the
impact of frozen water in the soil by replacing Ks in Eq. [A47] by(1 − q f )2b+3Ks (where b is the Brooks–Corey parameter and q f is
the frozen water content) and additionally correcting the grid box
mean water table depth (zÑ) to account for the fact that the profile
of soil moisture above the water table will not follow the sameequilibrium profile (assumed in Eq. [A29]) when ice is present.
The second option for calculating Fsat is the use of the probability
distributed model of Moore (1985) (see Eq. [38]).

A3.8 CLM 4.5
In CLM 4.5 (Lawrence et al., 2011; Oleson et al., 2013), thefraction of the saturated area, Fsat is calculated asp. 41 of 53

Fsat = Fmax exp(- 0.5 f d zÑ )[A49]

where Fmax is the maximum value of Fsat and fd is the decay factor[L−1]; Fmax is defined as the fraction of subgrid cells from a highresolution digital elevation map in a grid cell whose topographic
index (the ratio of the upstream area to the slope; Niu et al., 2005)
is larger than or equal to the grid cell mean topographic index. It
is the value of the discrete cumulative distribution function of the
topographic index when the grid cell mean water table depth is zero.

A3.9 CABLE
The original CABLE model (version 2, Kowalczyk et al.,2013) generates surface runoff from excess infiltration only when
the first three soil layers are at least 95% saturated. There is noother surface runoff generation process. However, Decker (2015)
implemented subgrid-scale soil water content variability, explicit
runoff generation, and groundwater in CABLE. The fraction of
the saturated area, Fsat, follows Entekhabi and Eagleson (1989)
and is defined as

Fsat = ò¥kp

f ( s )d s

[A50]

where f (s) is the Gamma distribution of Eq. [42], and kp is definedby Eq. [A53]. One can solve this integral assuming a constant a c
in Eq. [42], which then links s = (q − q r)/(q s − q r), the mean of thevertically averaged relative saturation over the grid cell, with l s dueto the properties of the Gamma distribution (Eq. [42]):
ls =

1
2s

[A51]

For a c = 1/2,
æ k ÷öç p÷[A52]
Fsat = 1- erf çç÷ççè l s ÷÷ø
Finally, k p is parameterized by an empirical formulation as
k p = C v sslope[A53]with Cv as a fitting parameter and sslope as the standard deviationof the subgrid slope.
Haverd et al. (2016) and Cuntz and Haverd (2018) recently
implemented a new soil (Haverd and Cuntz, 2010) and snow
model with physically accurate freeze–thaw processes within
CABLE, which is currently combined with the developments of
Decker (2015).

A3.10 SSiB
In the SSiB model, the normalized runoff Iroff is spatially distributed as a function of the fractional area of grid area x (0 < x < 1)
(Sato et al., 1989):
I roff ( x ) = aSSiB exp(-bSSiB x ) + c SSiB

[A54]

where aSSiB, bSSiB, and cSSiB are constants. This distribution has
also been applied to the convective precipitation. The constants (a,VZJ | Advancing Critical Zone Scienceb, and c) were obtained by comparison with the observational data
and are normalized so that1

ò0 I roff ( x )d x = 1[A55]

Based on the above two spatial distributions, the saturation
fraction, Fsat, could be obtained (Sato et al., 1989) asæ K D t ö÷ c
1
ç
s
SSiB
Fsat =log çç
[A56]
÷÷bSSiBçè Pdrop aSSiB ÷÷ø aSSiB
where Ks is the saturated hydraulic conductivity [L T−1], Dt is
the time interval [T], and Pdrop is the precipitation reaching thesurface after interception by the canopy [L T−1]. The spatial distribution of convective precipitation has also been applied to obtainthe value of Pdrop (Sellers et al., 1996).

A3.12 JSBACH
JSBACH uses the original Arno scheme (Dümenil and Todini,1992) to determine surface runoff and infiltration. Accordingly,
the saturated fraction is estimated for the grid box fraction forwhich the soil water capacity of the root zone is less than or equal
to the grid box mean root zone soil moisture. The shape parameterb is determined using Eq. [45].

66Appendix A4

Use of Pedotransfer Functions for the Estimationof Saturated Hydraulic Conductivity

For soils with low organic C content, Noah-MP and CLM 4.5
use the PTFs developed at point scale by Clapp and Hornberger(1978) and Cosby et al. (1984) to estimate Ks. In ISBA-SURFEX,
Ks is related to the soil textural properties (clay and sand) using theNoilhan and Lacarrère (1995) continuous relationships derived
from theses PTFs. These PTFs are basically class PTFs (Van Looy
et al., 2017) that give average values of Brooks–Corey parameters
and the measured Ks for each textural class of the USDA classification and are also developed at point scale.
Soil classes in CLSM are parameterized by Campbell (1974)
equations, and the corresponding hydraulic parameters including
Ks are based on lookup tables for 12 different soil textural classes orusing the PTFs of Wösten et al. (2001), which are both developed
at the point scale. Campbell’s method used the Brooks–Coreyparameterization of the soil water retention curve and a single
point measurement of K at a given water content to calculate thecomplete hydraulic conductivity function.
None of the models directly considers the effect of soil structure on saturated hydraulic conductivity. Only OLAM-SOIL
considers implicitly the impact of structural properties on Ks bylinearly interpolating between the measured Ks value and the valueof the hydraulic conductivity obtained at a pressure head of about−6 cm as proposed by Weynants et al. (2009). Based on the work of
Jarvis (2007), this value was considered to delineate the saturation
range that is controlled by structural properties.
p. 42 of 53

In JSBACH, Ks values are assigned to 11 textural classes basedon data presented by Beringer et al. (2001). However, the originof the tabulated Ks values, which appear to be mean measured
values for a specific textural class, is not clear. Most likely, they
were derived from the dataset of Clapp and Hornberger (1978).
JULES uses sand and clay content to estimate the Ks for themain soil column. Below the soil column is an additional “deep water
store,” within which Ks decreases exponentially with depth, with adampening factor set equal to 3 as proposed by Clark and Gedney(2008). The Ks values can also be defined for each soil layer and canaccount for the presence of soil organic matter (Chadburn et al.,2015a). Rahman and Rosolem (2017) incorporated the effect of preferential flow into JULES (but note that this has not yet been adopted
in the current “official” UKMO version of JULES) to allow simulation of highly fractured unsaturated chalk soils. Their bulk hydraulicconductivity (BC) model introduces only two additional parameters(namely the macroporosity factor and the soil wetness thresholdparameter for fracture flow activation) and uses the Ks from the chalkmatrix. The BC model was implemented into JULES and applied
to a study area encompassing the Kennet catchment in the southern
UK, and the model performance at the catchment scale was evaluatedagainst independent datasets (e.g., runoff and latent heat flux). The
results demonstrated that the inclusion of the BC model in JULES
improved the simulations of land surface water and energy fluxes
across the chalk-dominated Kennet catchment. This simple approachto account for soil structure has potential for large-scale land surface modeling applications. ORCHIDEE uses van Genuchten soil
parameters (Ks, n, a, q r, and q s) for each USDA class. In addition, adecay of Ks with depth is imposed (as also in JULES), as describedabove. With respect to the horizontal variability, ORCHIDEE uses
an exponential PDF to describe horizontal heterogeneity.

Alaoui, A., J. Lipiec, and H.H. Gerke. 2011. A review of the changes in thesoil pore system due to soil deformation: A hydrodynamic perspective.
Soil Tillage Res. 115-116:1–15. doi:10.1016/j.still.2011.06.002

Acknowledgments

Assouline, S., K. Narkis, R. Gherabli, and G. Sposito. 2016. Combined
effect of sodicity and organic matter on soil properties underlong-term irrigation with treated wastewater. Vadose Zone J. 15(4).
doi:10.2136/vzj2015.12.0158

This review emerged from discussions among the GEWEX-SoilWat community
held at the workshop in Leipzig (Germany) in 2016 and is a product of the International Soil Modeling Consortium (ISMC) effort. We would like to thank Randal
Koster (NASA-GMAO) for discussion of CLSM parts.

References
Abramopoulos, F., C. Rosenzweig and B. Choudhury. 1988. Improved
ground hydrology calculations for global climate models (GCMs):
Soil water movement and evapotranspiration. J. Climate 1:921–941.
doi:10.1175/1520-0442(1988)001<0921:ighcfg>2.0.co;2
Adem, H.H., and S.K. Vanapalli. 2015. Review of methods for predicting
in situ volume change movement of expansive soil over time. J. Rock
Mech. Geotech. Eng. 7:73–86. doi:10.1016/j.jrmge.2014.11.002
Ahuja, L.R. 1974. Applicability of the Green–Ampt approach towater infiltration through surface crust. Soil Sci. 118:283–288.
doi:10.1097/00010694-197411000-00001
Ahuja, L.R. 1983. Modeling infiltration into crusted soils by theGreen–Ampt approach. Soil Sci. Soc. Am. J. 47:412–418.
doi:10.2136/sssaj1983.03615995004700030004x
Ahuja, L.R., D.K. Cassel, R.R. Bruce, and B.B. Barnes. 1989. Evaluation of spatial distribution of hydraulic conductivity using effective porosity data.
Soil Sci. 148:404–411. doi:10.1097/00010694-198912000-00002
Ahuja, L.R., F. Fiedler, G.H. Dunn, J.G. Benjamin, and A. Garrison.
1998. Changes in soil water retention curves due to tillage
and natural reconsolidation. Soil Sci. Soc. Am. J. 62:1228–1233.
doi:10.2136/sssaj1998.03615995006200050011x

VZJ | Advancing Critical Zone ScienceAnwar, S.A., A.S. Zakey, S.M. Robaa, and M.M. Abdel Wahab. 2018. The
influence of two land-surface hydrology schemes on the regional
climate of Africa using the RegCM4 model. Theor. Appl. Climatol.
doi:10.1007/s00704-018-2556-8
Archer, N.A.L., M. Bonell, N. Coles, A.M. MacDonald, C.A. Auton, and R.
Stevenson. 2013. Soil characteristics and landcover relationships on soil
hydraulic conductivity at a hillslope scale: A view towards local flood
management. J. Hydrol. 497:208–222. doi:10.1016/j.jhydrol.2013.05.043
Archer, N.A.L., J.N. Quinton, and T.M. Hess. 2012. Patch vegetation andwater redistribution above and below ground in south-east Spain.
Ecohydrology 5:108–120. doi:10.1002/eco.210
Arnold, J.G., K.N. Potter, K.W. King, and P.M. Allen. 2005. Estimation of soil
cracking and the effect on surface runoff in a Texas Blackland Prairie
watershed. Hydrol. Processes 19:589–603. doi:10.1002/hyp.5609
Assouline, S. 2004. Rainfall-induced soil surface sealing: A critical review
of observations, conceptual models, and solutions. Vadose Zone J.
3:570–591. doi:10.2136/vzj2004.0570
Assouline, S. 2006a. Modeling the relationship between soil bulk density
and the hydraulic conductivity function. Vadose Zone J. 5:697–705.
doi:10.2136/vzj2005.0084
Assouline, S. 2006b. Modeling the relationship between soil bulk
density and the water retention curve. Vadose Zone J. 5:554–563.
doi:10.2136/vzj2005.0083
Assouline, S. 2013. Infiltration into soils: Conceptual approaches and solutions. Water Resour. Res. 49:1755–1772. doi:10.1002/wrcr.20155
Assouline, S. and M. Ben-Hur. 2006. Effects of rainfall intensity and slope
gradient on the dynamics of interrill erosion during soil surface sealing. doi:10.1016/j.catena.2006.02.005
Assouline, S. and Y. Mualem. 2002. Infiltration during soil sealing: The
effect of areal heterogeneity of soil hydraulic properties. Water Resour
Res 38(12):1286. doi:10.1029/2001WR001168
Assouline, S., and Y. Mualem. 2006. Runoff from heterogeneous small barecatchments during soil surface sealing. Water Resour. Res. 42:W12405.
doi:10.1029/2005WR004592
Assouline, S., and K. Narkis. 2011. Effects of long-term irrigation withtreated wastewater on the hydraulic properties of a clayey soil. Water
Resour. Res. 47:W08530. doi:10.1029/2011WR010498

Assouline, S., and D. Or. 2008. Air entry–based characteristic length for estimation of permeability of variably compacted earth materials. Water
Resour. Res. 44:W11403. doi:10.1029/2008WR006937
Assouline, S., and D. Or. 2013. Conceptual and parametric representation of soil hydraulic properties: A review. Vadose Zone J. 12(4).
doi:10.2136/vzj2013.07.0121
Assouline, S., J.S. Selker, and J.-Y. Parlange. 2007. A simple accurate
method to predict time of ponding under variable intensity rainfall.
Water Resour. Res. 43:W03426. doi:10.1029/2006WR005138
Assouline, S., S.E. Thompson, L. Chen, T. Svoray, S. Sela, and G.G. Katul.
2015. The dual role of soil crusts in desertification. J. Geophys. Res.
Biogeosci. 120:2108–2119. doi:10.1002/2015JG003185
Bachmann, J., and R.R. van der Ploeg. 2002. A review on recent developments in soil water retention theory: Interfacial tension and temperature effects. J. Plant Nutr. Soil Sci. 165:468–478. doi:10.1002/15222624(200208)165:4<468::AID-JPLN468>3.0.CO;2-G
Balsamo, G., A. Beljaars, K. Scipal, P. Viterbo, B. van den Hurk, M. Hirschi,and A.K. Betts. 2009. A revised hydrology for the ECMWF model:
Verification from field site to terrestrial water storage and impact
in the integrated forecast system. J. Hydrometeorol. 10:623–643.
doi:10.1175/2008JHM1068.1
Barry, D.A., J.-Y. Parlange, R. Haverkamp, and P.J. Ross. 1995. Infiltration
under ponded conditions: 4. Explicit predictive infiltration formula.
Soil Sci. 160:8–17. doi:10.1097/00010694-199507000-00002

p. 43 of 53

Barry, D.A., J.-Y. Parlange, L. Li, D.-S. Jeng, and M. Crapper. 2005.
Green Ampt approximations. Adv. Water Resour. 28:1003–1009.
doi:10.1016/j.advwatres.2005.03.010
Barry, D.A., J.-Y. Parlange, M.-C. Liu, G.C. Sander, M.B. Parlange, D.A. Lockington, et al. 2007. Infiltration and ponding. In: J.A. Filar, editor, Encyclopedia of life support systems: Mathematical sciences. Mathematicalmodeling. Vol. II. EOLSS Publ., Oxford, UK.
Batjes, N.H. 2006. ISRIC-WISE derived soil properties on a 5 by 5 arcminutes global grid. Rep. 2006/02. ISRIC World Soil Inf., Wageningen,the Netherlands.
Baumhardt, R.L., M.J.M. Römkens, F.D. Whisler, and J.-Y. Parlange. 1990.
Modeling infiltration into a sealing soil. Water Resour. Res. 26:2497–
2505. doi:10.1029/WR026i010p02497
Bear, J. 1972. Dynamics of fluids in porous media. Elsevier, New York.
Bechtold, M., G.J.M. De Lannoy, R.D. Koster, R.H. Reichle, S. Mahanama, W.
Bleuten, et al. PEAT-CLSM: A specific treatment of peatland hydrology
in the NASA Catchment Land Surface model. J. Adv. Model. Earth Syst.
doi:10.1029/2018MS001574
Beck, H.E., A. van Dijk, A. de Roo, E. Dutra, G. Fink, R. Orth, and J.
Schellekens. 2017. Global evaluation of runoff from 10 state-ofthe-art hydrological models. Hydrol. Earth Syst. Sci. 21:2881–2903.
doi:10.5194/hess-21-2881-2017
Bengough, A. 2012. Water dynamics of the root zone: Rhizosphere
biophysics and its control on soil hydrology. Vadose Zone J. 11(2).
doi:10.2136/vzj2011.0111
Ben Neriah, A., S. Assouline, U. Shavit, and N. Weisbrod. 2014. Impact of
ambient conditions on evaporation from porous media. Water Resour.
Res. 50:6696–6712. doi:10.1002/2014WR015523

Bouwer, H. 1964. Unsaturated flow in ground-water hydraulics. J. Hydraul.
Eng. 90(HY5):121–144.
Brakensiek, D.L., and W.J. Rawls. 1983. Agricultural management effects
on water processes: II. Green and Ampt parameters for crusting soils.
Trans. ASAE 26:1753–1757. doi:10.13031/2013.33838
Bresler, E. 1973. Simultaneous transport of solutes and water undertransient unsaturated flow conditions. Water Resour. Res. 9:975–986.
doi:10.1029/WR009i004p00975
Bresler, E., and R. Hanks. 1969. Numerical method for estimating simultaneous flow of water and salt in unsaturated soils. Soil Sci. Soc. Am. J.
33:827–832. doi:10.2136/sssaj1969.03615995003300060013x
Bresler, E., B.L. McNeal, and D.L. Carter. 1982. Saline and sodic soils: Principles–dynamics–modeling. Springer, Berlin.
doi:10.1007/978-3-642-68324-4
Bronick, C.J., and R. Lal. 2005. Soil structure and management: A review.
Geoderma 124:3–22. doi:10.1016/j.geoderma.2004.03.005
Bronswijk, J.J.B. 1991. Relation between vertical soil movements andwater-content changes in cracking clays. Soil Sci. Soc. Am. J. 55:1220–
1226. doi:10.2136/sssaj1991.03615995005500050004x
Brooks, R., and A. Corey. 1964. Hydraulic properties of porous media.
Hydrol. Pap. 3. Colorado State Univ., Fort Collins.
Brown, G.G., I. Barois, and P. Lavelle. 2000. Regulation of soil organic matter dynamics and microbial activity in the drilosphere and the role ofinteractions with other edaphic functional domains. Eur. J. Soil Biol.
36:177–198. doi:10.1016/S1164-5563(00)01062-1
Brutsaert, W. 2005. Hydrology: An introduction. Cambridge Univ. Press,
Cambridge, UK. doi:10.1017/CBO9780511808470

Beringer, J., A.H. Lynch, F.S. Chapin, M. Mack and G.B. Bonan.
2001. The representation of arctic soils in the land surface
model: The importance of mosses. J. Climate 14:3324–3335.
doi:10.1175/1520-0442(2001)014<3324:Troasi>2.0.Co;2

Buckingham, E. 1907. Studies on the movement of soil moisture. Bureauof Soils Bull. 38. US Gov. Print. Office, Washington, DC.

Best, M.J., P. Cox, and D. Warrilow. 2005. Determining the optimal soil temperature scheme for atmospheric modeling applications. BoundaryLayer Meteorol. 114:111–142. doi:10.1007/s10546-004-5075-3

Cai, X.T., Z.L. Yang, Y.L. Xia, M.Y. Huang, H.L. Wei, L.R. Leung, and M.B. Ek.
2014. Assessment of simulated water balance from Noah, Noah-MP,
CLM, and VIC over CONUS using the NLDAS test bed. J. Geophys. Res.
Atmos. 119:13751–13770. doi:10.1002/2014JD022113

Best, M.J., M. Pryor, D.B. Clark, G.G. Rooney, R.L.H. Essery, C.B. Menard, etal. 2011. The Joint UK Land Environment Simulator (JULES), model
description: 1. Energy and water fluxes. Geosci. Model Dev. 4:677–699.
doi:10.5194/gmd-4-677-2011
Betson, R.P., and J.B. Marius. 1969. Source areas of storm runoff. Water
Resour. Res. 5:574–582. doi:10.1029/WR005i003p00574
Beven, K. 1984. Infiltration into a class of vertically non-uniform soils.
Hydrol. Sci. J. 29:425–434. doi:10.1080/02626668409490960
Beven, K., and P. Germann. 1982. Macropores and water flow in soils. Water Resour. Res. 18:1311–1325. doi:10.1029/WR018i005p01311
Beven, K.J. and M.J. Kirkby. 1979. A physically based, variable contributing
area model of basin hydrology. Hydrol. Sci. Bull. 24:43–69.
Bierkens, M.F.P. 2015. Global hydrology 2015: State, trends, and directions.
Water Resour. Res. 51:4923–4947. doi:10.1002/2015WR017173
Bisal, F., and K.F. Nielsen. 1967. Effect of frost action on the size of soil aggregates. Soil Sci. 104:268–272. doi:10.1097/00010694-196710000-00006
Bonan, G.B. 1998. The land surface climatology of the NCAR Land Surface
Model coupled to the NCAR Community Climate Model. J. Climate
11:1307–1326. doi:10.1175/1520-0442(1998)011<1307:Tlscot>2.0.Co;2
Bonan, G.B., K.W. Oleson, M. Vertenstein, S. Levis, X.B. Zeng, Y.J. Dai, et al.
2002. The land surface climatology of the Community Land Model
coupled to the NCAR Community Climate Model. J. Climate 15:3123–
3149. doi:10.1175/1520-0442(2002)015<3123:Tlscot>2.0.Co;2
Boone, A., V. Masson, T. Meyers, and J. Noilhan. 2000. The influenceof the inclusion of soil freezing on simulations by a soil–vegetation–atmosphere transfer scheme. J. Appl. Meteorol. 39:1544–1569.
doi:10.1175/1520-0450(2000)039<1544:TIOTIO>2.0.CO;2

Busscher, W. 2011. Hardpan soils: Management. Springer, Dordrecht,
the Netherlands.

Campbell, G.S. 1974. Simple method for determining unsaturated
conductivity from moisture retention data. Soil Sci. 117:311–314.
doi:10.1097/00010694-197406000-00001
Carsel, R.F., and R.S. Parrish. 1988. Developing joint probability distributions of soil water retention characteristics. Water Resour. Res.
24:755–769. doi:10.1029/WR024i005p00755
Celia, M.A., E.T. Bouloutas, and R.L. Zarba. 1990. A general mass‐conservative numerical solution for the unsaturated flow equation. Water
Resour. Res. 26:1483–1496. doi:10.1029/WR026i007p01483
Chadburn, S., E. Burke, R. Essery, J. Boike, M. Langer, M. Heikenfeld, et al.
2015a. An improved representation of physical permafrost dynamics
in the JULES land-surface model. Geosci. Model Dev. 8:1493–1508.
doi:10.5194/gmd-8-1493-2015
Chadburn, S.E., E.J. Burke, R.L.H. Essery, J. Boike, M. Langer, M. Heikenfeld,
et al. 2015b. Impact of model developments on present and future
simulations of permafrost in a global land-surface model. Cryosphere9:1505–1521. doi:10.5194/tc-9-1505-2015
Chamberlain, E.J., and A.J. Gow. 1979. Effect of freezing and thawing
on the permeability and structure of soils. Eng. Geol. 13:73–92.
doi:10.1016/0013-7952(79)90022-X
Chen, J. and P. Kumar. 2001. Topographic influence on the seasonal and interannual variation of water and energy balance of basins in North America. J. Climate 14:1989–2014.
doi:10.1175/1520-0442(2001)014<1989:Tiotsa>2.0.Co;2Childs, E.C. 1969. An introduction to the physical basis of soil water phenomena. John Wiley & Sons, New York.

Boone, A., and P.J. Wetzel. 1999. A simple method for modeling sub-gridsoil texture variability for use in an atmospheric climate model. J. Meteorol. Soc. Jpn. 77:317–333. doi:10.2151/jmsj1965.77.1B_317

Childs, E.C., and M. Bybordi. 1969. The vertical movement of water instratified porous material: 1. Infiltration. Water Resour. Res. 5:446–459.
doi:10.1029/WR005i002p00446

Boussetta, S., G. Balsamo, A. Beljaars, A.-A. Panareda, J.-C. Calvet, C. Jacobs,
et al. 2013. Natural land carbon dioxide exchanges in the ECMWF
integrated forecasting system: Implementation and offline validation.
J. Geophys. Res. Atmos. 118:5923–5946. doi:10.1002/jgrd.50488

Choi, H.I., P. Kumar, and X.Z. Liang. 2007. Three-dimensional volumeaveraged soil moisture transport model with a scalable parameterization of subgrid topographic variability. Water Resour. Res. 43:W04414.
doi:10.1029/2006WR005134

VZJ | Advancing Critical Zone Science

p. 44 of 53

Chow, V.T., D.R. Maidment, and L.W. Mays. 1988. Applied hydrology.
McGraw-Hill, Singapore.
Chu, S.T., C.A. Onstad, and W.J. Rawls. 1986. Field evaluation of layered
Green–Ampt model for transient crust conditions. Trans. ASAE
29:1268–1272. doi:10.13031/2013.30307
Chu, X., and M.A. Marino. 2005. Determination of ponding condition andinfiltration in layered soils under unsteady rainfall. J. Hydrol. 313:195–207. doi:10.1016/j.jhydrol.2005.03.002
Clapp, R.B., and G.M. Hornberger. 1978. Empirical equations forsome soil hydraulic properties. Water Resour. Res. 14:601–604.
doi:10.1029/WR014i004p00601
Clark, D.B., and N. Gedney. 2008. Representing the effects of subgrid variability of soil moisture on runoff generation in a land surface model. J.
Geophys. Res. Atmos. 113:D10111. doi:10.1029/2007JD008940
Clark, D.B., L.M. Mercado, S. Sitch, C.D. Jones, N. Gedney, M.J. Best, etal. 2011. The Joint UK Land Environment Simulator (JULES), model
description: 2. Carbon fluxes and vegetation dynamics. Geosci. Model
Dev. 4:701–722. doi:10.5194/gmd-4-701-2011
Clark, M.P., Y. Fan, D.M. Lawrence, J.C. Adam, D. Bolster, D.J. Gochis,
et al. 2015. Improving the representation of hydrologic processes in earth system models. Water Resour. Res. 51:5929–5956.
doi:10.1002/2015WR017096
Clothier, B.E. 2001. Infiltration. In: K.A. Smith and C.E. Mullins, editors, Soiland environmental analysis: Physical methods. 2nd ed. Marcel Dekker,
New York. p. 239–280.
Colman, E.A., and G.B. Bodman. 1945. Moisture and energy conditions during downward entry of water into moist and layered soils. Soil Sci. Soc.
Am. Proc. 9:3–11. doi:10.2136/sssaj1945.036159950009000C0001x
Corradini, C., A. Flammini, R. Morbidelli, and R.S. Govindaraju. 2011. A
conceptual model for infiltration in two-layered soils with a more
permeable upper layer: From local to field scale. J. Hydrol. 410:62–72.
doi:10.1016/j.jhydrol.2011.09.005
Corradini, C., R.S. Govindaraju, and R. Morbidelli. 2002. Simplified modelling of areal average infiltration at the hillslope scale. Hydrol. Processes
16:1757–1770. doi:10.1002/hyp.394
Corradini, C., F. Melone, and R.E. Smith. 1994. Modeling infiltration during complex rainfall sequences. Water Resour. Res. 30:2777–2784.
doi:10.1029/94WR00951
Corradini, C., F. Melone, and R.E. Smith. 1997. A unified model for infiltration and redistribution during complex rainfall patterns. J. Hydrol.
192:104–124. doi:10.1016/S0022-1694(96)03110-1
Cosby, B.J., G.M. Hornberger, R.B. Clapp, and T.R. Ginn. 1984. A statistical exploration of the relationships of soil moisture characteristics
to the physical properties of soils. Water Resour. Res. 20:682–690.
doi:10.1029/WR020i006p00682
Cosh, M.H., J.R. Stedinger, and W. Brutsaert. 2004. Variability of surface
soil moisture at the watershed scale. Water Resour. Res. 40:W12513.
doi:10.1029/2004WR003487
Craig, J.R., G. Liu, and E.D. Soulis. 2010. Runoff–infiltration partitioning
using an upscaled Green–Ampt solution. Hydrol. Processes 24:2328–
2334. doi:10.1002/hyp.7601
Crockford, H., S. Topalidis, and D.P. Richardson. 1991. Water repellency in adry sclerophyll eucalypt forest: Measurements and processes. Hydrol.
Processes 5:405–420. doi:10.1002/hyp.3360050408
Cueto-Felgueroso, L., and R. Juanes. 2009. A phase field model of unsaturated flow. Water Resour. Res. 45:W10409. doi:10.1029/2009WR007945
Cuntz, M., and V. Haverd. 2018. Physically accurate soil freeze–thaw
processes in a global land surface scheme. J. Adv. Model. Earth Syst.
10:54–77. doi:10.1002/2017MS001100
Cuntz, M., J. Mai, L. Samaniego, M. Clark, V. Wulfmeyer, O. Branch, et al.
2016. The impact of standard and hard-coded parameters on thehydrologic fluxes in the Noah-MP land surface model. J. Geophys. Res.
Atmos. 121:10676–10700. doi:10.1002/2016JD025097
Dagan, G., and E. Bresler. 1983. Unsaturated flow in spatially variable
fields: 1. Derivation of models of infiltration and redistribution. Water
Resour. Res. 19:413–420. doi:10.1029/WR019i002p00413
Dahan, O., B. Tatarsky, Y. Enzel, C. Kulls, M. Seely, and G. Benito.
2008. Dynamics of flood water infiltration and ground water recharge in hyperarid desert. Groundwater 46:450–461.
doi:10.1111/j.1745-6584.2007.00414.x

VZJ | Advancing Critical Zone ScienceDarcy, H. 1856. Les fontaines publiques de la ville de Dijon: Exposition etapplication des principes a suivre et des formulesa employer dans lesquestions de distribution d’eau. V. Dalmont, Paris.
Das, N.N., B.P. Mohanty, and E.G. Njoku. 2008. A Markov chain Monte Carloalgorithm for upscaled soil–vegetation–atmosphere-transfer modeling to evaluate satellite-based soil moisture measurements. Water
Resour. Res. 44:W05416. doi:10.1029/2007WR006472
Davidson, J., D. Nielsen, and J.W. Biggar. 1963. The measurement anddescription of water flow through Columbia silt loam and Hesperia
sandy loam. Hilgardia 34:601–617. doi:10.3733/hilg.v34n15p601
Davidson, S.E., and J.B. Page. 1956. Factors influencing swelling and shrinking in soils. Soil Sci. Soc. Am. J. 20:320–324.
doi:10.2136/sssaj1956.03615995002000030007x
de Boer, F. 2016. HiHydroSoil: A high resolution soil map of hydraulicproperties, Version 1.2. FutureWater Rep. 134. FutureWater, Wageningen, the Netherlands.
Decharme, B., A. Boone, C. Delire, and J. Noilhan. 2011. Local evaluationof the Interaction between Soil Biosphere Atmosphere soil multilayer
diffusion scheme using four pedotransfer functions. J. Geophys. Res.
116:D20126. doi:10.1029/2011jd016002
Decharme, B., E. Brun, A. Boone, C. Delire, P. Le Moigne, and S. Morin. 2016.
Impacts of snow and organic soils parameterization on northernEurasian soil temperature profiles simulated by the ISBA land surface
model. Cryosphere 10:853–877. doi:10.5194/tc-10-853-2016
Decharme, B., and H. Douville. 2006. Introduction of a sub-grid hydrology in the ISBA land surface model. Clim. Dyn. 26:65–78.
doi:10.1007/s00382-005-0059-7
Decharme, B., E. Martin, and S. Faroux. 2013. Reconciling soil thermal andhydrological lower boundary conditions in land surface models. J.
Geophys. Res. Atmos. 118:7819–7834. doi:10.1002/jgrd.50631
Decker, M. 2015. Development and evaluation of a new soil moisture and runoff parameterization for the CABLE LSM including
subgrid-scale processes. J. Adv. Model. Earth Syst. 7:1788–1809.
doi:10.1002/2015MS000507
De Lannoy, G.J.M., R.D. Koster, R.H. Reichle, S.P.P. Mahanama, and Q. Liu.
2014. An updated treatment of soil texture and associated hydraulicproperties in a global land modeling system. J. Adv. Model. Earth Syst.
6:957–979. doi:10.1002/2014MS000330
Delleur, J.W. 2006. The handbook of groundwater engineering. CRC Press,
Boca Raton, FL.
Demaria, E.M., B. Nijssen, and T. Wagener. 2007. Monte Carlo sensitivity analysis of land surface parameters using the Variable
Infiltration Capacity model. J. Geophys. Res. Atmos. 112:D11113.
doi:10.1029/2006JD007534
de Rosnay, P., J. Polcher, M. Bruen and K. Laval. 2002. Impact of a physically based soil water flow and soil–plant interaction representation
for modeling large-scale land surface processes. J. Geophys. Res.
107(D11). doi:10.1029/2001JD000634.
DiCarlo, D.A. 2004. Experimental measurements of saturation overshoot on infiltration. Water Resour. Res. 40: W04215.
doi:10.1029/2003WR002670
DiCarlo, D.A. 2013. Stability of gravity-driven multiphase flow in porous
media: 40 Years of advancements. Water Resour. Res. 49:4531–4544.
doi:10.1002/wrcr.20359
Dinka, T.M., and R.J. Lascono. 2012. Review paper: Challenges and limitations in studying the shrink–swell and crack dynamics of Vertisol soils.
Open J. Soil Sci. 2:82–90. doi:10.4236/ojss.2012.22012
Dirmeyer, P.A., and F.R.J. Zeng. 1999. Precipitation infiltration in thesimplified SiB land surface scheme. J. Meteorol. Soc. Jpn. 77:291–303.
doi:10.2151/jmsj1965.77.1B_291
Doerr, S.H., R.A. Shakesby, and R.P.D. Walsh. 2000. Soil water repellency:
Its causes, characteristics and hydro-geomorphological significance.
Earth Sci. Rev. 51:33–65. doi:10.1016/S0012-8252(00)00011-8
Doerr, S.H. and A.D. Thomas. 2000. The role of soil moisture in controlling
water repellency: New evidence from forest soils in Portugal. J. Hydrol.
231–232:134–147. doi:10.1016/S0022-1694(00)00190-6
d’Orgeval, T., J. Polcher, and P. de Rosnay. 2008. Sensitivity of the West African hydrological cycle in ORCHIDEE to infiltration processes. Hydrol.
Earth Syst. Sci. 12:1387–1401. doi:10.5194/hess-12-1387-2008

p. 45 of 53

Ducharne, A., R.D. Koster, M.J. Suarez, M. Stieglitz, and P. Kumar. 2000.
A catchment-based approach to modeling land surface processes in a general circulation model: 2. Parameter estimation andmodel demonstration. J. Geophys. Res. 105(D20):24823–24838.
doi:10.1029/2000JD900328
Ducharne, A., K. Laval, and J. Polcher. 1998. Sensitivity of the hydrologicalcycle to the parameterization of soil hydrology in a GCM. Clim. Dyn.
14:307–327. doi:10.1007/s003820050226
Ducharne, A., C. Ottlé, F. Maignan, N. Vuichard, J. Ghattas, F. Wang, etal. 2017. The hydrol module of ORCHIDEE: Scientific documentation[rev 3977] and on, work in progress, towards CMIP6v1. Tech. Rep. Inst.
Pierre Simon Laplace, Paris.
Dümenil, L., and E. Todini. 1992. A rainfall–runoff scheme for the use in the
Hamburg climate model. In: J.P. O’Kane, editor, Advances in theoreticalhydrology: A tribute to James Dooge. Elsevier, Amsterdam. p. 129–157.
doi:10.1016/B978-0-444-89831-9.50016-8
Dunkerley, D. 2012. Effects of rainfall intensity fluctuations on infiltration
and runoff: Rainfall simulation on dryland soils, Fowlers Gap, Australia.
Hydrol. Processes 26:2211–2224. doi:10.1002/hyp.8317
Dunkerley, D. 2017. An approach to analysing plot scale infiltration andrunoff responses to rainfall of fluctuating intensity. Hydrol. Processes
31:191–206. doi:10.1002/hyp.10990
Dunne, T. 1978. Field studies of hillslope flow processes and their significance. In: M.J. Kirkby, editor, Hillslope hydrology. Wiley-Interscience,
New York. p. 227–293.
Dunne, T., and R.D. Black. 1970. An experimental investigation of
runoff production in permeable soils. Water Resour. Res. 6:478–490.
doi:10.1029/WR006i002p00478
Edwards, L.M. 1991. The effect of alternate freezing and thawing on aggregate stability and aggregate size distributionof some Prince Edward Island soils. J. Soil Sci. 42:193–204.
doi:10.1111/j.1365-2389.1991.tb00401.x
Egorov, A.G., R.Z. Dautov, J.L. Nieber, and A.Y. Sheshukov. 2003. Stability
analysis of gravity-driven infiltrating flow. Water Resour. Res. 39:1266.
doi:10.1029/2002WR001886
Eigenbrod, K.D. 1996. Effects of cyclic freezing and thawing on volume
changes and permeabilities of soft fine-gained soils. Can. Geotech. J.
33:529–537. doi:10.1139/t96-079-301
Eliassi, M., and R.J. Glass. 2002. On the porous-continuum modeling ofgravity-driven fingers in unsaturated materials: Extension of standard
theory with a hold-back-pile-up effect. Water Resour. Res. 38(11):1234.
doi:10.1029/2001WR001131
Entekhabi, D. and P.S. Eagleson. 1989. Land Surface hydrology parameterization for atmospheric general circulation models
including subgrid scale spatial variability. J. Climate 2:816–831.
doi:10.1175/1520-0442(1989)002<0816:lshpfa>2.0.co;2.
Everett, D.H. 1955. A general approach to hysteresis: 4. An alternative
formulation of the domain model. Trans. Faraday Soc. 51:1551–1557.
doi:10.1039/tf9555101551
Eyring, V., S. Bony, G.A. Meehl, C.A. Senior, B. Stevens, R.J. Stouffer, and K.E.
Taylor. 2016. Overview of the Coupled Model Intercomparison Project
Phase 6 (CMIP6) experimental design and organization. Geosci. Model
Dev. 9:1937–1958. doi:10.5194/gmd-9-1937-2016
FAO. 1988. UNESCO soil map of the world, revised legend. World Resour.
Rep. 60. FAO, Rome.
FAO. 2003. The digitized soil map of the world and derived soil properties.
FAO Land and Water Digital Media Series 1. FAO, Rome.
FAO. 2009. Harmonized world soil database (version 1.1). FAO, Rome.
http://webarchive.iiasa.ac.at/Research/LUC/External-World-soildatabase/HTML/Farrell, D.A., and W.E. Larson. 1972. Dynamics of the soil-water system during
a rainstorm. Soil Sci. 113:88–95. doi:10.1097/00010694-197202000-00003
Feddes, R.A., E. Bresler, and S.P. Neuman. 1974. Field test of a modified
numerical model for water uptake by root systems. Water Resour. Res.
10:1199–1206. doi:10.1029/WR010i006p01199
Flowers, M.D., and R. Lal. 1998. Axle load and tillage effects on
soil physical properties and soybean grain yield on a Mollic Ochraqualf in northwest Ohio. Soil Tillage Res. 48:21–35.
doi:10.1016/S0167-1987(98)00095-6

VZJ | Advancing Critical Zone ScienceFouli, Y., B.J. Cade-Menun, and H.W. Cutforth. 2013. Freeze–thaw cycles
and soil water content effects on infiltration rate of three Saskatchewan soils. Can. J. Soil Sci. 93:485–496. doi:10.4141/cjss2012-060
Franzluebbers, A. 2002. Water infiltration and soil structure related
to organic matter and its stratification with depth. Soil Tillage Res.
66:197–205. doi:10.1016/S0167-1987(02)00027-2
Frauenfeld, B., and C. Truman. 2004. Variable rainfall intensity effects on
runoff and interrill erosion from two Coastal Plain Ultisols in Georgia.
Soil Sci. 169:143–154. doi:10.1097/01.ss.0000117784.98510.46
Freeze, R.A. 1980. A stochastic-conceptual analysis of rainfall–runoff processes on a hillslope. Water Resour. Res. 16:391–408.
doi:10.1029/WR016i002p00391
Frenkel, H., J.O. Goertzen, and J.D. Rhoades. 1978. Effects of clay type andcontent, exchangeable sodium percentage, and electrolyte concentration on clay dispersion and soil hydraulic conductivity. Soil Sci. Soc.
Am. J. 42:32–39. doi:10.2136/sssaj1978.03615995004200010008x
Furman, A., A. Warrick, D. Zerihun, and C. Sanchez. 2006. Modified Kostiakov infiltration function: Accounting for Initial
And Boundary Conditions. J. Irrig. Drain. Eng. 132:587–596.
doi:10.1061/(ASCE)0733-9437(2006)132:6(587)Gao, H., and M. Shao. 2015. Effects of temperature changes on soil hydraulicproperties. Soil Tillage Res. 153:145–154. doi:10.1016/j.still.2015.05.003
Gardner, R. 1955. Relations of temperature to moisture tension of soils.
Soil Sci. 79:257–266. doi:10.1097/00010694-195504000-00003
Gardner, W.R. 1960. Soil water relations in arid and semi-arid conditions. In:
Arid Zone Research XVI: Plant–water relationships in arid and semiarid conditions: Proceedings of the Madrid Symposium, Madrid, Spain.
24–30 Sept. 1959. UNESCO, Paris.
Garrote, L., and R.L. Bras. 1995. A distributed model for real-time flood
forecasting using digital elevation models. J. Hydrol. 167:279–306.
doi:10.1016/0022-1694(94)02592-Y
Gedney, N. and P.M. Cox. 2003. The sensitivity of global climate model simulations to the representationof soil moisture heterogeneity. J. Hydrometeorol. 4:1265–1275.
doi:10.1175/1525-7541(2003)004<1265:Tsogcm>2.0.Co;2
Getirana, A., A. Boone, C. Peugeot, S. Boussetta, S. Ait-Mesbah, M.
Anderson, et al. 2017. Streamflows over a West African basin fromthe ALMIP-2 model ensemble. J. Hydrometeorol. 18:1831–1845.
doi:10.1175/JHM-D-16-0233.1
Getirana, A.C.V., E. Dutra, M. Guimberteau, J. Kam, H.-Y. Li, B. Decharme,
et al. 2014. Water balance in the Amazon basin from a land
surface model ensemble. J. Hydrometeorol. 15:2586–2614.
doi:10.1175/JHM-D-14-0068.1
Ghanbarian, B., V. Taslimitehrani, G.Z. Dong, and Y.A. Pachepsky. 2015.
Sample dimensions effect on prediction of soil water retention
curve and saturated hydraulic conductivity. J. Hydrol. 528:127–137.
doi:10.1016/j.jhydrol.2015.06.024
Giraldez, J.V., and G. Sposito. 1985. Infiltration in swelling soils. Water
Resour. Res. 21:33–44. doi:10.1029/WR021i001p00033
Glass, R.J., T.S. Steenhuis, and J.-Y. Parlange. 1989. Mechanism for finger persistence in homogeneous, unsaturated, porous media: Theory and verification. Soil Sci. 148:60–70. doi:10.1097/00010694-198907000-00007
Gleeson, T., L. Smith, N. Moosdorf, J. Hartmann, H.H. Durr, A.H. Manning, etal. 2011. Mapping permeability over the surface of the Earth. Geophys.
Res. Lett. 38:L02401. doi:10.1029/2010gl045565
Global Soil Data Task. 2014. Global soil data products CD-ROM contents(IGBP-DIS). [Dataset.] Oak Ridge Nat. Lab. Distributed Active Archive
Ctr., Oak Ridge, TN. http://dx.doi.org/10.3334/ORNLDAAC/565.
Govindaraju, R.S., C. Corradini, and R. Morbidelli. 2006. A semi-analyticalmodel of expected areal-average infiltration under spatial heterogeneity of rainfall and soil saturated hydraulic conductivity. J. Hydrol.
316:184–194. doi:10.1016/j.jhydrol.2005.04.019
Govindaraju, R.S., C. Corradini, and R. Morbidelli. 2012. Local- and fieldscale infiltration into vertically non-uniform soils with spatially-variable surface hydraulic conductivities. Hydrol. Processes 26:3293–3301.
doi:10.1002/hyp.8454
Govindaraju, R.S., R. Morbidelli, and C. Corradini. 2001. Areal infiltration
modeling over soils with spatially correlated hydraulic conductivities.
J. Hydrol. Eng. 6:150–158. doi:10.1061/(ASCE)1084-0699(2001)6:2(150)

p. 46 of 53

Grant, S.A., and J. Bachmann. 2002. Effect of temperature on capillarypressure. In: P.A.C. Raats et al., editors, Environmental mechanics, water, mass and energy transfer in the biosphere. Geophys.
Monogr. Ser. 129. Am. Geophys. Union, Washington, DC. p. 199–212.
doi:10.1029/129GM18
Grant, S.A., and A. Salehzadeh. 1996. Calculation of temperature effects onwetting coefficients of porous solids and their capillary pressure functions. Water Resour. Res. 32:261–270. doi:10.1029/95WR02915
Green, W.H., and G. Ampt. 1911. Studies on soil physics. J. Agric. Sci. 4:1–24.
doi:10.1017/S0021859600001441
Greene, R.S.B., and P.B. Hairsine. 2004. Elementary processes of soil–waterinteraction and thresholds in soil surface dynamics: A review. Earth
Surf. Processes Landforms 29:1077–1091. doi:10.1002/esp.1103
Gudmundsson, L., L.M. Tallaksen, K. Stahl, D.B. Clark, E. Dumont, S.
Hagemann, et al. 2012. Comparing large-scale hydrological model
simulations to observed runoff percentiles in Europe. J. Hydrometeorol. 13:604–620. doi:10.1175/JHM-D-11-083.1
Habets, F., J. Noilhan, C. Golaz, J.P. Goutorbe, P. Lacarrère, E. Leblois,
et al. 1999. The ISBA surface scheme in a macroscale hydrological model applied to the Hapex-Mobilhy area: II. Simulationof streamflows and annual water budget. J. Hydrol. 217:97–118.
doi:10.1016/S0022-1694(99)00020-7
Habets, F., and G.M. Saulnier. 2001. Subgrid runoff parameterization. Phys.
Chem. Earth, Part B 26:455–459. doi:10.1016/S1464-1909(01)00034-X
Hagemann, S., and L.D. Gates. 2003. Improving a subgrid runoff parameterization scheme for climate models by the use of high resolution
data derived from satellite observations. Clim. Dyn. 21:349–359.
doi:10.1007/s00382-003-0349-x
Hagemann, S., and T. Stacke. 2015. Impact of the soil hydrology schemeon simulated soil moisture memory. Clim. Dyn. 44:1731–1750.
doi:10.1007/s00382-014-2221-6
Hailegeorgis, T.T., Y.S. Abdella, K. Alfredsen, and S. Kolberg. 2015. Evaluation of regionalization methods for hourly continuous streamflow
simulation using distributed models in boreal catchments. J. Hydrol.
Eng. 20(11). doi:10.1061/(ASCE)HE.1943-5584.0001218
Haines, W.B. 1930. Studies in the physical properties of soil: V. The
hysteresis effect in capillary properties, and the modes of moisture distribution associated therewith. J. Agric. Sci. 20:97–116.
doi:10.1017/S002185960008864X

equation based on physical parameters. Soil Sci. 149:292–300.
doi:10.1097/00010694-199005000-00006
Healy, R.W. 2008. Simulating water, solute, and heat transport in the subsurface with the VS2DI software package. Vadose Zone J. 7:632–639.
doi:10.2136/vzj2007.0075
Hengl, T., J.M. de Jesus, G.B.M. Heuvelink, M.R. Gonzalez, M. Kilibarda, A. Blagotic, et al. 2017. SoilGrids250m: Global gridded
soil information based on machine learning. PLoS One 12.
doi:10.1371/journal.pone.0169748
Hengl, T., J.M. de Jesus, R.A. MacMillan, N.H. Batjes, G.B.M. Heuvelink, E.
Ribeiro, et al. 2014. SoilGrids1km: Global soil information based on
automated mapping. PLoS One 9. doi:10.1371/journal.pone.0105992
Herbst, M., B. Diekkrüger, and J. Vanderborght. 2006. Numericalexperiments on the sensitivity of runoff generation to the spatial variation of soil hydraulic properties. J. Hydrol. 326:43–58.
doi:10.1016/j.jhydrol.2005.10.036
Hillel, D. 1980. Fundamentals of soil physics. Academic Press, San Diego.
Hillel, D. 1998. Environmental soil physics. Academic Press, San Diego.
Hillel, D., and W.R. Gardner. 1969. Steady infiltration into crust topped profiles. Soil Sci. 108:137–142. doi:10.1097/00010694-196908000-00010
Hillel, D., and W.R. Gardner. 1970. Transient infiltration into crust topped
profiles. Soil Sci. 109:69–76. doi:10.1097/00010694-197002000-00001
Hirmas, D.R., D. Giménez, A. Nemes, R. Kerry, N.A. Brunsell, and C.J.
Wilson. 2018. Climate-induced changes in continental-scale soil
macroporosity may intensify water cycle. Nature 516:100–103.
doi:10.1038/s41586-018-0463-x
Hogarth, W.L., J. Hopmans, J.-Y. Parlange, and R. Haverkamp. 1988. Application of a simple soil-water hysteresis model. J. Hydrol. 98:21–29.
doi:10.1016/0022-1694(88)90203-X
Hogue, T.S., L.A. Bastidas, H.V. Gupta, and S. Sorooshian. 2006. Evaluating model performance and parameter behavior for varying levelsof land surface model complexity. Water Resour. Res. 42:W08430.
doi:10.1029/2005WR004440
Hopmans, J. 1989. Stochastic description of field-measured infiltration
data. Trans. ASAE 32:1987–1993. doi:10.13031/2013.31252
Hopmans, J.W., and J.H. Dane. 1986. Temperature dependenceof soil hydraulic properties. Soil Sci. Soc. Am. J. 50:4–9.
doi:10.2136/sssaj1986.03615995005000010001x

Hamza, M.A., and W.K. Anderson. 2005. Soil compaction in cropping systems: A review of the nature, causes and possible solutions. Soil Tillage
Res. 82:121–145. doi:10.1016/j.still.2004.08.009

Hopmans, J.W., J.-Y. Parlange, and S. Assouline. 2007. Infiltration. In: J.W.
Delleur, editor The handbook of groundwater engineering. 2nd ed.
CRC Press, Boca Raton, FL.

Hanks, R.J., and S.A. Bowers. 1962. Numerical solution of the moistureflow equation for infiltration into layered soils. Soil Sci. Soc. Am. Proc.
26:530–534. doi:10.2136/sssaj1962.03615995002600060007x

Horton, R.E. 1933. The role of infiltration in the hydrologic cycle. Eos Trans.
AGU 14:446–460.

Hanks, R.J., A. Klute, and E. Bresler. 1969. A numeric method for estimating
infiltration, redistribution, drainage, and evaporation of water fromsoil. Water Resour. Res. 5:1064–1069. doi:10.1029/WR005i005p01064
Hannes, M., U. Wollschläger, T. Wöhling, and H.-J. Vogel. 2016. Revisiting hydraulic hysteresis based on long-term monitoring of
hydraulic states in lysimeters. Water Resour. Res. 52:3847–3865.
doi:10.1002/2015WR018319

Horton, R.E. 1940. The infiltration‐theory of surface‐runoff. Eos Trans. AGU
21:541. doi:10.1029/TR021i002p00541-1
Horton, R.E. 1941. An approach toward a physical interpretation of infiltration-capacity. Soil Sci. Soc. Am. J. 5(C):399–417.
doi:10.2136/sssaj1941.036159950005000C0075x
Huang, H.-C., Y.-C. Tan, C.-W. Liu, and C.-H. Chen. 2005. A novel hysteresis model in unsaturated soil. Hydrol. Processes 19:1653–1665.
doi:10.1002/hyp.5594

Hansen, S., P. Abrahamsen, C.T. Petersen, and M. Styczen. 2012. Daisy:
Model use, calibration, and validation. Trans. ASABE 55:1315–1333.
doi:10.13031/2013.42244

Huang, M., X. Liang, and Y. Liang. 2003. A transferability study of model
parameters for the variable infiltration capacity land surface scheme. J.
Geophys. Res. Atmos. 108(D22):8864. doi:10.1029/2003JD003676

Haverd, V., and M. Cuntz. 2010. Soil-Litter-Iso: A one-dimensionalmodel for coupled transport of heat, water and stable isotopes insoil with a litter layer and root extraction. J. Hydrol. 388:438–455.
doi:10.1016/j.jhydrol.2010.05.029

Huang, S.C., R. Kumar, M. Florke, T. Yang, Y. Hundecha, P. Kraft, et al.
2017. Evaluation of an ensemble of regional hydrological models
in 12 large-scale river basins worldwide. Clim. Change 141:381–397.
doi:10.1007/s10584-016-1841-8

Haverd, V., M. Cuntz, L.P. Nieradzik, and I.N. Harman. 2016. Improved
representations of coupled soil–canopy processes in the CABLE
land surface model (Subversion revision 3432). Geosci. Model Dev.
9:3111–3122. doi:10.5194/gmd-9-3111-2016

Hugelius, G., C. Tarnocai, G. Broll, J.G. Canadell, P. Kuhry, and D.K.
Swanson. 2013. The Northern Circumpolar Soil Carbon Database:Spatially distributed datasets of soil coverage and soil carbon storage in the northern permafrost regions. Earth Syst. Sci. Data 5:3–13.
doi:10.5194/essd-5-3-2013

Haverkamp, R., J.-Y. Parlange, R. Cuenca, P.J. Ross, and T.S. Steenhuis. 1998.
Scaling of the Richards equation and its application to watershed
modelling. In: G. Sposito, editor, Scale dependence and scale invariance in hydrology. Cambridge Univ. Press, Cambridge, UK. p. 190–223.
doi:10.1017/CBO9780511551864.008

Huyakorn, P.S., E.P. Springer, V. Guvanasen, and T.D. Wadsworth. 1986. A
three-dimensional finite-element model for simulating water flow invariably saturated porous media. Water Resour. Res. 22:1790–1808.
doi:10.1029/WR022i013p01790

Haverkamp, R., J.-Y. Parlange, J. Starr, G. Schmitz, and C. Fuentes. 1990. Infiltration under ponded conditions: 3. A predictive

Ibrahim, H.A., and W. Brutsaert. 1968. Intermittent infiltration into soils withhysteresis. J. Hydraul. Div., Proc. Am. Soc. Civ. Eng. 94(HY1):113–137.

VZJ | Advancing Critical Zone Science

p. 47 of 53

IGBP. 2000. Global Soil Data Task (IGBP-DIS, ISO-image of CD). International Geosphere–Biosphere Program, PANGAEA, Alfred Wegener Inst.,
Univ. of Bremen, Bremen, Germany. doi:10.1594/PANGAEA.869912
Imeson, A.C., J.M. Verstraten, E.J. van Mulligen, and J. Sevink. 1992.
The effects of fire and water repellency on infiltration andrunoff under Mediterranean type forest. Catena 19:345–361.
doi:10.1016/0341-8162(92)90008-Y
Jame, Y.-W. 1977. Heat and mass transfer in freezing unsaturated soil. Univ.
of Saskatchewan, Saskatoon, SK.
Jana, R., and B. Mohanty. 2011. Enhancing PTFs with remotely sensed data
for multi-scale soil water retention estimation. J. Hydrol. 399:201–211.
doi:10.1016/j.jhydrol.2010.12.043
Jana, R., and B. Mohanty. 2012a. On topographic controls of soil hydraulic parameter scaling at hillslope scales. Water Resour. Res. 48:2518.
doi:10.1029/2011WR011204
Jana, R.B., and B.P. Mohanty. 2012b. A comparative study of multiple approaches to soil hydraulic parameter scaling applied at the hillslope
scale. Water Resour. Res. 48:W02520. doi:10.1029/2010WR010185
Jana, R.B., and B.P. Mohanty. 2012c. A topography-based scaling algorithm
for soil hydraulic parameters at hillslope scales: Field testing. Water
Resour. Res. 48:W02519. doi:10.1029/2011WR011205
Jana, R.B., B.P. Mohanty, and E.P. Springer. 2008. Multiscale Bayesianneural networks for soil water content estimation. Water Resour. Res.
44:W08408. doi:10.1029/2008WR006879
Jansson, P.E. 2012. CoupModel: Model use, calibration, and validation.
Trans. ASABE 55:1337–1344. doi:10.13031/2013.42245
Jarvis, N.J. 2007. A review of non-equilibrium water flow and solute transport in soil macropores: Principles, controlling factors
and consequences for water quality. Eur. J. Soil Sci. 58:523–546.
doi:10.1111/j.1365-2389.2007.00915.x
Jorda, H., M. Bechtold, N. Jarvis, and J. Koestel. 2015. Using boosted
regression trees to explore key factors controlling saturated andnear-saturated hydraulic conductivity. Eur. J. Soil Sci. 66:744–756.
doi:10.1111/ejss.12249
Jury, W.A., W.R. Gardner, and W.H. Gardner. 1991. Soil physics. 5th ed. John
Wiley & Sons, New York.
Jury, W.A., Z. Wang, and A. Tuli. 2003. A conceptual model of unstableflow in unsaturated soil during redistribution. Vadose Zone J. 2:61–67.
doi:10.2136/vzj2003.6100
Kale, R., and B. Sahoo. 2011. Green–Ampt infiltration models for varied
field conditions: A revisit. Water Resour. Manage. 25:3505–3536.
doi:10.1007/s11269-011-9868-0
Keune, J., F. Gasper, K. Goergen, A. Hense, P. Shrestha, M. Sulis, and S.
Kollet. 2016. Studying the influence of groundwater representations on land surface-atmosphere feedbacks during the European
heat wave in 2003. J. Geophys. Res. Atmos. 121:13301–13325.
doi:10.1002/2016JD025426
Kim, D., R.L. Ray, and M. Choi. 2017. Simulations of energy balance components at snow-dominated montane watershed by land surface models.
Environ. Earth Sci. 76:337. doi:10.1007/s12665-017-6655-0
Kim, J., and B. Mohanty. 2017. A physically-based hydrological connectivity algorithm for describing spatial patterns of soil moisture in theunsaturated zone: Hydrologic connectivity in soil moisture. J. Geophys.
Res. Atmos. 122:2096–2114. doi:10.1002/2016JD025591
Kim, K.H., and W.P. Miller. 1996. Effect of rainfall electrolyte concentration and slope on infiltration and erosion. Soil Technol. 9:173–185.
doi:10.1016/S0933-3630(96)00011-6
King, P.M. 1981. Comparison of methods for measuring severity of water
repellence of sandy soils and assessment of some factors that affect its
measurement. Aust. J. Soil Res. 19:275–285. doi:10.1071/SR9810275
Kishné, A.S., C.L.S. Morgan, Y. Ge, and W.L. Miller. 2010. Antecedent soil
moisture affecting surface cracking of a Vertisol in field conditions.
Geoderma 157:109–117. doi:10.1016/j.geoderma.2010.03.020
Kollet, S.J., and R.M. Maxwell. 2006. Integrated surface-groundwaterflow modeling: A free-surface overland flow boundary condition ina parallel groundwater flow model. Adv. Water Resour. 29:945–958.
doi:10.1016/j.advwatres.2005.08.006

grated, distributed watershed model. Water Resour. Res. 44:W02402.
doi:10.1029/2007wr006004
Kollet, S.J., and R.M. Maxwell. 2008b. Demonstrating fractal scalingof baseflow residence time distributions using a fully-coupled
groundwater and land surface model. Geophys. Res. Lett. 35:L07402.
doi:10.1029/2008GL033215
Koren, V., J. Schaake, K. Mitchell, Q.Y. Duan, F. Chen, and J.M. Baker. 1999. A
parameterization of snowpack and frozen ground intended for NCEP
weather and climate models. J. Geophys. Res. Atmos. 104:19569–
19585. doi:10.1029/1999JD900232
Koster, R.D., M.J. Suarez, A. Ducharne, M. Stieglitz, and P. Kumar. 2000. A
catchment-based approach to modeling land surface processes in a
general circulation model: 1. Model structure. J. Geophys. Res. Atmos.
105:24809–24822. doi:10.1029/2000JD900327
Kostiakov, A.N. 1932. On the dynamics of the coefficients of water percolation in soils and on the necessity of studying it from a dynamic pointof view for purpose of amelioration. In: O. Fauser, editor, Transactionsof the 6th Commission of the International Society of Soil Science,
Russia. Vol. A. Martinus Nijhoff, Groningen, the Netherlands. p. 17–21.
Kowalczyk, E.A., L. Stevens, R.M. Law, M.R. Dix, Y.P. Wang, I.N. Harman, et al.
2013. The land surface model component of ACCESS: Description andimpact on the simulated surface climatology. Aust. Meteorol. Oceanogr. J. 63:65–82. doi:10.22499/2.6301.005
Kowalczyk, E.A., Y.P. Wang, R.M. Law, H.L. Davies, J.L. McGregor, and G.
Abramowitz. 2006. The CSIRO Atmosphere Biosphere Land Exchange(CABLE) model for use in climate models and as an offline model. Pap.
013. CSIRO Mar. Atmos. Res., Aspendale, VIC, Australia.
Kutilek, M., V. Kuraz, and M. Krejca. 1993. Measurement time and spatial
variability of field infiltration. Int. Agrophys. 7:133–140.
Kværnø, S.H., and L. Øygarden. 2006. The influence of freeze–thaw cycles
and soil moisture on aggregate stability of three soils in Norway. Catena 67:175–182. doi:10.1016/j.catena.2006.03.011
Laliberte, G.E., A.T. Corey, and R.H. Brooks. 1966. Properties of unsaturated
porous media. Hydrol. Pap. 17. Colorado State Univ., Fort Collins.
Lange, B., P. Lüescher, and P.F. Germann. 2009. Significance of tree roots
for preferential infiltration in stagnic soils. Hydrol. Earth Syst. Sci.
13:1809–1821. doi:10.5194/hess-13-1809-2009
Lawrence, D.M., K.W. Oleson, M.G. Flanner, P.E. Thornton, S.C. Swenson, P.J.
Lawrence, et al. 2011. Parameterization improvements and functionaland structural advances in Version 4 of the Community Land Model. J.
Adv. Model. Earth Syst. 3:M03001. doi:10.1029/2011MS00045
Lawrence, D.M., and A.G. Slater. 2008. Incorporating organic soil into a global
climate model. Clim. Dyn. 30:145–160. doi:10.1007/s00382-007-0278-1
Lehmann, P., and D. Or. 2012. Hydromechanical triggering of landslides:
From progressive local failures to mass release. Water Resour. Res.
48:W03535. doi:10.1029/2011WR010947
Lehrsch, G.A. 1998. Freeze–thaw cycles increase near-surface aggregatestability. Soil Sci. 163:63–70. doi:10.1097/00010694-199801000-00009
Lehrsch, G.A., R.E. Sojka, D.L. Carter, and P.M. Jolley. 1991. Freezing effects on aggregate stability affected by texture, mineralogy, and organic matter. Soil Sci. Soc. Am. J. 55:1401–1406.
doi:10.2136/sssaj1991.03615995005500050033x
Letts, M.G., N.T. Roulet, N.T. Comer, M.R. Skarupa, and D.L. Verseghy.
2000. Parametrization of peatland hydraulic properties for the
Canadian Land Surface Scheme. Atmos.-Ocean 38:141–160.
doi:10.1080/07055900.2000.9649643
Le Vine, N., A. Butler, N. McIntyre, and C. Jackson. 2016. Diagnosinghydrological limitations of a land surface model: Application of JULES
to a deep-groundwater chalk basin. Hydrol. Earth Syst. Sci. 20:143–159.
doi:10.5194/hess-20-143-2016
Levy, G.J., H.J.C. Smith, and M. Agassi. 1989. Water temperature effect on
hydraulic conductivity and infiltration rate of soils. S. Afr. J. Plant Soil
6:240–244. doi:10.1080/02571862.1989.10634520
Lewis, M. 1937. The rate of infiltration of water in irrigation‐practice. EosTrans. AGU 18:361–368.
Li, Q., S.F. Sun, and Y.K. Xue. 2010. Analyses and development of a hierarchy of frozen soil models for cold region study. J. Geophys. Res. Atmos.
115:D03107. doi:10.1029/2009jd012530

Kollet, S.J., and R.M. Maxwell. 2008a. Capturing the influence ofgroundwater dynamics on land surface processes using an inte-

VZJ | Advancing Critical Zone Sciencep. 48 of 53

Liang, X., D.P. Lettenmaier, E.F. Wood, and S.J. Burges. 1994. A simple hydrologically based model of land surface water and energy fluxes for
GSMs. J. Geophys. Res. 99(D7):14415–14428. doi:10.1029/94JD00483
Liang, X., and Z. Xie. 2001. A new surface runoff parameterization with
subgrid-scale soil heterogeneity for land surface models. Adv. Water
Resour. 24:1173–1193. doi:10.1016/S0309-1708(01)00032-X
Lin, H.S., K.J. Mclnnes, L.P. Wilding, and C.T. Hallmark. 1998. Macroporosity
and initial moisture effects on infiltration rates in Vertisols and vertic
intergrades. Soil Sci. 163:2–8. doi:10.1097/00010694-199801000-00002
Loague, K., and G.A. Gander. 1990. R-5 revisited: 1. Spatial variabilityof infiltration on a small rangeland catchment. Water Resour. Res.
26:957–971. doi:10.1029/WR026i005p00957
Logsdon, S.D., and D.B. Jaynes. 1996. Spatial variability of hydraulic conductivity in a cultivated field at different times. Soil Sci. Soc. Am. J. 60.
doi:10.2136/sssaj1996.03615995006000030003x
MacDonald, M.K., J.W. Pomeroy, and R.L.H. Essery. 2017. Water andenergy fluxes over northern prairies as affected by chinook winds
and winter precipitation. Agric. For. Meteorol. 2017:372–385.
doi:10.1016/j.agrformet.2017.10.025
Mantoglou, A., and L.W. Gelhar. 1987a. Capillary tension head variance,
mean soil moisture content, and effective specific soil moisture capacity of transient unsaturated flow in stratified soils. Water Resour. Res.
23:47–56. doi:10.1029/WR023i001p00047
Mantoglou, A., and L.W. Gelhar. 1987b. Effective hydraulic conductivityof transient unsaturated flow in stratified soils. Water Resour. Res.
23:57–67. doi:10.1029/WR023i001p00057
Materia, S., P.A. Dirmeyer, Z.C. Guo, A. Alessandri, and A. Navarra. 2010.
The sensitivity of simulated river discharge to land surface representation and meteorological forcings. J. Hydrometeorol. 11:334–351.
doi:10.1175/2009JHM1162.1
McGhie, D.A., and A.M. Posner. 1981. The effect of plant top material on
the water repellence of fired sands and water repellent soils. Aust. J.
Agric. Res. 32:609–620. doi:10.1071/AR9810609
McNeal, B.L., and N.T. Coleman. 1966. Effect of solution composition
on soil hydraulic conductivity. Soil Sci. Soc. Am. J. 30:308–312.
doi:10.2136/sssaj1966.03615995003000030007x
Mein, R.G., and C.L. Larson. 1973. Modeling infiltration during a steadyrain. Water Resour. Res. 9:384–394. doi:10.1029/WR009i002p00384
Meng, H., J.D. Salas, T. Green, and L. Ahuja. 2006. Scaling analysis of space–time infiltration based on the universal multifractal model. J. Hydrol.
322:220–235. doi:10.1016/j.jhydrol.2005.03.016
Mertens, J., D. Raes, and J. Feyen. 2002. Incorporating rainfall intensity
into daily rainfall records for simulating runoff and infiltration into soil
profiles. Hydrol. Processes 16:731–739. doi:10.1002/hyp.1005
Mezencev, V. 1948. Theory of formation of the surface runoff. Meteorol.
Hidrol. 3:33–40.
Miller, D.E., and W.H. Gardner. 1962. Water infiltration into stratified soil. Soil Sci.
Soc. Am. J. 26:115–119. doi:10.2136/sssaj1962.03615995002600020007x
Miller, E.E., and R.D. Miller. 1956. Physical theory for capillary flow phenomena. J. Appl. Phys. 27:324–332. doi:10.1063/1.1722370
Milly, P.C.D., and P.S. Eagleson. 1988. Effect of storm scale on surface runoff
volume. Water Resour. Res. 24:620–624. doi:10.1029/WR024i004p00620
Milovac, J., J. Ingwersen, and K. Warrach-Sagi. 2014. Soil texture forcing
data for the whole world for the Weather Research and Forecasting(WRF) model of the University of Hohenheim (UHOH) based on the
Harmonized World Soil Database (HWSD) at 30 arc-second horizontal resolution. World Data Ctr. for Climate, German Climate Computing Ctr., Hamburg.
Mis, J. 1982. Formulation and solution of fundamental problems of vertical infiltration. Vodohosp. Cas. 30:304–311.
Mohammed, A., B. Kurylyk, E. Cey, and M. Hayashi. 2018. Snowmelt
infiltration and macropore flow in frozen soils: Overview, knowledge and a conceptual framework. Vadose Zone J. 17:180084.
doi:10.2136/vzj2018.04.0084
Mohanty, B.P. 1999. Scaling hydraulic properties of a macroporous soil.
Water Resour. Res. 35:1927–1931. doi:10.1029/1999WR900050
Mohanty, B.P., and J. Zhu. 2007. Effective hydraulic parameters inhorizontally and vertically heterogeneous soils for steady-state

VZJ | Advancing Critical Zone Science

land–atmosphere interaction. J. Hydrometeorol. 8:715–729.
doi:10.1175/JHM606.1
Montzka, C., M. Herbst, L. Weihermüller, A. Verhoef, and H. Vereecken.
2017. A global data set of soil hydraulic properties and sub-grid variability of soil water retention and hydraulic conductivity curves. Earth
Syst. Sci. Data 9:529–543. doi:10.5194/essd-9-529-2017
Moore, I.D. 1981a. Effect of surface sealing on infiltration. Trans. ASAE
24:1546–1552. doi:10.13031/2013.34488
Moore, I.D. 1981b. Infiltration equations modified for surface effect. J. Irrig.
Drain. Div., Am. Soc. Civ. Eng. 107:71–86.
Moore, R.J. 1985. The probability-distributed principle and runoff
production at point and basin scales. Hydrol. Sci. J. 30:273–297.
doi:10.1080/02626668509490989
Morbidelli, R., C. Corradini, and R.S. Govindaraju. 2006. A field-scale
infiltration model accounting for spatial heterogeneity of rainfall andsoil saturated hydraulic conductivity. Hydrol. Processes 20:1465–1481.
doi:10.1002/hyp.5943
Morbidelli, R., C. Saltalippi, A. Flammini, and R.S. Govindaraju. 2018.
Role of slope on infiltration: A review. J. Hydrol. 557:878–886.
doi:10.1016/j.jhydrol.2018.01.019
Mualem, Y. 1974. A conceptual model of hysteresis. Water Resour. Res.
10:514–520. doi:10.1029/WR010i003p00514
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
doi:10.1029/WR012i003p00513
Mualem, Y., and S. Assouline. 1989. Modeling soil seal as a nonuniform layer. Water Resour. Res. 25:2101–2108. doi:10.1029/WR025i010p02101
Mualem, Y., and S. Assouline. 1996. Soil sealing, infiltration and runoff. In:
A.S. Issar and S.D. Resnick, editors, Runoff, infiltration and subsurface
flow of water in arid and semi-arid regions. Water Science and Technology Library. Kluwer Acad., Dordrecht, the Netherlands. p. 131–181.
doi:10.1007/978-94-017-2929-1_4
Mualem, Y., and S. Assouline. 1992. Flow processes in sealing soils: Conceptions and solutions. In: M.E. Sumner and B.A. Stewart, editors, Soil
crusting: Chemical and physical processes. Lewis Publ., Boca Raton, FL.
Mualem, Y., and A. Beriozkin. 2009. General scaling rules of the hysteretic
water retention function based on Mualem’s domain theory. Eur. J. Soil
Sci. 60:652–661. doi:10.1111/j.1365-2389.2009.01130.x
Mualem, Y., and G. Dagan. 1975. A dependent domain modelof capillary hysteresis. Water Resour. Res. 11:452–460.
doi:10.1029/WR011i003p00452
Mueller, A., E. Dutra, H. Cloke, A. Verhoef, G. Balsamo, and F. Pappenberger.
2016. Water infiltration and redistribution in land surface models. Eur.
Ctr. for Medium-Range Weather Forecasts Tech. Mem. Research Dep.,
Univ. of Reading, Reading, UK.
Neuman, S.P. 1973. Saturated–unsaturated seepage by finite elements. J.
Hydraul. Div., Am. Soc. Civ. Eng. 99:2233–2250.
Neuman, S.P. 1976. Wetting front pressure head in the infiltration model of Green and Ampt. Water Resour. Res. 12:564–566.
doi:10.1029/WR012i003p00564
Neuman, S.P., R.A. Feddes, and E. Bresler. 1974. Finite element simulation offlow in saturated–unsaturated soils considering water uptake by plants.
3rd Annu. Rep. Hydraulic Engineering Lab., Technion, Haifa, Israel.
Nielsen, D., J. Biggar, and K. Erh. 1973. Spatial variability of field-measured soilwater properties. Hilgardia 42:215–259. doi:10.3733/hilg.v42n07p215
Nielsen, D.R., J.W. Hopmans, and K. Reichardt. 1998. An emerging technology for scaling field soil-water behavior. In: G. Sposito, editor, Scaledependence and scale invariance in hydrology. Cambridge Univ. Press,
Cambridge, UK. p. 136–166. doi:10.1017/CBO9780511551864.006
Nimmo, J.R. 1992. Semiempirical model of soil water hysteresis. Soil Sci. Soc.
Am. J. 56:1723–1730. doi:10.2136/sssaj1992.03615995005600060011x
Nimmo, J.R., and E.E. Miller. 1986. The temperature dependence of isothermal moisture vs. potential characteristics of soils. Soil Sci. Soc. Am.
J. 50:1105–1113. doi:10.2136/sssaj1986.03615995005000050004x
Nimmo, J.R., K.S. Perkins, K.M. Schmidt, D. Miller, J.D. Stock, and K. Singha.
2009. Hydrologic characterization of desert soils with varying degrees
of pedogenesis: 1. Field experiments evaluating plant-relevant soil
water behavior. Vadose Zone J. 8:480–495. doi:10.2136/vzj2008.0052

p. 49 of 53

Nissen, T.M., and M.M. Wander. 2003. Management and soil-quality effects
on fertilizer-use efficiency and leaching. Soil Sci. Soc. Am. J. 67:1524–
1532. doi:10.2136/sssaj2003.1524
Niu, G.-Y., and Z.-L. Yang. 2007. An observation-based formulation of snow cover fraction and its evaluation over large North
American river basins. J. Geophys. Res. D Atmos. 112: D21101.
doi:10.1029/2007JD008674
Niu, G.Y., Z.L. Yang, R.E. Dickinson, and L.E. Gulden. 2005. A simple TOPMODEL-based runoff parameterization (SIMTOP) for use in global climate models. J. Geophys. Res. 110:D21106. doi:10.1029/2005JD006111
Niu, G.-Y., Z.-L. Yang, K.E. Mitchell, F. Chen, M.B. Ek, M. Barlage, et al. 2011. The
community Noah land surface model with multiparameterization options(Noah-MP): 1. Model description and evaluation with local-scale measurements. J. Geophys. Res. 116:D12109. doi:10.1029/2010JD015139
Noilhan, J., and P. Lacarrère. 1995. GCM gridscale evaporation frommesoscale modelling. J. Clim. 8:206–223. doi:10.1175/15200442(1995)008<0206:GGSEFM>2.0.CO;2
Nuth, M., and L. Laloui. 2008. Advances in modelling hysteretic water
retention curve in deformable soils. Comput. Geotech. 35:835–844.
doi:10.1016/j.compgeo.2008.08.001
Ogden, F.L., W. Lai, R.C. Steinke, J. Zhu, C.A. Talbot, and J.L. Wilson. 2015.
A new general 1-D vadose zone flow solution method. Water Resour.
Res. 51:4282–4300. doi:10.1002/2015WR017126
Oleson, K.W., D.M. Lawrence, G.B. Bonan, B. Drewniak, M. Huang, C.D.
Koven, et al. 2013. Technical description of version 4.5 of the Community Land Model (CLM). NCAR Tech. Note NCAR/TN-503+STR. Natl. Ctr.
for Atmos. Res., Boulder, CO.
Oleson, K.W., G.-Y. Niu, Z.-L. Yang, D.M. Lawrence, P.E. Thornton, P.J. Lawrence, et al. 2008. Improvements to the Community Land Model andtheir impact on the hydrological cycle. J. Geophys. Res. 113:G01021.
doi:10.1029/2007JG000563
Or, D. 2008. Scaling of capillary, gravity and viscous forces affecting
flow morphology in unsaturated porous media. Adv. Water Resour.
31:1129–1136. doi:10.1016/j.advwatres.2007.10.004
Or, D., F.J. Leij, V. Snyder, and T.A. Ghezzehei. 2000. Stochastic model for
posttillage soil pore space evolution. Water Resour. Res. 36:1641–1652.
doi:10.1029/2000WR900092
Pachepsky, Y., and R.L. Hill. 2017. Scale and scaling in soils. Geoderma
287:4–30. doi:10.1016/j.geoderma.2016.08.017
Paniconi, C., and M. Putti. 1994. A comparison of Picard and Newton
iteration in the numerical solution of multidimensional variably saturated flow problems. Water Resour. Res. 30:3357–3374.
doi:10.1029/94WR02046
Paniconi, C., and M. Putti. 2015. Physically based modeling in catchmenthydrology at 50: Survey and outlook. Water Resour. Res. 51:7090–7129.
doi:10.1002/2015WR017780
Parhi, P.K., S.K. Mishra, and R. Singh. 2007. A modification to Kostiakov
and Modified Kostiakov infiltration models. Water Resour. Manage.
21:1973–1989. doi:10.1007/s11269-006-9140-1
Park, S.-E., A. Bartsch, D. Sabel, W. Wagner, V. Naeimi, and Y. Yamaguchi. 2011. Monitoring freeze/thaw cycles using ENVISAT ASAR Global Mode. Remote Sens. Environ. 115:3457–3467.
doi:10.1016/j.rse.2011.08.009
Parlange, J.-Y. 1976. Capillary hysteresis and relationship between
drying and wetting curves. Water Resour. Res. 12:224–228.
doi:10.1029/WR012i002p00224
Parlange, J.-Y. 1980. Water transport in soils. Annu. Rev. Fluid Mech.
12:77–102. doi:10.1146/annurev.fl.12.010180.000453
Parlange, J.-Y., R. Haverkamp, and J. Touma. 1985. Infiltration under
ponded conditions: 1. Optimal analytical solution and comparison with experimental observations. Soil Sci. 139:305–311.
doi:10.1097/00010694-198504000-00003
Parlange, J.-Y., and D.E. Hill. 1976. Theoretical analysis of wetting front instability in soils. Soil Sci. 122:236–239. doi:10.1097/00010694-197610000-00008
Parlange, J.-Y., W.L. Hogarth, and M.B. Parlange. 1984. Optimal analysis of the effect of a surface crust. Soil Sci. Soc. Am. J. 48:494–497.
doi:10.2136/sssaj1984.03615995004800030004x
Parlange, J.-Y., I.G. Lisle, R.D. Braddock, and R.E. Smith. 1982. The
three-parameter infiltration equation. Soil Sci. 133:337–341.
doi:10.1097/00010694-198206000-00001

VZJ | Advancing Critical Zone Science

Parlange, J.-Y., T.S. Steenhuis, R. Haverkamp, D.A. Barry, P.J. Culligan, W.L.
Hogarth, et al. 1999. Soil properties and water movement. In: M.B.
Parlange and J.W. Hopmans, editors, Vadose zone hydrology: Cutting
across disciplines. Oxford Univ. Press, New York. p. 99–129.
Parlange, M.B., A.T. Cahill, D.R. Nielsen, J.W. Hopmans, and O. Wendroth.
1998. Review of heat and water movement in field soils. Soil Tillage
Res. 47:5–10. doi:10.1016/S0167-1987(98)00066-X
Peck, A.J. 1983. Field variability of soil physical properties. In: D. Hillel, editor,
Advances in irrigation. Vol. 2. Academic Press, San Diego. p. 189–221.
Philip, J.R. 1957. The theory of infiltration: 1. The infiltration equation and its
solution. Soil Sci. 83:345–358. doi:10.1097/00010694-195705000-00002
Philip, J.R. 1967. Sorption and infiltration in heterogeneous media. Aust. J.
Soil Res. 5:1–10. doi:10.1071/SR9670001
Philip, J.R. 1969. Theory of infiltration. Adv. Hydrosci. 5:215–296.
Philip, J.R. 1970. Flow in porous media. Annu. Rev. Fluid Mech. 2:177–204.
doi:10.1146/annurev.fl.02.010170.001141
Philip, J.R. 1975. Stability analysis of infiltration. Soil Sci. Soc. Am. Proc.
39:1042–1049. doi:10.2136/sssaj1975.03615995003900060013x
Philip, J.R., and D.A. de Vries. 1957. Moisture movement in porous
materials under temperature gradients. Eos Trans. AGU 38:222–232.
doi:10.1029/TR038i002p00222
Philip, J.R., and J.H. Knight. 1974. On solving the unsaturated flow
equation: 3. New quasi-analytical technique. Soil Sci. 117:1–13.
doi:10.1097/00010694-197401000-00001
Poesen, J., J. Nachtergaele, G. Verstraeten, and C. Valentin. 2003. Gully
erosion and environmental change: Importance and research needs.
Catena 50:91–133. doi:10.1016/S0341-8162(02)00143-1
Poulovassilis, A. 1962. Hysteresis of pore water: An applicationof the concept of independent domains. Soil Sci. 93:405–412.
doi:10.1097/00010694-196206000-00007
Prasad, S.N., and M.J.M. Römkens. 1982. An approximate integral solutionof vertical infiltration under changing boundary conditions. Water
Resour. Res. 18:1022–1028. doi:10.1029/WR018i004p01022
Pruess, K. 1991. TOUGH2: A general-purpose numerical simulator for
multiphase fluid and heat flow. Lawrence Berkeley Lab., Berkeley, CA.
doi:10.2172/5212064
Puigdefábregas, J. 2005. The role of vegetation patterns in structuring
runoff and sediment fluxes in drylands. Earth Surf. Processes Landforms 30:133–147. doi:10.1002/esp.1181
Qi, J., P.A. Vermeer, and G. Cheng. 2006. A review of the influence of freeze–
thaw cycles on soil geotechnical properties. Permafrost Periglacial
Processes 17:245–252. doi:10.1002/ppp.559
Qu, W., H.R. Bogena, J.A. Huisman, J. Vanderborght, M. Schuh, E. Priesack,and H. Vereecken. 2015. Predicting subgrid variability of soil watercontent from basic soil information. Geophys. Res. Lett. 42:789–796.
doi:10.1002/2014GL062496
Quirk, J.P. 1994. Interparticle forces: A basis for the interpretation of soil physical behaviour. Adv. Agron. 53:121–183.
doi:10.1016/S0065-2113(08)60614-8
Quirk, J.P., and R.K. Schofield. 1955. The effect of electrolyte
concentration of soil permeability. J. Soil Sci. 6:163–178.
doi:10.1111/j.1365-2389.1955.tb00841.x
Raats, P.A.C. 1973. Unstable wetting fronts in uniform
and nonuniform soils. Soil Sci. Soc. Am. J. 37:681–685.
doi:10.2136/sssaj1973.03615995003700050017x
Raats, P.A.C. 2001. Developments in soil-water physics since the mid 1960s.
Geoderma 100:355–387. doi:10.1016/S0016-7061(01)00028-3
Raats, P.A.C., and A. Klute. 1969. One-dimensional simultaneous
motion of the aqueous phase and the solid phase of saturated and partially saturated porous media. Soil Sci. 107:329–333.
doi:10.1097/00010694-196905000-00003
Raats, P.A.C., D. Smiles, and A.W. Warrick, editors. 2002. Environmental mechanics: Water, mass and energy transfer in the biosphere. Geophys,
Monogr. Ser. 129. Am. Geophys. Union, Washington, DC.
Rahman, M., and R. Rosolem. 2017. Towards a simple representation ofchalk hydrology in land surface modelling. Hydrol. Earth Syst. Sci.
21:459–471. doi:10.5194/hess-21-459-2017
Rahmati, M., L. Weihermüller, J. Vanderborght, Y.A. Pachepsky, L. Mao,
S.H. Sadeghi, et al. 2018. Development and analysis of the Soil Water

p. 50 of 53

Infiltration Global database. Earth Syst. Sci. Data 10:1237–1263.
doi:10.5194/essd-10-1237-2018
Regalado, C.M., and A. Ritter. 2008. Cokriging spatial interpolationof soil water dependent repellency parameters determinedwith two different tests. Soil Sci. Soc. Am. J. 72:1683–1693.
doi:10.2136/sssaj2008.0008
Rengasamy, P., and K. Olsson. 1991. Sodicity and soil structure. Aust. J. Soil
Res. 29:935–952. doi:10.1071/SR9910935
Reynolds, C.A., T.J. Jackson, and W.J. Rawls. 2000. Estimating soil waterholding capacities by linking the Food and Agriculture Organizationsoil map of the world with global pedon databases and continuous pedotransfer functions. Water Resour. Res. 36:3653–3662.
doi:10.1029/2000WR900130
Reynolds, W.D., and D.E. Eldrick. 1985. In situ measurements offield-saturated hydraulic conductivity, sorptivity and a-parameter using the Guelph permeameter. Soil Sci. 140:292–302.
doi:10.1097/00010694-198510000-00008
Richards, L.A. 1931. Capillary conduction of liquids through porous mediums. Physics 1:318–333. doi:10.1063/1.1745010
Ritsema, C.J., and L.W. Dekker. 1998. Three-dimensional patterns ofmoisture, water repellency, bromide and pH in a sandy soil. J. Contam.
Hydrol. 31:295–313. doi:10.1016/S0169-7722(97)00067-3
Roeckner, E., G. Bäuml, L. Bonaventura, R. Brokopf, M. Esch, M. Giorgetta,
et al. 2003. The atmospheric general circulation model ECHAM 5: I.
Model description. MPI Rep. 349. Max Planck Inst., Hamburg.
Römkens, M.J.M., R.L. Baumhardt, J.-Y. Parlange, F.D. Whisler, M.B. Parlange,
and S.N. Prasad. 1986a. Effect of rainfall characteristics on seal hydraulicconductance. In: F. Callebaut et al., editors, Assessment of soil surface
sealing and crusting. Ghent Univ. Press, Ghent, Belgium. p. 228–235.
Römkens, M.J.M., R.L. Baumhardt, M. Parlange, F.D. Whisler, J.-Y. Parlange,
and S.N. Prasad. 1986b. Rain-induced surface seals: Their effect on
ponding and infiltration. Ann. Geophys. 4:417–424.
Römkens, M.J.M., and S.N. Prasad. 1992. A spectral series approach toinfiltration for crusted and non-crusted soils. In: M.E. Sumner and B.A.
Stewart, editors, Soil crusting: Chemical and physical processes. Lewis
Publ., Boca Raton, FL. p. 151–178.
Ross, P.J., R. Haverkamp, and J.-Y. Parlange. 1996. Calculating parameters
for infiltration equations from soil hydraulic functions. Transp. PorousMedia 24:315–339. doi:10.1007/BF00154096
Rubin, J., and R. Steinhardt. 1964. Soil water relations during rain
infiltration: III. Water uptake at incipient ponding. Soil Sci. Soc. Am. J.
28:614–619. doi:10.2136/sssaj1964.03615995002800050013x
Russo, D. 2005. Salinity/Physical effects. In: D. Hillel, editor, Encyclopediaof soils in the environment 3. Elsevier Ltd, Oxford, UK. p. 442–453.
doi:10.1016/B0-12-348530-4/00529-4
Russo, D., I. Russo, and A. Laufer. 1997. On the spatial variability of parameters of the unsaturated hydraulic conductivity. Water Resour. Res.
33:947–956. doi:10.1029/96WR03947
Sadeghi, M., B. Ghahraman, A.N. Ziaei, K. Davary, and K. Reichardt. 2012.
Additional scaled solutions to Richards’ equation for infiltration anddrainage. Soil Tillage Res. 119:60–69. doi:10.1016/j.still.2011.12.004
Sahoo, A.K., P.A. Dirmeyer, P.R. Houser, and M. Kafatos. 2008. A study ofland surface processes using land surface models over the Little
River Experimental Watershed, Georgia. J. Geophys. Res. 113:D20121.
doi:10.1029/2007JD009671
Samaniego, L., R. Kumar, and S. Attinger. 2010. Multiscale parameterregionalization of a grid-based hydrologic model at the mesoscale.
Water Resour. Res. 46:W05523. doi:10.1029/2008WR007327
Samaniego, L., R. Kumar, S. Thober, O. Rakovec, M. Zink, N. Wanders, et al.
2017. Toward seamless hydrologic predictions across spatial scales.
Hydrol. Earth Syst. Sci. 21:4323–4346. doi:10.5194/hess-21-4323-2017
Sato, N., P.J. Sellers, D.A. Randall, E.K. Schneider, J. Shukla, J.L. KinterIII, et al. 1989. Effects of implementing the simple biospheremodel in a general circulation model. J. Atmos. Sci. 46: 2757–2782.
doi:10.1175/1520-0469(1989)046<2757:Eoitsb>2.0.Co;2
Saulnier, G.M., and R. Datin. 2004. Analytical solution to a bias in the
TOPMODEL framework balance. Hydrol. Processes 18:1195–1218.
doi:10.1002/hyp.1346
Scanlon, B. 2004. Review of HYDRUS-1D. Southw. Hydrol. 3(4):37.

VZJ | Advancing Critical Zone ScienceSchaake, J.C., V.I. Koren, Q.Y. Duan, K. Mitchell, and F. Chen. 1996.
Simple water balance model for estimating runoff at different
spatial and temporal scales. J. Geophys. Res. 101(D3):7461–7475.
doi:10.1029/95JD02892
Schaap, M.G., and F.J. Leij. 1998. Database-related accuracy anduncertainty of pedotransfer functions. Soil Sci. 163:765–779.
doi:10.1097/00010694-199810000-00001
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA:
A computer program for estimating soil hydraulic parameterswith hierarchical pedotransfer functions. J. Hydrol. 251:163–176.
doi:10.1016/S0022-1694(01)00466-8
Schlüter, S., J. Vanderborght, and H.J. Vogel. 2012. Hydraulic non-equilibrium during infiltration induced by structural connectivity. Adv. Water
Resour. 44:101–112. doi:10.1016/j.advwatres.2012.05.002
Scott, P.S., G.J. Farquhar, and N. Kouwen. 1983. Hysteretic effects on net
infiltration. In: Advances in infiltration: Proceedings of the National
Conference on Advances in Infiltration, Chicago. 12–13 Dec. 1983.
ASAE Publ. 83-11. Am. Soc. Agric. Eng., St. Joseph, MI. p. 163–170.
Selker, J.S., and S. Assouline. 2017. An explicit, parsimonious, andaccurate estimate for ponded infiltration into soils using theGreen and Ampt approach. Water Resour. Res. 53:7481–7487.
doi:10.1002/2017WR021020
Selker, J.S., J. Duan, and J.-Y. Parlange. 1999. Green and Ampt infiltration into soils of variable pore size with depth. Water Resour. Res.
35:1685–1688. doi:10.1029/1999WR900008
Sellers, P.J., Y. Mintz, Y.C. Sud, and A. Dalcher. 1986. A simple biospheremodel (Sib) for use within general-circulation models. J. Atmos. Sci.
43:505–531. doi:10.1175/1520-0469(1986)043<0505:Asbmfu>2.0.Co;2.
Sellers, P.J., D.A. Randall, G.J. Collatz, J.A. Berry, C.B. Field, D.A. Dazlich,
et al. 1996. A revised land surface parameterization (SiB2) for
atmospheric GCMS: I. Model formulation. J. Climate 9: 676–705.
doi:10.1175/1520-0442(1996)009<0676:Arlspf>2.0.Co;2.
Seneviratne, S.I., T. Corti, E.L. Davin, M. Hirschi, E.B. Jaeger, I. Lehner,
et al. 2010. Investigating soil moisture–climate interactions
in a changing climate: A review. Earth Sci. Rev. 99:125–161.
doi:10.1016/j.earscirev.2010.02.004
Sharma, M.L., G.A. Gander, and C.G. Hunt. 1980. Spatial variability of infiltration in a watershed. J. Hydrol. 45:101–122.
doi:10.1016/0022-1694(80)90008-6
Sharma, M.L., and E. Seely. 1979. Spatial variability and its effects on areal
infiltration. In: Proceedings Hydrology and Water Resources Symposium,
Perth, WA. 10–12 Sept. 1979. Inst. Eng., Perth, WA, Australia. p. 69–73.
Sharma, S.K., B.P. Mohanty, and J. Zhu. 2006. Including topography andvegetation attributes for developing pedotransfer functions. Soil Sci.
Soc. Am. J. 70:1430–1440. doi:10.2136/sssaj2005.0087
She, H.Y., and B.E. Sleep. 1998. The effect of temperature on capillary pressure–saturation relationships for air–water and perchloroethylene–water systems. Water Resour. Res. 34:2587–2597. doi:10.1029/98WR01199
Shellito, P.J., E.E. Small, and M.H. Cosh. 2016. Calibration of Noah soil
hydraulic property parameters using surface soil moisture from SMOS
and basinwide in situ observations. J. Hydrometeorol. 17:2275–2292.
doi:10.1175/JHM-D-15-0153.1
Shi, Y.N., K.J. Davis, F.Q. Zhang, and C.J. Duffy. 2014. Evaluation of
the parameter sensitivities of a coupled land surface hydrologic
model at a critical zone observatory. J. Hydrometeorol. 15:279–299.
doi:10.1175/JHM-D-12-0177.1
Shouse, P.J., and B.P. Mohanty. 1998. Scaling of near-saturated hydraulicconductivity measured using disc infiltrometers. Water Resour. Res.
34:1195–1205. doi:10.1029/98WR00318
Silva, V.R., D.J. Reinert, and J.M. Reichert. 2000. Soil density, chemical
attributes and maize root distribution as affected by grazing and soil
management. (In Portuguese, with English abstract.) Rev. Bras. Cienc.
Solo 24:191–199. doi:10.1590/S0100-06832000000100021
Šimůnek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003.
Review and comparison of models for describing non-equilibrium
and preferential flow and transport in the vadose zone. J. Hydrol.
272:14–35. doi:10.1016/S0022-1694(02)00252-4
Šimůnek, J., and M.Th. van Genuchten. 2008. Modeling nonequilibrium
flow and transport processes using HYDRUS. Vadose Zone J. 7:782–
797. doi:10.2136/vzj2007.0074

p. 51 of 53

Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2008. Development andapplications of the HYDRUS and STANMOD software packages andrelated codes. Vadose Zone J. 7:587–600. doi:10.2136/vzj2007.0077

Topp, G.C. 1971. Soil-water hysteresis: The domain theory extended
to pore interaction conditions. Soil Sci. Soc. Am. J. 35:219–225.
doi:10.2136/sssaj1971.03615995003500020017x

Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2016. Recent developments and applications of the HYDRUS computer software packages.
Vadose Zone J. 15(7). doi:10.2136/vzj2016.04.0033

Tyler, S.W., and S.W. Wheatcraft. 1990. Fractal processes in soil water retention. Water Resour. Res. 26:1047–1054. doi:10.1029/WR026i005p01047

Sisson, J.B., and P.J. Wierenga. 1981. Spatial variability of steady-state infiltration rates as a stochastic process. Soil Sci. Soc. Am. J. 45:699–704.
doi:10.2136/sssaj1981.03615995004500040005x
Sivapalan, M., and E.F. Wood. 1986. Spatial heterogeneity and scale in
the infiltration response of catchments. In: V.K. Gupta, et al., editors,
Scale problems in hydrology. D. Reidel, Hingham, MA. p. 81–106.
doi:10.1007/978-94-009-4678-1_5
Sivapalan, M., and R.A. Woods. 1995. Evaluation of the effects of generalcirculation models subgrid variability and patchiness of rainfall andsoil-moisture on land-surface water-balance fluxes. Hydrol. Processes
9:697–717. doi:10.1002/hyp.3360090515
Skaggs, R.W. 1982. Infiltration. In: C.T. Haan et al., editors, Hydrologic
modeling of small watersheds. ASAE Monogr. 5. Am. Soc. Agric. Eng.,
St. Joseph, MI. p. 121–166.
Smiles, D.E. 1974. Infiltration into a swelling material. Soil Sci. 117:140–147.
doi:10.1097/00010694-197403000-00002
Smith, R.E. 1972. The infiltration envelope: Results from a theoretical infiltrometer. J. Hydrol. 17:1–22. doi:10.1016/0022-1694(72)90063-7
Smith, R.E. 1990. Analysis of infiltration through a twolayer soil profile. Soil Sci. Soc. Am. J. 54:1219–1227.
doi:10.2136/sssaj1990.03615995005400050004x
Smith, R.E., C. Corradini, and F. Melone. 1999. A conceptual model for
infiltration and redistribution in crusted soils. Water Resour. Res.
35:1385–1393. doi:10.1029/1998WR900046
Smith, R.E., and D.C. Goodrich. 2000. Model For rainfall excess patterns on randomly heterogeneous areas. J. Hydrol. Eng. 5:355–362.
doi:10.1061/(ASCE)1084-0699(2000)5:4(355)Smith, R.E., and R.H.B. Hebbert. 1979. A Monte Carlo analysis of thehydrologic effects of spatial variability of infiltration. Water Resour. Res.
15:419–429. doi:10.1029/WR015i002p00419
Smith, R.E., K.R.J. Smettem, P. Broadbridge, and D.A. Woolhiser. 2002. Infiltration theory for hydrologic applications. Water Resour. Monogr. 15.
Am. Geophys. Union, Washington, DC. doi:10.1029/WM015
Soet, M., R.J. Ronda, J.N.M. Stricker, and A.J. Dolman. 2000. Land surface
scheme conceptualisation and parameter values for three sites
with contrasting soils and climate. Hydrol. Earth Syst. Sci. 4:283–294.
doi:10.5194/hess-4-283-2000
Sposito, G. 1975. Steady vertical flows in swelling soils. Water Resour. Res.
11:461–464. doi:10.1029/WR011i003p00461
Stange, C., and R. Horn. 2005. Modeling the soil water retention curvefor conditions of variable porosity. Vadose Zone J. 4:602–613.
doi:10.2136/vzj2004.0150
Strudley, M.W., T.R. Green, and J.C. Ascough. 2008. Tillage effects on soil
hydraulic properties in space and time: State of the science. Soil Tillage Res. 99:4–48. doi:10.1016/j.still.2008.01.007
Sun, S.F., and Y.K. Xue. 2001. Implementing a new snow scheme insimplified simple biosphere model. Adv. Atmos. Sci. 18:335–354.
doi:10.1007/BF02919314
Swartzendruber, D. 1987. Rigorous derivation and interpretation of theGreen and Ampt equation. In: Y.S. Fok, editor, Infiltration development
and application. Water Resour. Res. Ctr., Univ. of Hawaii, Honolulu.
Talsma, T., and J.-Y. Parlange. 1972. One dimensional vertical infiltration.
Aust. J. Soil Res. 10:143–150. doi:10.1071/SR9720143
Täumer, K., H. Stoffregen, and G. Wessolek. 2005. Determination of
repellency distribution using soil organic matter and water content.
Geoderma 125:107–115. doi:10.1016/j.geoderma.2004.07.004
Thomas, P.J., J.C. Baker, L.W. Zelazny, and D.R. Hatch. 2000. Relationshipof map unit variability to shrink–swell indicators. Soil Sci. Soc. Am. J.
64:262–268. doi:10.2136/sssaj2000.641262x
Tillotson, P.M., and D.R. Nielsen. 1984. Scale factors in soil science. Soil Sci.
Soc. Am. J. 48:953–959. doi:10.2136/sssaj1984.03615995004800050001x

VZJ | Advancing Critical Zone Science

Unger, P.W. 1991. Overwinter changes in physical properties of no-tillage soil. Soil Sci. Soc. Am. J. 55:778–782.
doi:10.2136/sssaj1991.03615995005500030024x
van Dam, J.C., P. Groenendijk, R.F.A. Hendriks, and J.G. Kroes. 2008. Advances of modeling water flow in variably saturated soils with SWAP.
Vadose Zone J. 7:640–653. doi:10.2136/vzj2007.0060
van Dam, J.C., J.M.H. Hendrickx, H.C. Van Ommen, M.H. Bannink, M.Th.
van Genuchten, and L.W. Dekker. 1990. Water and solute movement
in a coarse-textured water-repellent field soil. J. Hydrol. 120:359–379.
doi:10.1016/0022-1694(90)90159-U
van den Hurk, B., and P. Viterbo. 2003. The Torne-Kalix PILPS 2(e) experiment as a test bed for modifications to the ECMWF land surface scheme.
Global Planet. Change 38:165–173. doi:10.1016/S0921-8181(03)00027-4
Vandervaere, J.-P., M. Vauclin, R. Haverkamp, C. Peugeot, J.-L.
Thony, and M. Gilfedder. 1998. Prediction of crust-induced
surface runoff with disc infiltrometer data. Soil Sci. 163:9–21.
doi:10.1097/00010694-199801000-00003
van Genuchten, M.Th. 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J.
44:892–898. doi:10.2136/sssaj1980.03615995004400050002x
Van Looy, K., J. Bouma, M. Herbst, J. Koestel, B. Minasny, U. Mishra,
et al. 2017. Pedotransfer functions in Earth system science:
Challenges and perspectives. Rev. Geophys. 55:1199–1256.
doi:10.1002/2017RG000581
Vaught, R., K.R. Brye, and D.M. Miller. 2006. Relationships among coefficient of linear extensibility and clay fractions in expansive, stoney soils.
Soil Sci. Soc. Am. J. 70:1983–1990. doi:10.2136/sssaj2006.0054
Vereecken, H., R. Kasteel, J. Vanderborght, and T. Harter. 2007. Upscaling
hydraulic properties and soil water flow processes in heterogeneous
soils: A review. Vadose Zone J. 6:1–28. doi:10.2136/vzj2006.0055
Vereecken, H., J. Maes, and J. Feyen. 1990. Estimating unsaturated hydraulic conductivity from easily measured soil properties. Soil Sci. 149:1–12.
doi:10.1097/00010694-199001000-00001
Vereecken, H., J. Maes, J. Feyen, and P. Darius. 1989. Estimating the soil-moisture retention characteristics from texture, bulk-density, and carbon content. Soil Sci. 148:389–403.
doi:10.1097/00010694-198912000-00001
Vereecken, H., M. Weynants, M. Javaux, Y. Pachepsky, M.G. Schaap, andM.Th. van Genuchten. 2010. Using pedotransfer functions to estimate
the van Genuchten–Mualem soil hydraulic properties: A review. Vadose Zone J. 9:795–820. doi:10.2136/vzj2010.0045
Verhoef, A., and G. Egea-Cegarra. 2013. Soil water and its management. In: P.J. Gregory and S. Nortcliff, editors, Soil conditions
and plant growth. Blackwell Publ., Chichester, UK. p. 269–322.
doi:10.1002/9781118337295.ch9
Vieira, S.R., D.R. Nielsen, and J.W. Biggar. 1981. Spatial variability offield-measured infiltration rate. Soil Sci. Soc. Am. J. 45:1040–1048.
doi:10.2136/sssaj1981.03615995004500060007x
Viterbo, P. and A.C.M. Beljaars. 1995. An improved land-surface parameterization scheme in the ECMWF model and its validation. J. Climate
8:2716–2748. doi:10.1175/1520-0442(1995)008<2716:ailsps>2.0.co;2
Vogel, H.-J., I. Cousin, O. Ippisch, and P. Bastian. 2006. The dominant roleof structure for solute transport in soil: Experimental evidence andmodelling of structure and transport in a field experiment. Hydrol.
Earth Syst. Sci. 10:495–506. doi:10.5194/hess-10-495-2006
Vogel, H.-J., and O. Ippisch. 2008. Estimation of a critical spatial discretization limit for solving Richards’ equation at large scales. Vadose Zone J.
7:112–114. doi:10.2136/vzj2006.0182
Vogel, H.J., and K. Roth. 1998. A new approach for determining effective soil hydraulic functions. Eur. J. Soil Sci. 49:547–556.
doi:10.1046/j.1365-2389.1998.4940547.x
Vrugt, J.A., C.J.F. ter Braak, C.G.H. Diks, B.A. Robinson, J.M. Hyman
and D. Higdon. 2009. Accelerating Markov chain Monte Carlosimulation by differential evolution with self-adaptive randomized

p. 52 of 53

subspace sampling. Int. J. Nonlinear Sci. Numer. Simul. 10:273–290.
doi:10.1515/IJNSNS.2009.10.3.273
Walko, R.L., and R. Avissar. 2008. The Ocean–Land–Atmosphere Model(OLAM): I. Shallow-water tests. Mon. Weather Rev. 136:4033–4044.
doi:10.1175/2008MWR2522.1
Walko, R.L., L.E. Band, J. Baron, T.G.F. Kittel, R. Lammers, T.J. Lee, etal. 2000. Coupled atmosphere–biophysics–hydrology models
for environmental modeling. J. Appl. Meteorol. 39:931–944.
doi:10.1175/1520-0450(2000)039<0931:CABHMF>2.0.CO;2
Wang, L.L., D.H. Chen, and H.J. Bao. 2016. The improved Noah land surface
model based on storage capacity curve and Muskingum method
and application in GRAPES model. Atmos. Sci. Lett. 17:190–198.
doi:10.1002/asl.642
Wang, Y.P., E. Kowalczyk, R. Leuning, G. Abramowitz, M.R. Raupach, B.
Pak, et al. 2011. Diagnosing errors in a land surface model (CABLE)in the time and frequency domains. J. Geophys. Res. 116:G01034.
doi:10.1029/2010jg001385
Wang, Z., J. Feyen, M.Th. van Genuchten, and D.R. Nielsen. 1998. Air entrapment effects on infiltration rate and flow instability. Water Resour.
Res. 34:213–222. doi:10.1029/97WR02804
Wang, Z., A. Tuli, and W.A. Jury. 2003a. Unstable flow during redistribution inhomogeneous soil. Vadose Zone J. 2:52–60. doi:10.2136/vzj2003.5200
Wang, Z., L.S. Wu, T. Harter, J.H. Lu, and W.A. Jury. 2003b. A field studyof unstable preferential flow during soil water redistribution. Water
Resour. Res. 39:1075. doi:10.1029/2001WR000903
Warrick, A.W. 2003. Soil water dynamics. Oxford Univ. Press, New York.
Warrick, A.W., and A.A. Hussen. 1993. Scaling of Richards’ equation for infiltration and drainage. Soil Sci. Soc. Am. J. 57:15–18.
doi:10.2136/sssaj1993.03615995005700010004x
Warrick, A.W., and D.R. Nielsen. 1980. Spatial variability of soil physicalproperties in the field. In: D. Hillel, editor, Applications of soil physics.
Academic Press, San Diego. p. 319–344.
Warrick, A.W., and T.C.J. Yeh. 1990. One-dimensional, steady vertical flow in a layered soil profile. Adv. Water Resour. 13:207–210.
doi:10.1016/0309-1708(90)90042-3
Warrilow, D.A., A.B. Sangster, and A. Slingo. 1986. Modelling of land surface processes and their influence on European climate. Dyn. Climatol.
DCTN38. Meteorol. Office, Bracknell, UK.
Weynants, M., H. Vereecken, and M. Javaux. 2009. Revisiting Vereecken
pedotransfer functions: Introducing a closed-form hydraulic model.
Vadose Zone J. 8:86–95. doi:10.2136/vzj2008.0062
Whisler, F.D., A.A. Curtis, A. Niknam, and M.J.M. Römkens. 1979. Modelinginfiltration as affected by soil crusting. In: H.J. Morel Seytoux, editor,
Surface and subsurface hydrology: Proceedings of the Fort Collins
Third International Hydrology Symposium on Theoretical and Applied
Hydrology, Fort Collins, CO. 27–29 July 1977. Water Resour. Publ., Fort
Collins, CO. p. 400–413.
White, I., and M.J. Sully. 1992. On the variability and use of the hydraulicconductivity alpha parameter in stochastic treatments of unsaturated
flow. Water Resour. Res. 28:209–213. doi:10.1029/91WR02198
Williams, P.J., and T.P. Burt. 1974. Measurement of hydraulic conductivityof frozen soils. Can. Geotech. J. 11:647–650. doi:10.1139/t74-066
Witter, J.V., P.D. Jungerius, and M.J. ten Harkel. 1991. Modelling watererosion and the impact of water repellency. Catena 18:115–124.
doi:10.1016/0341-8162(91)90011-L
Wood, E.F., D.P. Lettenmaier, X. Liang, D. Lohmann, A. Boone, S. Chang, etal. 1998. The Project for Intercomparison of Land-surface Parameterization Schemes (PILPS) Phase 2(c) Red–Arkansas River basin experiment:
1. Experiment description and summary intercomparisons. Global
Planet. Change 19:115–135. doi:10.1016/S0921-8181(98)00044-7
Woolhiser, D.A., R.E. Smith, and J.-V. Giraldez. 1996. Effects of spatial variability of saturated hydraulic conductivity on Hortonian overland flow.
Water Resour. Res. 32:671–678. doi:10.1029/95WR03108
Wösten, J.H.M., A. Lilly, A. Nemes, and C. Le Bas. 1999. Development anduse of a database of hydraulic properties of European soils. Geoderma
90:169–185. doi:10.1016/S0016-7061(98)00132-3
Wösten, J.H.M., Y.A. Pachepsky, and W.J. Rawls. 2001. Pedotransfer
functions: Bridging the gap between available basic soil data

VZJ | Advancing Critical Zone Scienceand missing soil hydraulic characteristics. J. Hydrol. 251:123–150.
doi:10.1016/S0022-1694(01)00464-4
Xue, Y., P.J. Sellers, J.L. Kinter and J. Shukla. 1991. A simplified biosphere model for global climate studies. J. Climate 4:345–364.
doi:10.1175/1520-0442(1991)004<0345:Asbmfg>2.0.Co;2Xue, Y.K., F.J. Zeng, and C.A. Schlosser. 1996. SSiB and its sensitivity to soil
properties: A case study using HAPEX-Mobilhy data. Global Planet.
Change 13:183–194. doi:10.1016/0921-8181(95)00045-3
Yager, R.R. 2001. The power average operator. IEEE Trans. Syst. Man Cybern.
A 31:724–731. doi:10.1109/3468.983429
Yang, K., T. Koike, B.S. Ye, and L. Bastidas. 2005. Inverse analysis of the roleof soil vertical heterogeneity in controlling surface soil state and energy partition. J. Geophys. Res. 110:D08101. doi:10.1029/2004jd005500
Yang, Z.-L., and R.E. Dickinson. 1996. Description of the Biosphere–Atmosphere Transfer Scheme (BATS) for the Soil Moisture Workshop andevaluation of its performance. Global Planet. Change 13:117–134.
doi:10.1016/0921-8181(95)00041-0
Yang, Z.-L., and G.-Y. Niu. 2003. The versatile integrator of surface andatmosphere processes: 1. Model description. Global Planet. Change
38:175–189. doi:10.1016/S0921-8181(03)00028-6
Young, I.M., and J.W. Crawford. 2004. Interactions and self-organization in the soil–microbe complex. Science 304:1634–1637.
doi:10.1126/science.1097394
Zaslavsky, D. 1964. Theory of unsaturated flow into a non-uniform soil
profile. Soil Sci. 97:400–410. doi:10.1097/00010694-196406000-00006
Zhan, X.W., Y.K. Xue, and G.J. Collatz. 2003. An analytical approach for estimating CO2 and heat fluxes over the Amazonian region. Ecol. Modell.
162:97–117. doi:10.1016/S0304-3800(02)00405-2
Zhang, Y.Q., H.X. Zheng, F.H.S. Chiew, J. Pena-Arancibia, and X.Y. Zhou.
2016. Evaluating regional and global hydrological models against
streamflow and evapotranspiration measurements. J. Hydrometeorol.
17:995–1010. doi:10.1175/JHM-D-15-0107.1
Zhao, W., and A. Li. 2015. A review on land surface processesmodelling over complex terrain. Adv. Meteorol. 17:607181.
doi:10.1155/2015/607181
Zhou, X.Y., Y.Q. Zhang, Y.P. Wang, H.Q. Zhang, J. Vaze, L. Zhang, et al. 2012.
Benchmarking global land surface models against the observed meanannual runoff from 150 large basins. J. Hydrol. 470–471:269–279.
doi:10.1016/j.jhydrol.2012.09.002
Zhu, J., and B.P. Mohanty. 2002a. Spatial averaging of van Genuchten
hydraulic parameters for steady-state flow in heterogeneous soils: A
numerical study. Vadose Zone J. 1:261–272. doi:10.2136/vzj2002.2610
Zhu, J., and B.P. Mohanty. 2002b. Upscaling of soil hydraulic properties forsteady state evaporation and infiltration. Water Resour Res 38(9):1178.
doi:10.1029/2001WR000704.
Zhu, J., and B.P. Mohanty. 2003. Effective hydraulic parameters for steadystate vertical flow in heterogeneous soils. Water Resour. Res. 39:1227.
doi:10.1029/2002WR001831
Zhu, J., and B.P. Mohanty. 2004. Soil hydraulic parameter upscaling forsteady-state flow with root water uptake. Vadose Zone J. 3:1464–1470.
doi:10.2136/vzj2004.1464
Zhu, J., and B.P. Mohanty. 2006. Effective scaling factor for transient infiltration in heterogeneous soils. J. Hydrol. 319:96–108.
doi:10.1016/j.jhydrol.2005.07.004
Zhu, J., B.P. Mohanty, and N.N. Das. 2006. On the effective averaging
schemes of hydraulic properties at the landscape scale. Vadose Zone J.
5:308–318. doi:10.2136/vzj2005.0035
Zhu, J., B.P. Mohanty, A.W. Warrick, and M.Th. van Genuchten. 2004.
Correspondence and upscaling of hydraulic functions for steadystate flow in heterogeneous soils. Vadose Zone J. 3:527–533.
doi:10.2136/vzj2004.0527
Zobler, L. 1986. A world soil file for global climate modeling. Tech. Mem.
87802. NASA Goddard Space Flight Ctr., Greenbelt, MD.
Zobler, L. 1999. Global soil types, 1-degree grid (Zobler). [Dataset.]
Oak Ridge Natl. Lab. Distributed Active Archive Ctr., Oak Ridge, TN.
doi:10.3334/ORNLDAAC/418

p. 53 of 53