Journal of Hydrology 565 (2018) 570–580

Contents lists available at ScienceDirect

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

Research papers

Unsaturated hydraulic properties of heterogeneously packed sands: A porescale computational study

T⁎

Mojtaba G. Mahmoodlua,b, , Martinus Th. van Genuchtenb,c, Thomas Sweijenb, Amir Raoofb
aDepartment of Watershed and Rangeland Management, Gonbad Kavous University, Gonbad Kavous, Iran
Department of Earth Sciences, Utrecht University, Utrecht, The NetherlandscCenter for Environmental Studies, CEA, São Paulo University, Rio Claro, SP, Brazilb

A R T I C LE I N FO

A B S T R A C T

This manuscript was handled by C. Corradini,
Editor-in-Chief, with the assistance of Corrado
Corradini, Associate Editor

The macroscopic eﬀects of soil textural heterogeneity and ﬁne-scale soil layering on unsaturated ﬂow remain
poorly understood. In this study we used the Discrete Element Method (DEM) to numerically generate artiﬁcial
particle packings of heterogeneous sands from which the pore structure can be extracted. Packings were generated from known grain size distributions following sand grain mixing and ﬁne-scale proﬁle layering. Five sandswith diﬀerent mean grain sizes were used to investigate the eﬀects of heterogeneity on their pore structure,including changes in the permeability and soil water retention (SWR). Heterogeneous media were created bymixing a relatively ﬁne sand with coarser sands. In addition, we created layered media having sharp as well astransitional interfaces between the two sand samples. Mixing ﬁne and coarse sand caused reductions in the
average pore body and pore throat sizes of the coarser sands, and hence also in their intrinsic permeability. The
layered media with a transitional interface showed lower porosities at the interface because of penetration ofsmall particles into the larger pores, while the porosities increased at sharp interfaces. The nonlinear relationships between permeability and the average pore body and throat radii were explored using diﬀerent unimodal
sands. Mixing ﬁne and coarse sand caused a decrease in the capillary pressure at a given water content of the
new medium, and hence larger values of the van Genuchten α parameter. The SWRCs of the layered soils with asharp interface were best described using bimodal functions. Sharp interfaces caused a non-monotonic change inthe drainage curve due to discrepancies between the pore throat radii of the two adjacent sands. This change wasinﬂuenced by d50 of the sands used for the layered media. In contrast to the sharp interface, a monotonic(smooth) change in the drainage curve was observed for two sands having small diﬀerences in their d50 values.
The study gave considerable insight into how particle heterogeneity and layering aﬀects the hydraulic propertiesof unsaturated media.

Keywords:
Soil water retention curve (SWRC)
Intrinsic permeability
Heterogeneous sands
Discrete element method
Pore structure

1. Introduction
Reliable estimates of the unsaturated soil hydraulic properties are
required for eﬀective management and simulations of many near-surface hydrological processes such as inﬁltration, runoﬀ, groundwater
recharge, irrigation and optimal water use in agricultural operations,and contaminant transport in the vadose zone (Zhu and Sun, 2012; Di
Prima et al., 2016). Applications also include such diverse problems as
slope stability analyses (Mukhlisin et al., 2011; Antinoro et al., 2017),
ﬂow and seepage into and through earth dams (Alonso and Cardoso,
2010; Al-Ismaily et al., 2015) and inﬁltration into and through waste
disposal sites (Lopes et al., 2013).
The hydraulic properties of variably-saturated porous media include
the soil water retention curve relating capillary pressure and ﬂuid⁎

saturation, and the unsaturated soil hydraulic conductivity function. A
large number of experimental protocols have been developed over theyears to estimate either one or both of these functions using suchmethods as hanging water columns or sand boxes, pressure cells,pressure plate extractors, suction tables, soil freezing, falling head
methods, inﬁltrometer and permeameter devices, laboratory one-stepand multistep outﬂow methods, centrifugation and evaporation
methods, and other approaches (Dane and Hopmans, 2002; Nimmo
et al., 2002; Scanlon et al., 2002; Bittelli and Flury, 2009; Zhuang et al.,2017).
Notwithstanding the plethora of available methods, a major challenge for accurate delineation of ﬁeld-scale hydraulic properties at bothvery small and large spatial scales remains the overwhelming heterogeneity of the subsurface (Chen et al., 2012; Hohenbrink and Lischeid,

Corresponding author at: Department of Watershed and Rangeland Management, Gonbad Kavous University, Gonbad Kavous, Iran.
E-mail address: m.gharehmahmoodlu@uu.nl (M.G. Mahmoodlu).

https://doi.org/10.1016/j.jhydrol.2018.07.060
Received 3 January 2018; Received in revised form 22 June 2018; Accepted 24 July 2018
Available online 29 July 20180022-1694/ © 2018 Published by Elsevier B.V.

Journal of Hydrology 565 (2018) 570–580

M.G. Mahmoodlu et al.

reﬂecting grain heterogeneity. The coeﬃcients of uniformity (Cu),curvature (Cc), and sorting (So) were estimated for all sands and their
mixtures as follows (Dubois et al., 2009; Katukiza et al., 2014):

2015; Sakaki and Smits, 2015; Morbidelli et al., 2016). At relatively
small scales, heterogeneities are manifested by non-uniform particleand pore-size distributions, soil textural layering and soil structure(Mohanty and Zhu, 2007; Coppola et al., 2009; Huang et al., 2011;
Arora, 2012; Sakaki and Smits, 2015). The totality of micropores (deﬁned here as those smaller than 75 μm) and macropores (> 75 μm),
including their geometries and spatial distributions, are the primary
determinants of two-phase ﬂuid ﬂow and contaminant transport inporous media. Fluid ﬂow and capillary displacement processes arestrongly aﬀected by pore throat sizes since they provide the constrictions along the ﬂow paths (Perrier et al., 1996; Raoof et al., 2010;
Falode and Manuel, 2014; Mahmoodlu et al., 2016). Precise studies ofthe eﬀects of soil textural heterogeneity and soil layering or stratiﬁcation on the unsaturated soil hydraulic properties require detailed analyses of the geometry of sand grains and related pore structures in terms
of prevailing pore body and pore throat sizes. Quantiﬁcation of thesegeometry eﬀects oﬀers an opportunity to relate the soil water retention
properties to soil structural components (e.g., Perrier, et al., 1996).
Although the particle size distribution is relatively easy to measureusing a range of techniques, estimation of grain positions and pore
structures cannot be done in the same way (Gee and Or, 2002;
Mehlhorn et al., 2008; Mahmoodlu et al., 2016). X-ray microtomography (which uses x-ray radiation to create cross-sections of aphysical object) together with image processing (methods to convert animage into digital form and perform some operations on it) can be used
to obtain more directly details of pore structure geometries as deﬁned
by the size, arrangement, and connectivity of pores (e.g., Peth et al.,2010; Cnudde and Boone, 2013; Vaz et al., 2014). Unfortunately, X-ray
microtomography is still relatively expensive (Mehlhorn et al., 2008)
and image processing relatively time consuming.
In this study, we used an alternative approach to investigate theeﬀects of small-scale heterogeneities on the soil hydraulic properties.
Following Keller et al. (2013) and Mahmoodlu et al. (2016), we considered idealized grain geometries to obtain the pore structures of soils
having predeﬁned unimodal and bimodal particle size distributions(having a single peak and two peaks, respectively), as well as of layeredsoil samples. Particle packings were generated using the discrete element method (DEM), which allows one to simulate the movement ofgrains during compaction of the particles (Cundall and Strack, 1979).
Computations were carried out using the open source Yade-DEM software (Šmilauer et al., 2015) which also allows us to compute the hydraulic properties as well as the soil water retention curve (SWRC).
Once a packing of spheres was constructed, the pore space (consistingof pore bodies and pore throats) was extracted using regular triangulation (Chareyre et al., 2012), which leads to an assembly of grainbased tetrahedra. Using this assembly, the intrinsic permeability can be
determined (e.g. Tong et al., 2012; Chareyre et al., 2012) as well as the
SWRC following approaches by Yuan et al. (2015) and Sweijen et al.
(2016).
Speciﬁc objectives of our study were to: (1) use the discrete element
model to extract pore structures of diﬀerent heterogeneous particle
packings; (2) explore the eﬀects of vertical heterogeneity (i.e., soil
textural layering) and heterogeneous sphere packings on the pore
structures; and (3) investigate the eﬀects of these heterogeneities on thesoil water retention and the intrinsic permeability, including the van
Genuchten hydraulic parameters α and n.

Cu =

d60
d10

(1)

Cc =

2
d30
d10 × d60

(2)

So =

d75
d25

(3)

where dx denotes the grain diameter (in mm), for which x% of thegrains are ﬁner.
Results of the coeﬃcients of uniformity, curvature, and sorting are
given in Table 1. Because the coeﬃcient of uniformity was less than 4for the ﬁve main sands, the samples were considered to be uniformly
graded by containing essentially identical particle sizes (Purushothama,2007). The coeﬃcient of curvature was used as another useful measureof the shape of the particle size distribution. The sorting coeﬃcient,which also describes the distribution of grain size, was additionallyconsidered as a measure of the uniformity of the soil samples. The data
in Table 1 show that the three coeﬃcients (Cu, Cc and So) all indicate
similar narrow degrees of sorting of the sands. The degree of sorting of
the mixtures declined towards the coarser sands (i.e., S4, S5) because of
the creation of increasingly more bimodal grain size distributions.
3. Pore scale analysis
In this section we describe how the homogeneous and heterogenous
sand media were constructed and analyzed in terms of their hydraulicproperties using the discrete element method. We ﬁrst provide detailsof the packing procedures and resulting pore structures as deﬁned by
the size, arrangement and connectivity of the pores. We next investigate
the pore body radius distributions of the main sands, followed by astudy of heterogeneous sphere packings obtained by mixing ﬁne sandS1 with the other main sands (samples S2 to S5 in Table 1), as well as ofthe eﬀects of soil textural layering of the sand packings. The discreteelement method is subsequently used to evaluate the intrinsic permeability and soil water retention properties of the diﬀerent sand packings(i.e., the single, mixed and layered media).
3.1. Sand packings
To limit errors in the simulations while minimizing simulation
times, we decided to use for our calculations 2000 particles from agiven particle size distribution. In view of our interest to investigate theimplications of both grain size heterogeneity and soil layering, we
considered two diﬀerent ways of spherical particle packings. For theindividual sand samples (ﬁve sand samples, i.e., S1, S2, S3, S4 and S5)
and the mixtures of two sands (resulting in four sand samples; S1&S2,
S1&S3, S1&S4 and S1&S5), particles were inserted in a large box suchthat clouds of particles were generated (following procedures by
Mahmoodlu et al. 2016).
We assumed normal elasticity, tangential elasticity, and sliding as
the main processes taking place at the contact points of spheres. Normal
elasticity was calculated using a linear contact law (Belheine et al.,2009; Mahmoodlu et al., 2016). When two particles are in contact andpushed toward each other by the surrounding medium, they are ﬂattened at their contact point. The tendency of both particles to keep theirinitial shapes may cause an elastic force fn:

2. Sand properties
As input to the discrete element model (DEM), we used the synthetic
particle size distributions of ﬁve sands, each one having a relative
narrow particle size distribution as shown in Fig. 1. Particle sizes varied
from 76 to 1190 µm (Table 1). Sand S1, which had the ﬁnest texture(Fig. 1a, Table 1), was used to create heterogeneous soil samples by
mixing it with the other sands (Fig. 1b), thus allowing us to explore themixing eﬀects on grain size and the resulting pore size distributions

fn = kn δn(4)where kn is the stiﬀness of the particles and δn denotes the normal
displacement.
A tangential force, ft, may arise as a result of shear at the contact571

Journal of Hydrology 565 (2018) 570–580

M.G. Mahmoodlu et al.

Fig. 1. Synthetic particle size distributions of the ﬁve main sands (a) and their mixtures (b).

d50 of sand S2), a very sharp interface resulted (sand samples S1-S2, S1S3, S1-S4 and S1-S5 presented in Table 2). On the other hand, using acoarser sand at the bottom and ﬁner sand S1 at the top produced amuch smoother interface because of the intrusion by gravity of somesmaller particles from the top layer into the coarser pores of the bottom
layer (Fig. 2b, d and sand samples S2-S1, S3-S1, S4-S1 and S5-S1 presented in Table 2). To keep the length of both layers in the sandpackings with two layers (i.e., those having an interface) constant, themass fraction of Sand S1 in the sand packings made up of two layers
was variable among the diﬀerent samples (approximately 70%, 80%,90%, and 95% for the packings of S1-S2, S1-S3, S1-S4, and S1-S5, respectively). For instance, the packing of S1-S5 required in this way only
5% of S5. Hence, large fractions of S1 in all sand packings consisting of
two layers produced very minor discrepancies between the pore radii
and/or pore throat distributions of the S1-S2 and S1-S5 sand packings.
Compaction was accounted for in the pore-scale calculations by repeated repacking and rearrangement of the spherical particles to decrease porosities to the prescribed values, with concomitant changes in
the pore structures.

point. When the tangential force reaches a certain threshold value, thecontact point will break and cause sliding of the two particles. This
process is described byft ⩽ fn tan ϕ(5)where ϕ denotes the friction angle used to estimate the threshold value.
For each particle, a force balance is made by including forces due tocontacts with other particles, the boundary conditions, and the eﬀects
of gravity. The resulting net force gives the acceleration of the particle
using Newton’s second law (Cundall and Strack, 1979; Belheine et al.,2009).
The grain-scale parameters needed for this approach were takenfrom previous DEM studies. However, mechanical parameterization of
the particles could be done arbitrarily since the porosity was used as aproxy for the state-of-stress. Thus, a conﬁning stress was applied suchthat a particular porosity was obtained. Values of the particle stiﬀness,friction angle, and Poisson ratio were ﬁxed at 9.6 × 108 Pa, 30°, and0.04, respectively (Belheine et al., 2009). For our simulations we allowed the six surrounding boundaries of the domain to conﬁne anddrag the particles until a porosity value of 0.45 was reached. After thispoint, only the top boundary was allowed to move down to further
conﬁne the particle packings until a target porosity of 0.35 resulted.
For the sand packings consisting of two layers, diﬀerent interfaces
between the two sand layers were created during packing. For thesecases we ﬁrst inserted a number of particles into a large box such that acloud of particles was generated for the ﬁrst (bottom) layer, similarly as
for the previous (mixed) samples. But now ﬁve of the boundaries of thedomain were kept constant, while the top boundary was allowed tomove down to conﬁne the packings until the target porosity of 0.35 wasreached (the average porosity of the coarse sand). Once reached, we
froze the bottom layer and then started packing the top layer in thesame way as the ﬁrst layer. The resulting interface between the twolayers was very much dependent upon the grain diameter (or d50) of theﬁrst layer. If we used a relatively ﬁne sand as the ﬁrst layer (Fig. 2a, c),and a coarser sand as the top layer (i.e., d50 of sand S1 being less than

3.2. Pore structure
A regular triangulation method was used to extract the geometry of
the various pore structures (Chareyre et al., 2012). The method proceeds by dividing the pore space into tetrahedrons, with each tetrahedron serving as a pore unit. Each tetrahedron is then surrounding one
pore body and has four vertices that are located at the centers of thespherical particles (Fig. 3). The sides, or facets, of the tetrahedrons are
the locations at which two pore bodies are connected. The facet is thenarrowest transect between two connecting pore units and is referred toas the pore throat.
Although several studies suggest a distribution in pore connectivity
values due to the natural structure of porous media (e.g., Raoof et al.,2010; Vasilyev et al., 2012), our approach assumes that four porethroats are connected to each pore body. Each pore throat in this way is

Table 1
Selected properties of the sands used in this study.
SampleMinimum

Maximum

d10

d25

d30

d50

d60

d75

Cu

Cc

So

µm
S1
S2
S3
S4
S5
S1&S2
S1&S3
S1&S4
S1&S5

76
176
276
476
650
76
76
76
76

525
625
725
925
1190
625
725
925
1190

188
290
388
588
790
239
288
388
489

230
330
430
630
835
280
330
430
532.5

245
345
445
645
845
295
345
445
545

283
385
485
685
885
334
384
484
584

300
400
500
700
900
350
400
500
600

330
435
530
730
932
382.5
430
530
631

1.60
1.38
1.29
1.19
1.14
1.46
1.39
1.29
1.23

1.06
1.03
1.02
1.01
1.00
1.04
1.03
1.02
1.01

1.20
1.15
1.11
1.08
1.06
1.17
1.14
1.11
1.09

Cu: Uniformity coeﬃcient; Cc: Curvature coeﬃcient; So: Sorting coeﬃcient.
572

Journal of Hydrology 565 (2018) 570–580

M.G. Mahmoodlu et al.

Fig. 2. Sphere packings of two layers with diﬀerent mean particle diameters (d50) with d50(S1) < d50(S2) < d50(S3). a): a sharp interface with S1 at the bottom and S2
at the top, b) a transitional interface with S2 at bottom and S1 at the top, c) a sharp interface with S1 at bottom and S3 at the top, and d) a transitional interface with
S3 at the bottom and S1 at the top.

one to solve for the pressure distribution under fully water saturated
conditions. The approach assumes that one pressure value can be associated with each pore unit. The ﬂux, qij [L3T−1], of water from agiven pore unit i to pore unit j is then given by (Vasilyev et al., 2012):

then surrounded by three particles to which an inscribed circle can beﬁtted (Fig. 3). The radius of this circle, Rij, is taken as the pore throat
radius (Torskaya et al., 2014). Inside a pore unit, an inscribed sphere
exists that touches the four surrounding particles (e.g. Sweijen et al.,2016). The approach using inscribed spheres is often referred to as thepore ﬁnite volume method (Gladkikh and Bryant, 2005; Unsal et al.,2009; Mousavi and Bryant, 2012; Prodanović et al., 2015), while thecorresponding radius is denoted as Ri.

qij = kij(pi −pj )lij(6)where pi and pj [ML−1T−2] are the pressures in pore unit i, and j, respectively, lij [L] is the length between the centers of pore units i and j,
and kij [M−1L5T] denotes the hydraulic conductivity of the facet (porethroat) between pores i and j. The latter is given by:

3.3. Calculation of the permeability
Following a study by Chareyre et al. (2012), the pore ﬁnite volume
method as implemented in the DEM was used to compute the permeability of granular materials in the sphere packings. This method allows
573

Journal of Hydrology 565 (2018) 570–580

M.G. Mahmoodlu et al.

Table 2
van Genuchten soil hydraulic parameters of the various samples.
Sampleϴr (m3m−3)ϴs (m3m−3)

w1 (−)

α1 (kPa−1)

α2 (kPa−1)

n1 (−)n2 (−)

Description

S1
S2
S3
S4
S5
S1&S2
S1&S3
S1&S4
S1&S5
S1-S2
S1-S3
S1-S4
S1-S5
S2-S1
S3-S1
S4-S1
S5-S1

0.060
0.076
0.063
0.063
0.062
0.064
0.063
0.072
0.062
0.072
0.065
0.062
0.052
0.172
0.200
0.205
0.223

0.350
0.350
0.350
0.350
0.350
0.350
0.350
0.350
0.350
0.350
0.350
0.350
0.350
0.350
0.350
0.350
0.350

1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
1.00
0.50
0.50
0.50
1.00
1.00
1.00
1.00

0.141
0.189
0.231
0.302
0.389
0.162
0.198
0.280
0.359
0.168
0.240
0.326
0.428
0.171
0.197
0.232
0.272

–
–
–
–
–
–
–
–
–
–
0.141
0.143
0.144
–
–
–

15.3
14.9
16.1
19.1
29.5
20.0
18.4
21.4
18.2
10.0
18.6
11.8
22.3
68.
384.
347.
371.

–
–
–
–
–
–
–
–
–
–
11.1
11.8
9.9
–
–
–

individual
individual
individual
individual
individual
mixed (0.5–0.5)
mixed (0.5–0.5)
mixed (0.5–0.5)
mixed (0.5–0.5)
sharp interface
sharp interface
sharp interface
sharp interface
transition interface
transition interface
transition interface
transition interface

kij =Aijf Rij2

conditions:

2μ

(7)

Pc = Pair −Pw

where Aijf [L2] is the smallest transect of the pore throat, μ [ML−1] isTo estimate the capillary pressure - saturation curve for drainage,
the non-wetting air pressure was increased in steps to allow air to moveinto increasingly smaller pore throats. Air will enter a water-saturated
pore unit if the capillary pressure is larger than the entry pressure of thepore throat (Pe), which then can be obtained using the Young-Laplace
equation:

the dynamic viscosity, and Rij [L] is the hydraulic radius of the facet,
i.e.,

Rij =

Θij
ωij

(8)

in which Θij [L ] is the pore throat volume, and ωij [L ] the surfaceof the pore-throat. We refer readers to Chareyre et al. (2012) for moredetails. Under steady-state conditions, and assuming incompressible
ﬂow, we further use the continuity equation given by (Zhang et al.,2015):
3

2

Pe =

j=1

2γ
Rij

(11)

where γ is the surface tension and Rij the radius of the inscribed circle in
the pore throat (Fig. 3). Following a study by Zhang et al. (2015), thecapillary pressure during imbibition must be lower than the entry
pressure before water will invade an air-saturated pore unit. At eachpressure step, the equilibrium positions of the air-water interfaces
within the network were determined. This information was then used to
calculate the average saturation of the network at any given capillarypressure (Raoof and Hassanizadeh, 2012).
To estimate the complete SWRC, the simulation procedure was asfollows (Sweijen et al., 2016 and 2017; Yuan et al., 2015): After ﬁrst
selecting a capillary pressure value, we determined if air could invade awater-saturated pore body from a neighboring air-invaded pore body.
For invasion to occur, both pore bodies had to be connected to a corresponding reservoir. Disconnected (e.g. trapped) water-saturated pores
were used estimate residual water saturations. An algorithm was used
next to keep track of disconnected air and water regions to ensure thatno displacement between the disconnected pores could occur. However,in real porous media with hydrophilic surfaces, water may also be
stored in edges and corners as well as in in the form of adsorbed ﬁlms
on the grains.

4∑ qij = 0

(10)

(9)

Eqs. (6), (7) and (9) were solved using an implicit numerical scheme
to obtain pressure values (Raoof et al., 2010). The total ﬂux through themedium was subsequently used to compute the permeability of themedium.
3.4. Capillary pressure-saturation calculations
The capillary pressure-saturation (or soil water retention, SWRC)
curves were estimated using a module inside Yade-DEM that was implemented by Yuan et al. (2015) and Sweijen et al. (2016). We assumed
that the top boundary of the simulated domain was an air reservoir withpressure Pair, and the bottom boundary a water reservoir at a ﬁxed
pressure, Pw. The capillary pressure, Pc, was set equal to the pressure
diﬀerence between the two reservoirs, which is valid for quasi-static

Fig. 3. Schematic of (a) a spherical packing with (b) a pore body with radius Ri and (c) a pore throat with radius Rij as used in this study (Mahmoodlu et al., 2016).
574

Journal of Hydrology 565 (2018) 570–580

M.G. Mahmoodlu et al.

the sand mixtures (Fig. 4c, d) and the two-layered packings (Fig. 4e, f).
Results show relatively uniform distributions of the pore body and porethroat radii of the main sands due to fairly uniform increases in thegrain diameters (Fig. 4a, b). The plots indicate that particle mixing
aﬀected the pore- and throat-size distributions at the given porosity of
0.35 (Fig. 4c, d). Mixing Sand S1 with the coarser sand grains (i.e.,
Sands S4 and S5) created mixtures with much larger pore body andthroat radii as compared to Sand S1 as such. Linear relationships wereobtained between the average pore throat and pore body radii during
mixing, similarly as in our earlier study on compaction (Mahmoodlu
et al., 2016).
Since SWRCs are dominated by pore throat radii, at least in the wet
range, one may expect fairly uniform changes in the SWRC of each sand(S1 to S5) when the grain diameters increase, but less uniform changes
when considering the mixed media (Mahmoodlu et al., 2016). We obtained a linear relationship between the average pore throat and pore
body radii when mixing sand S1 with the other sands. To explore theeﬀect of particle diameter on pore size during mixing of the sands, theaverage pore body and pore throat radii of all sands were calculated
and then plotted versus d50. The overall linear correlation indicates thatthe average pore body and throat radii can be estimated reasonably
well at a given porosity value. These correlations can be useful togenerate pore network models of soils with known porosity values. Our
calculations showed that the relationship between the average pore
body and/or throat radius with d50 is nonlinear if the diﬀerence between d50 of the two main sands involved is very large. Similar conclusions were established in our earlier study on grain compaction(Mahmoodlu et al., 2016).

Once the complete SWRCs were obtained, we analyzed the curves interms of the standard unimodal van Genuchten function (1980) givenbyθ w (pc ) = θrw +θsw−θrw[1 + (αpc )n]1 − 1/ n

pc < 0(12)

where θ w is the volumetric water content, θrw denotes the residual
water content, θsw is the saturated water content (assumed in this study

to be equal to the porosity, ϕ), α and n are semi-empirical shapeparameters (further referred to as van Genuchten parameters), and Pc
denotes capillary pressure. Eq. (12) was ﬁtted to the SWRC data using
the RETC software (van Genuchten et al., 1991), which is a nonlinear
least-squares optimization program commonly used to for analyzing
observed unsaturated soil hydraulic data.
The SWRC in some cases, especially for the layered samples, could
be described better using a bimodal retention function. For those cases
we used the bimodal formulation introduced by Durner (1994):k

1⎤
Se w = ∑ wi ⎡⎢ 1 + (αi p )ni ⎥
c
⎦
⎣
i=1

(13)

where k is the number of subsystems that form the total pore-size
distribution (in this study, k = 2) and wi are weighting factors for the
sub-curves, subject to 0 < wi less than 1 and w1 + w2 = 1. Similarly as
for the unimodal curve, the parameters of the two sub-curves of Eq.
(13) were subjected to the constraints αi > 0, mi > 0, ni > 1.
4. Results
4.1. Eﬀect of heterogeneous packing on pore structures

4.2. Porosity and intrinsic permeability
Fig. 4a, b shows cumulative frequency distributions of the pore bodyand pore throat radii for Sands S1 to S5 (the main sands), as well as of

At the pore scale, porosity represents the ratio of the total volume ofFig. 4. Cumulative pore body and throat radii distributions for the ﬁve main sands (a and b), their mixtures with S1 (c and d), and the layered soil packings (e and f).
Each mixture was obtained by mixing 50% of sand S1 with 50% of one of the other sands. All simulations had a target porosity of 0.35.
575

Journal of Hydrology 565 (2018) 570–580

M.G. Mahmoodlu et al.

Fig. 5. Normalized values of the porosity at three locations in the S1-S2 layered packings. Zones A, B, and C are located in the middle of the ﬁrst layer, the interface ofthe two layers, and the middle of the top layer, respectively. Soil packing S1-S3 (a) represents two sands with a sharp interface and S3-S1 (b) a transitional interface.
ϕ* is the normalized porosity (porosity of each zone divided to the bulk porosity).

S1-S5 in order to keep the same thickness of each layer. However, thebulk permeability of a layered pack soil increased towards the coarse
sand.
The relationships between permeability and pore size of the ﬁve
main sands became slightly nonlinear towards coarse sand, as well as
for the mixtures of S1 with the other sands (Fig. 6b, c, e, f, i, g). The
relationships between permeability and pore size were even morenonlinear when the permeability at a given porosity of 0.35 is plotted asa function of either pore body radius or pore throat radius, as shown by
the results in Fig. 6b, c, e, f, i, g. Although, the permeability versus pore
body or throat radii showed nonlinear behavior for the sand with aninterface, the plots suggest linear slopes which changes after the second(S1-S3) packing (Fig. 6i, g).
The eﬀect of an interface on the bulk permeability of the layered
sand samples was also investigated. Fig. 1S (Supplementary material)
shows that a sharp interface causes an increase in the bulk permeabilityof the sand packing. This diﬀerence became more visible for thepackings of S1 with sand having larger d50 values. Fig. 7a depictschanges in the ratio of the bulk permeability to the average permeability of two sands, versus the ratio of d50 of sand S1 (when at thebottom of the sample) over d50 of the second sand (S2 to S5 at the top).
Results show that the ratio of the bulk permeability to average permeability is higher for the packings having a sharp interface (Fig. 7a).
This ratio diminished for the packing of sand S1 with a coarser sand. In
contrast with the interface eﬀect on the bulk permeability, the permeability of mixtures of two sands increased for the packing of S1 with acoarser sand (Fig. 7b).

pore bodies and pore throats to the total volume of the porous medium.
Hence, changes in porosity of a soil proﬁle are directly aﬀected by thechanges in the pore sizes of the soil. Fig. 5 depicts changes in theporosity of a layered sand pack within three diﬀerent zones (A, B, andC) along the sample proﬁle. As expected, the porosity was found toincrease at or near a sharp interface (Fig. 5a) as compared to the porosities of the two main sands individually (i.e., S1 and S3 in this case).
However, notice that the porosity at the transitional interface (Fig. 5b)
remained the same as that of the ﬁner sand (here Sand S1). This webelieve is due to the inﬁltration of small particles into the larger poresof the underlaying coarser layer, leading to smaller average pore sizes
at the interface.
Soil permeability is determined by both the bulk geometric properties of a sample (e.g., porosity and the pore size distribution), as wellas its topological properties such as pore space connectivity and tortuosity. We investigated the eﬀects of pore structure on the permeability of each individual sand (i.e., Sands S1 to S5), of the mixtures ofsand S1 (at a 0.5 fraction) with the other sands (S2 to S5), as well as ofthe layered soil packings. Results for the individual sands (Fig. 6a) showthat the permeability of the main sands (S1 to S5) increased considerably and nonlinearly towards the coarsest sand with its largest d50
value (i.e., sand S5 with a d50 of 885 µm). This conﬁrms that the permeability is strongly aﬀected by pore size (Håkansson and Lipiec, 2000;
Assouline, 2006, Mahmoodlu et al., 2016).
To investigate the eﬀect of grain heterogeneity on the permeability,
sand S1 was mixed with the other sands (S2, S3, S4 and S4). For all
calculations, the ratio of mixing was kept at 50% of S1 and 50% of theother sands. As shown in Fig. 6d, mixing S1 with a coarser sand (S2 toS5) caused an increase in the fraction of smaller pores of the coarser
sand and consequently a reduction in the permeability of sand S3. Bycontrast, the permeability of S2 increased due to an increase in thefraction of larger pores. We found that mixing sand S1 (having a d50 of
283 µm) with S2 (having a d50 of 385 µm) caused a small increase in thepermeability due to small diﬀerences in d50 values of the two sands,
while mixing the coarser sands (S3, S4, and S5 with d50 values of 485,
685 and 885 µm, respectively) increased the intrinsic permeability
substantially. Our results clearly show that the permeability reductionsof the coarser sand are dominated very much by the relative amounts of
the ﬁner sands (and hence the smaller pore sizes).
Changes in the bulk permeability, as a function d50, of the layered
sand packs having a sharp interface did not show a clear pattern(Fig. 6h). This may be because of the increased fraction of S1 in thelayered sands when coarser sands were used. As noted earlier, thefraction of S1 varied from 70% for packing S1-S2 to 95% for packing

4.3. Soil water retention curveThis section shows how sample heterogeneity aﬀects the soil water
retention curve (SWRC), which has a ﬁrst-order eﬀect on predictions ofwater ﬂow in variably saturated media. We investigated the eﬀect ofpore body and pore throat size as well as of the type of interface (sharp
or transitional) between two sands in the layered samples on the SWRC.
We ﬁrst studied the eﬀect of grain size distribution on the SWRC.
Fig. 8 depicts the evolution of water saturation, Sw (given by θw/θsw ) as afunction of capillary pressure, Pc, for the ﬁve main sands. Results indicate that the capillary pressure at a given saturation was always
higher for the ﬁner sands due to the creation of smaller pore throats
regulating higher capillary entry pressures and controlling capillary
displacement (see S1 in Fig. 8).
For the heterogeneity in terms of the grain size distribution simulations, combining S1 (fraction of 0.5) with a coarser sand (fraction of576

Journal of Hydrology 565 (2018) 570–580

M.G. Mahmoodlu et al.

Fig. 6. Changes in the intrinsic permeability with d50 and the average pore body (APBR) and average pore throat (APTR) radii of the ﬁve main sands (a, b, c left
column), of the mixtures of 50% of S1 with 50% of S2, S3 and S4 with S5 (d, e, f, middle column), and of the layered samples with a sharp interface (h, i, g, rightcolumn).

0.5) caused a decrease in the capillary pressure at a given volumetricwater content. This hence produced larger values of the van Genuchten
α parameter in Eq. (3), and hence lower entry pressures (Fig. 9). The
SWRCs of the combination of sand S1 with the other sands changed
more non-uniformly towards the coarser sands. This is because of thecreation of larger pore throat radii and their controlling eﬀect on the
SWRC of the mixture of S1 with coarser sand.
We next investigated the eﬀects of sharp versus transitional interfaces on the SWRCs (Fig. 10). Calculated SWRCs of the layered sandpackings with a sharp interface showed a clear eﬀect of the interface
when d50 of two sands diﬀered substantially from each other. By
comparison, the SWRC of the sample with S1 (at the bottom) and S2 (atthe top) showed a more monotonic change and did not exhibit a clear
discontinuity at the interface of the two sands (Fig. 10). This is because
the capillary pressures at which the water-ﬁlled pores start to drain are
directly proportional to the equivalent radii of the pore throat (e.g.,
Durner, 1994). The radii of the throats for sands S1 and S2 were fairly
close to each other, which leads to a monotonic change in the SWRC.
Since the radii of the pore throats of S1 and the other of sands (e.g., S3,
S4, and S5) were not in a same range, their mixtures should show a

Fig. 8. Soil water retention curves of the ﬁve main sands.

more deﬁnite jump at the interface and less smooth behavior of theSWRCs.
As shown in Fig. 10, the drainage curves of the soil samples with atransitional interface (coarse sand at the bottom) showed relatively
high residual water saturations, which increased for a layered sandFig. 7. Eﬀects of sharp and transitional interfaces on the intrinsic permeability of the layered samples of S1 and the other sands. Results are plotted as a function ofd50(S1)/d50(Sx), where d50 is the mean diameter of each sand, and Sx represents Sands S2 to S5. Kc is the intrinsic permeability of the layered soil sample, Km intrinsic
permeability of the sand mixture, and Kave the average intrinsic permeability of the two sands involved. Calculations for all soil packings assumed a porosity of 0.35.
577

Journal of Hydrology 565 (2018) 570–580

M.G. Mahmoodlu et al.

Fig. 9. Eﬀect of mixing on the soil water retention curve.

the sample.
Results of both main sands and their mixtures with S1 showed anincrease in the van Genuchten α value towards sand having a larger
value of d50 (Table 2). The ratio of the van Genuchten parameter α ofthe sand mixture (αmix) to the averages value of α (αave) is plotted inFig. 11a versus the ratio of d50 of S1 to d50 of other sands. The plotshows a linear increase in the ratio of αmix to αave when S1 is mixed withcoarser material (i.e., the ratio of the d50 values of the two sands).
Fig. 11b, c further depicts for the sharp interfaces the ratio of thepack having coarser sand (S5) at the bottom. Moreover, the SWRCs of
the layered sands with a sharp interface showed very high van Genuchten n values and residual water saturations (Table 2). In our simulations, air always invaded the soil sample from the top boundary. A
possible explanation for the high residual saturations could be the highcapillary pressure needed to drain the ﬁner layer at the top (here S1).
However, as soon as an air pathway reaches the interface, air will moveinto the coarser sand and thus bypass several upstream saturated pores.
This may result in a constant, but non-uniform water saturation along

Fig. 10. Eﬀect of sharp and transition interfaces on the soil water retention curve.
578

Journal of Hydrology 565 (2018) 570–580

M.G. Mahmoodlu et al.

Fig. 11. Eﬀect of heterogeneity on the van Genuchten parameter α. Plot a) is for mixtures of S1 with the other sands, b) and c) represent sharp interfaces (with S1
located at the bottom), while d) belongs to the transitional interfaces (with the coarser sands at the bottom and S1 at the top). Results are given as a function ofd50(S1)/d50(Sx), where d50 is the mean diameter of each sand, and Sx represents Sands S2 to S5. Calculations for all soil packings assumed a porosity of 0.35.

van Genuchten parameter α1 and αave versus the ratio of d50 of S1 to d50of the other sands. We used α1 here for the coarser sand in the layered
sample, while αave is again the average of the α values of the two sands
involved. Similarly as for mixing scenario, a linear relationship existed
between the ratios of α1 to αave versus the ratio of d50 of S1 to d50 of theother sands. However, we found a nonlinear relationship between theratios α1/αave (for the layered-samples) to the ratio of αmix to αave (forthe sand mixtures), whereas the ratio α2/αave increased for two sandswith a large ratio of d50. In contrast with the sharp interfaces, the ratio
α1/αave versus the ratio d50(S1)/d50(Sx) (where Sx denotes again S2 to S5)
did not show linear behavior, but more of an inverse normal distribution.

of two unimodal sands. In general, the ratio of bulk permeability to theaverage permeability of two sands used in a two-layered packing increased more for the packings having a sharp interface.
Calculated SWRCs showed that mixtures of the ﬁner sand with acoarser sand caused a decrease in the capillary pressure of the coarser
sand. This in turn produced larger values of the van Genuchten parameter α and hence lower entry pressures. We found the ratio of the α (ofthe mixtures) to the average value of α to increase with larger diﬀerences in the d50 values of the two sands involved.
Sharp interfaces can cause a non-monotonic change in the drainage
curve due to discrepancies between the pore throat radii of the twosands in the sample. This change is inﬂuenced by d50 or the grain size
distribution of the sands used in the two-layered packing. In contrast
with sharp interfaces, transitional interfaces produced monotonic
changes in the drainage curves for two sands having a small diﬀerencein their d50 values.
High residual water saturations were found for the drainage curvesof sand samples having a transitional interface (with coarse sand at thebottom). Residual saturations increased with the d50 value of the sandat the bottom of the sample. Similarly as for the mixture scenario, theratio of the van Genuchten parameter α to the average value of α of thetwo sands involved, increased with larger diﬀerence in the d50 values ofthe sands.
Finally, we note that our study involved artiﬁcially created porous
media made up of spherical particles. The results give much insight in
how particle mixing and layering aﬀect the unsaturated soil hydraulicproperties. We acknowledge that comparisons with actual experimental
data are very much needed to more deﬁnitely conﬁrm our ﬁndings.

5. Conclusions
In this study we used the discrete element method (DEM) to explore
the eﬀects of soil heterogeneity on pore structure and the unsaturated
soil hydraulic properties (the intrinsic permeability and the soil water
retention curve, SWRC) of ﬁve diﬀerent sands and their mixtures. We
found, as expected, that mixing of the ﬁner sand with any of the coarse
sands causes a reduction in the average pore body and pore throat sizes,
and subsequently will diminish the intrinsic permeability values of thecoarser sands.
For the layered porous media with a transitional interface, the local
porosity of the soil decreased. This suggests that inﬁltration and sedimentation of small particles into larger pores of the coarser sands (typical of a graded bed in geological formations) decreases the soil poresizes at the interface, and consequently reduces the local porosity of thesoil at the interface. In contrast, and interestingly, the local porosity of
the layered media with sharp interfaces increased due to deposition oflarger grains on top of the ﬁne grains. Results conﬁrmed the development of a discontinuity between two layers through the creation ofrelatively large pores at the interface.
Non-linear relationships between permeability and the average pore
body and throat radii were found when mixing diﬀerent unimodal
sands. One likely reason is variation in the soil pore sizes of the mixture

Acknowledgments
The authors would like to thank Hamed Aslannejad (Utrecht
University) for his technical assistance. Comments by four anonymousreferees and the Associate Editor of Journal of Hydrology further
helped to improve the manuscript.
579

Journal of Hydrology 565 (2018) 570–580

M.G. Mahmoodlu et al.

Appendix A. Supplementary data

26-28 Nov. 2008. Schriftenreihe Geotechnik 21. Bauhaus-Universität, Weimar,
Germany, pp. 53–60.
Mohanty, B.P., Zhu, J., 2007. Eﬀective hydraulic parameters in horizontally and vertically heterogeneous soils for steady-state land-atmosphere interaction. J.
Hydrometeorol. 8 (4), 715–729.
Morbidelli, R., Saltalippi, C., Flammini, A., Corradini, C., Brocca, L., Govindaraju, R.S.,
2016. An investigation of the eﬀects of spatial heterogeneity of initial soil moisture
content on surface runoﬀ simulation at a small watershed scale. J. Hydrol. 539,
589–598.
Mousavi, M.A., Bryant, S.L., 2012. Connectivity of pore space as a control on two-phaseﬂow properties of tight-gas sandstones. Transport Porous Med. 94 (2), 537–554.
https://doi.org/10.1007/s11242-012-0017-x.
Mukhlisin, M., Baidillah, M.R., Taha, M.R., El-Shaﬁe, A., 2011. Eﬀect of soil water retention model on slope stability analysis. Int. J. Phys. Sci. 6 (19), 4629–4635.
Nimmo, J.R., Perkins, K.S., Lewis, A.M., 2002. Steady-State Centrifuge, Methods of Soil
Analysis. In: Inc, S.S.S.A. (Ed.), Part 4. Madison, WI, pp. 903–916.
Perrier, E., Rieu, M., Sposito, G., Marsily, G., 1996. Models of the water retention curvefor soils with a fractal pore size distribution. Water Resour. Res. 32 (10), 3025–3031.
Peth, S., Nellesen, J., Fischer, G., Horn, R., 2010. Non-invasive 3D analysis of local soil
deformation under mechanical and hydraulic stresses by mCT and digital image
correlation. Soil Tillage Res. 111, 3–18. https://doi.org/10.1016/j.still.2010.02.007.
Prodanović, M., Mehmani, A., Sheppard, A.P., 2015. Imaged-based multiscale networkmodelling of microporosity in carbonates. Geological Society, London, Special
Publications 406 (1), 95–113.
Purushothama, R.P., 2007. Soil Mechanics and Foundation Engineering, Second edition.
Pearson Education India, 840 pages.
Raoof, A., Hassanizadeh, S.M., 2012. A new formulation for pore-network modeling of
two-phase ﬂow. Water Resour. Res. 48 (1), 1–13.
Raoof, A.S.M., Hassanizadeh, S.M., Leijnse, A., 2010. Upscaling transport of adsorbing
solutes in porous media: pore-network modeling. Vadose Zone J. 9 (3), 624–636.
https://doi.org/10.2136/vzj2010.0026.
Sakaki, T., Smits, K.M., 2015. Water retention characteristics and pore structure of binary
mixtures. Vadose Zone J. 14 (2). https://doi.org/10.2136/vzj2014.06.0065.
Scanlon, B.R., Andrasky, B.J., Bilskie, J., 2002. Miscellaneous methods for measuringmatric or water potential. In: Dane, J.H., Topp, G.C. (Eds.), Methods of soil analysis.
Part 4. Physical methods. SSSA Book Ser. 5. SSSA, Madison, WI, pp. 643–670.
Šmilauer, V., Catalano, E., Chareyre, B., Dorofeenko, S., Duriez, J., Gladky, A., Kozicki, J.,
Modenese, C., Scholtès, L., Sibille, L., Stránský, J., Thoeni. K., 2015. Yade
Documentation (V.Šmilauer, ed.), The Yade Project, 1st ed.
Sweijen, T., Nikooee, E., Hassanizadeh, S.M., Chareyre, B., 2016. The eﬀects of swellingand porosity change on capillarity: DEM coupled with a pore unit assembly method.
Transport Porous Med. 113 (1), 207–226.
Sweijen, T., Aslannejad, H., Hassanizadeh, S.M., 2017. Capillary pressure-saturation relationships for porous granular materials: pore morphology method vs. pore unit
assembly method. Adv. Water Resour. 107, 22–31.
Tong, A.T., Catalano, E., Chareyre, B., 2012. Pore-scale ﬂow simulations: model predictions compared with experiments on bi-dispersed granular assemblies. Oil Gas Sci.
Technol. 67, 743–752. https://doi.org/10.2516/ogst/2012032.
Torskaya, T., Shabro, V., Torres-Verdín, C., Salazar-Tio, R., Revil, A., 2014. Grain shapeeﬀects on permeability, formation factor, and capillary pressure from pore-scalemodeling. Transport Porous Med. 102 (1), 71–90.
Unsal, E., Mason, G., Morrow, N.R., Ruth, D.W., 2009. Bubble snap-oﬀ and capillary backpressure during counter-current spontaneous imbibition into model pores. Langmuir25, 3387–3395.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898.
van Genuchten, M.Th., Leij, F.J., Yates, S.R., 1991. The RETC code for quantifying thehydraulic functions of unsaturated soils. Report No. EPA/600/2-91/065. R.S. Kerr
Environmental Research Laboratory, U. S. Environmental Protection Agency, Ada,
OK. 85 p. https://www.pc-progress.com/en/Default.aspx?retc).
Vasilyev, L., Raoof, A., Nordbotten, J.M., 2012. Eﬀect of mean network coordination
number on dispersivity characteristics. Transport Porous Med. 95 (2), 447–463.
Vaz, C.M.P., Tuller, M., Lasso, P.R.O., Crestana, S., 2014. New perspectives for the application of high-resolution benchtop X-ray microCT for quantifying void, solid andliquid phases in soils. In: W.G. Teixeira et al., editors, Application of soil physics inenvironmental analyses: Measuring, modelling and data integration. Progr. Soil Sci.
Springer Int. Publ., Cham, Switzerland. 261-281. doi:10.1007/978-3-319-060132_12.
Yuan, C., Chareyre, B., Darve, F., 2015. Pore-scale simulations of drainage in granularmaterials: ﬁnite size eﬀects and the representative elementary volume. Adv. Water
Resour. 95, 109–124.
Zhang, Q., Raoof, A., Hassanizadeh, S.M., 2015. Pore-scale study of ﬂow rate on colloid
attachment and remobilization in a saturated micromodel. J. Environ. Qual. 44 (5),1376–1383.
Zhu, J., Sun, D., 2012. Soil hydraulic properties for moisture redistribution in a largescale heterogeneous landscape. Hydrol. Sci. J. 57 (6), 1196–1206.
Zhuang, L., Hassanizadeh, S.M., Qin, C.Z., de Waal, A., 2017. Experimental investigationof hysteretic dynamic capillarity eﬀect in unsaturated ﬂow. Water Resour. Res. 53(11), 9078–9088.

Supplementary data associated with this article can be found, in the
online version, at https://doi.org/10.1016/j.jhydrol.2018.07.060.
References
Al-Ismaily, S.S., Al-Maktoumi, A.K., Kacimov, A.R., Al-Saqri, S.M., Al-Busaidi, H.A., 2015.
Impact of a recharge dam on the hydropedology of arid zone soils in Oman: anthropogenic formation factor. J. Hydrol. Eng. 20 (4). https://doi.org/10.1061/(ASCE)HE.1943-5584.0000886.
Alonso, E.E., Cardoso, R., 2010. Behavior of materials for earth and rockﬁll dams: perspective from unsaturated soil mechanics. Front. Archit. Civil Eng. China 4 (1), 1–39.
Antinoro, C., Arnone, E., Noto, L.V., 2017. The use of soil water retention curve models inanalyzing slope stability in diﬀerently structured soils. Catena 150, 133–145.
Arora, B., 2012. Modeling Biogeochemistry and Flow within Heterogeneous Formations
in Variably-Saturated Media. Doctoral Dissertation. A&M University, Texas.
Assouline, S., 2006. Modeling the relationship between soil bulk density and the hydraulic conductivity function. Vadose Zone J. 5, 697–705.
Belheine, N., Plassiard, J.P., Donzé, F.V., Darve, F., Seridi, A., 2009. Numerical simulationof drained triaxial test using 3D discrete element modeling. Comput. Geosci. 36,
320–331.
Bittelli, M., Flury, M., 2009. Errors in water retention curves determined with pressure
plates. Soil Sci. Soc. Am. J. 73, 1453–1460.
Chareyre, B., Cortis, A., Catalano, E., Barthélemy, E., 2012. Pore-scale modeling of viscous ﬂow and induced forces in dense sphere packings. Transport Porous Med. 94 (2),595–615.
Chen, P.Y., Chen, C.H., Hsu, N.S., Wu, C.M., Wen, J.C., 2012. Inﬂuence of heterogeneity
on unsaturated hydraulic properties: (1) local heterogeneity and scale eﬀect. Hydrol.
Process. 26 (23), 3593–3603.
Cnudde, V., Boone, M.N., 2013. High-resolution X-ray computed tomography in geosciences: a review of the current technology and applications. Earth Sci. Rev. 123,
1–17. https://doi.org/10.1016/j.earscirev.2013.04.003.
Coppola, A., Basile, A., Comegna, A., Lamaddalena, N., 2009. Monte Carlo analysis ofﬁeld water ﬂow comparing uni-and bimodal eﬀective hydraulic parameters forstructured soil. J. Contam. Hydrol. 104 (1), 153–165.
Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies.
Geotechnique 29, 47–65. https://doi.org/10.1680/geot.1979.29.1.47.
Dane, J.H., Hopmans, J.W., 2002. Hanging water column. p. 680-684. In: J.H. Dane andG.C. Topp (ed.) Methods of soil analysis. Part 4. SSSA Book Ser. 5. SSSA, Madison, WI.
Di Prima, S., Bagarello, V., Bautista, I., Burguet, M., Cerdà, A., Iovino, M., Prosdocimi, M.,
2016. Comparing a simple methodology to evaluate hydrodynamic parameters with
rainfall simulation experiments. EGU General Assembly 2016, held 17-22 April, 2016
in Vienna Austria, p.289.
Dubois, V., Abriak, N.E., Zentar, R., Ballivy, G., 2009. The use of marine sediments as apavement base material. Waste Manage. 29 (2), 774–782.
Durner, W., 1994. Hydraulic conductivity estimation for soils with heterogeneous pore
structure. Water Resour. Res. 30 (2), 211–223.
Falode, O., Manuel, E., 2014. Wettability eﬀects on capillary pressure, relative permeability, and irredcucible saturation using porous plate. J. Pet. Eng. ID 465418, 12p.
https://doi.org/10.1155/2014/465418.
Gee, G., Or, D., 2002. Particle size analysis. In: J. J. Dane and G. C. Topp (eds.), Methodsof Soil Analysis, Part 4. Physical Methods, 255-294, Soil Sci. Soc. Am., Madison, WI.
Gladkikh, M., Bryant, S., 2005. Prediction of imbibition in unconsolidated granular materials. J. Colloid Interface Sci. 288, 526–539.
Håkansson, I., Lipiec, J., 2000. A review of the usefulness of relative bulk density values
in studies of soil structure and compaction. Soil Till. Res. 53 (2), 71–85.
Hohenbrink, T.L., Lischeid, G., 2015. Does textural heterogeneity matter? quantifying
transformation of hydrological signals in soils. J. Hydrol. 523, 725–738.
Huang, M., Lee Barbour, S.L., Elshorbagy, A., Zettl, J.D., Cheng Si, B., 2011. Inﬁltrationand drainage processes in multi-layered coarse soils. Can. J. Soil Sci. 91 (2), 169–183.
https://doi.org/10.4141/CJSS09118.
Katukiza, A.Y., Ronteltap, M., Niwagaba, C.B., Kansiime, F., Lens, P.N.L., 2014. Grey
water treatment in urban slums by a ﬁltration system: optimisation of the ﬁltration
medium. J. Environ. Manage. 146, 131–141.
Keller, T., Lamandé, M., Peth, S., Berli, M., Delenne, J.Y., Baumgarten, W., Rabbel, W.,
Radjai, F., Rajchenbach, J., Selvadurai, A.P.S., Or, D., 2013. An interdisciplinary
approach towards improved understanding of soil deformation during compaction.
Soil Till. Res. 128, 61–80.
Lopes, R.L., Jucá, J.F.T., Mariano, M.O.H., 2013. Water inﬁltration in ﬁnal cover layer of
landﬁlls in northeast region of Brazil. Adv. Unsaturated Soils. CRC Press 243–248.
https://doi.org/10.1201/b14393-33.
Mahmoodlu, M.G., Raoof, A., Sweijen, T., van Genuchten, M.T., 2016. Eﬀects of sandcompaction and mixing on pore structure and the unsaturated soil hydraulic properties. Vadose Zone J. 15 (8). https://doi.org/10.2136/vzj2015.10.0136.
Mehlhorn, T., Prohaska, S., Homberg, U., Slowik, V., 2008. Modelling and analysis ofparticle and pore structures in soils. In: Witt, K.J. (Ed.), Workshop Internal Erosion.

580