Published August 11, 2016

Original Research

Effects of Sand Compaction
and Mixing on Pore Structure
and the Unsaturated Soil
Hydraulic Properties
M.G. Mahmoodlu,* A. Raoof, T. Sweijen,
and M. Th. van Genuchten

Core Ideas• A discrete element model was used
to extract the pore structure of particle packings.
• Soil grain compaction and size mixing affect pore structure.
• Size mixing had a greater effectthan compaction on the pore size
distribution.
• Compaction caused the SWRCs tochange uniformly with porosity.
• Compaction decreased the van
Genuchten a parameter but had
less effect on n.

The hydraulic properties of unsaturated porous media very much dependon their pore structure as defined by the size, arrangement, and connectivity of pores. Several empirical and quasi-empirical approaches have beenused over the years to derive pore structure information from the particle
size distribution. In this study, we used the discrete element method to simulate the pore structure of various sands as affected by compaction andparticle mixing processes. We used five sands with different mean grain sizes
to investigate the effects of different sand mixing ratios and degrees of compaction on pore structure as well as on the intrinsic permeability and thesoil water retention curve. Average pore body and pore throat sizes were
found to be determined mostly by the smaller particles as represented by
the effective diameter D10. The effects of compaction on the average pore
body and pore throat radii were used to simulate expected decreases in
the permeability. We obtained mostly linear relationships between permeability and the average pore body and throat radii when mixing different
unimodal sands. The intrinsic permeability of the coarser sands was found tobe far more sensitive to porosity than the finer sands. Simulations of unsaturated conditions showed that the van Genuchten hydraulic parameter aincreased nonlinearly with increasing grain size and mean pore body size ofthe sand mixtures. Compaction caused a linear decrease in a with decreasing porosity and pore body size. However, no clear correlation between thevan Genuchten parameter n and porosity or D10 was found for the different
compaction and mixing simulations.
Abbreviations: DEM, discrete element method; SWRC, soil water retention curve

M.G. Mahmoodlu, Dep. of Earth Sciences, Utrecht Univ., Utrecht, theNetherlands, and Dep. of Watershed
and Rangeland Management, Gonbad Kavous Univ., Iran; A. Raoof and
T. Sweijen, Dep. of Earth Sciences,
Utrecht Univ., Utrecht, the Netherlands;
M.Th. van Genuchten, Dep. of Earth
Sciences, Utrecht Univ., Utrecht, theNetherlands, and Dep. of Mechanical Engineering, Federal Universityof Rio de Janeiro, Rio de Janeiro, RJ,
Brazil.
*Corresponding
author(m.g.mahmoodlu@gmail.com).

Vadose Zone J.
doi:10.2136/vzj2015.10.0136
Supplemental material online.
Received 5 Oct. 2015.
Accepted 19 May 2016.
Vol. 15, Iss. 8, 2016
© Soil Science Society of America

5585 Guilford Rd., Madison, WI 53711 USA.
All rights reserved.

Temporal changes in the soil hydraulic properties due to compaction can affect a

range of near-surface and subsurface flow processes. Several natural and human-induced
processes such as the use of farm machinery in agricultural operations, grazing activities,
raindrop impact, timber harvesting in forests, fire, and various geotechnical engineering projects can lead to soil compaction (e.g., Richard et al., 2001a, 2001b; Assouline,
2006a; Keller et al., 2013; Kuncoro et al., 2014). Soil compaction alters pore spaces andconsequently affects the bulk density, porosity, and soil hydraulic properties. Many studies over the years have shown that soil compaction and related surface sealing problems at
the field scale lead to reductions in (i) infiltration and recharge rates, (ii) soil aeration andthe gas diffusivity, (iii) the efficiency of fertilizer use, (iv) seedling emergence, and (v) thegrowth and distribution of plant roots (e.g., Bradford et al., 1987; Richard et al., 2001a,
2001b; Hadas, 2004; Wall and Heiskanen, 2009; Glab and Kopec, 2009; Vahyala et al.,2013; Kuncoro et al., 2014; Gregorich et al., 2014; Masís-Meléndez et al., 2015; Berli etal., 2015; Sela et al., 2015). Compaction and the related processes of pore clogging and soil
deformation were comprehensively reviewed by Assouline (2006a) and Keller et al. (2013).
Changes in the soil physical and hydraulic properties may also occur in the absence of compaction. For example, fine sediment infiltration can lead to the deposition of small-sized

Vadose Zone Journal | Advancing Critical Zone Sciencesediments into relatively large near-surface soil pores, thus leading to low-permeability surface horizons (Cui et al., 2008; Evans
and Wilcox, 2014) and reducing the rates of water flow into andthrough the soil profile (Haghnazari et al., 2015). On the other
hand, reductions in water percolation due to fine sediment infiltration may be advantageous for certain hydraulic structures suchas dam reservoirs and qanat tunnels, as well as promoting water
harvesting in dry climates by increasing runoff.
Transport parameters determine flow and solute transport throughthe soil column (Mahmoodlu et al., 2014, 2015). Soil compaction and fine sediment infiltration tend to reduce the porosity,
hydraulic conductivity, and soil water retention. Pore-scale studies are important to understand the basic processes causing thesechanges in the soil hydraulic properties. Most previous studies on
soil compaction and pore structure have involved quasi-empiricalapproaches and the use of pedotransfer functions (e.g., Rawls etal., 1983; Rajkai et al., 1996) or expressing selected soil water
retention parameters as a function of bulk density or porosity(e.g., Baumhardt et al., 1990; Assouline, 2006a). An overview ofthe more refined analytical and numerical modeling approaches
related to soil deformation and compaction was provided by Keller
et al. (2013).
Studies of compaction and grain mixing effects on the unsaturated
soil hydraulic properties require precise studies of the geometry ofsand grains and related pore structures in terms of the prevailingpore body and pore throat sizes. Although the soil particle size distribution is relatively easy to measure using a range of techniques(Gee and Or, 2002), estimation of grain positions and pore structures cannot be done in the same way (Mehlhorn et al., 2008).
One direct method to measure pore structures is the use of X-ray
microtomography, which can provide detailed three-dimensional
images of the medium (e.g., Peth et al., 2010; Cnudde and Boone,
2013; Vaz et al., 2014). However, X-ray microtomography is still
relatively expensive and can serve only as a first step for identifying the pore structure. Image processing is the next step to obtainpore structure geometries, which often requires further analysis(Mehlhorn et al., 2008).
An alternative method is to consider idealized grain geometries
to obtain pore structures for a given particle size distribution(Mehlhorn et al., 2008; Keller et al., 2013). In this study, we usedthe discrete element method (DEM) to simulate the movementof grains and to obtain the resulting pore structures for a varietyof packings. For this purpose, we used the open source platformYade-DEM (Šmilauer et al., 2015), which also enables computation of the hydraulic properties as well as the soil water retention
curve (SWRC). The specific objectives of this study were to: (i)use the discrete element model to extract pore structures ofdifferent particle packings; (ii) explore the effects of soil graincompaction and mixing on pore structures; (iii) investigate theeffects of compaction and mixing on the intrinsic permeability;

VZJ | Advancing Critical Zone Scienceand (iv) evaluate the effects of compaction and mixing on the
SWRCs, including the van Genuchten hydraulic parameters aand n. In our approach, we kept the grain size distribution constant while increasing the level of compaction, thus reducing theporosity value and increasing the bulk density. For the mixing
simulations, on the other hand, we kept the total porosity constant while varying the particle size distribution.

66Sand Properties and Packings

As input to the discrete element models, we used the particle size
distributions of five filter sands (Fig. 1). The sands originated froma river bed in Papendrecht (the Netherlands) and were sieved todifferent retained sizes as given in Table 1. The mean particle size
varied from 350 to 1300 mm (Table 1). X-ray fluorescence analyses
showed that >96% of the sand grains consisted of silica (SiO2).
Microscopic images of the sands showed a good degree of roundness of the sand grains (Fig. 1), which is an advantage for using
the DEM. Coefficients of uniformity and curvature, given by
Cu = D60/D10 and Cc = (D30)2/(D10D60), respectively, were estimated for all sands (Table 1). Because Cu was <4 for all sands, thesamples can be considered to be uniformly graded by containing
essentially identical particle sizes. Sand S2 was the most uniform(Fig. 1), and hence used as the main sand for mixing with the othersamples. We considered the coefficient of curvature (Cc) as anotheruseful measure of the shape of the particle size distribution curve.
This coefficient is sometimes used to differentiate the degree ofsorting of a soil. Sand S1 showed a lesser degree of sorting, as
reflected in its bimodal grain size distribution (Fig. 1).
In view of its favorable Cu and Cc values, we used Sand S2 as themain sand for our study. We considered various combinations of
S2 with the other sands, as shown in Fig. 2, to explore the effectsof mixing. A total of 13 samples were considered: the five main

Fig. 1. Shape and diameter of the main sand particles (top) and their
particle size distributions (bottom).

p. 2 of 11

Table 1. The mean size, effective and mean diameters (D10 and D50,
respectively), coefficients of uniformity and curvature (C­u­ and Cc,
respectively), pH, and total organic C (TOC) in dry samples of thenatural sands used in this study.
Sand Mean size

D10

D50

Cu

Cc

pH

TOC
g m−3

—————— mm ——————
S1

200–500

290

390

1.140

1.063

S2

400–630

425

500

1.100

1.001

S3

500–1000

540

710

1.187

0.983

S4

700–1250

800

990

1.132

1.005

S5

1000–1600

1140

1360

1.102

1.020

7

<5.0 ´ 10−20

sands (S1–S5; Table 1) and eight combinations of S2 with the othersands (Fig. 2). The porosities of all main sand samples were set at
0.40. The different sand grain combinations were used to calculate values of D10 (referred to as the effective diameter), D30 (thegrain diameter corresponding to 30% w/w passing), D50 (median
diameter), and D60 (the grain diameter corresponding to 60% w/w
passing) of the mixtures. Figure 2 shows the calculated D10 and
D50 values of the main sand (S2) and the various combinations. As
expected, D10 and D50 both increased for the grain size distributions of the coarser sands and mixtures. We next used the particle
size distributions of the different sands (Table 1) and their mixtures (Fig. 2) to extract their pore structures.
To observe the evolution of soil properties during compaction,
we computed for each sample their pore structure geometries atfive different porosities. Porosity values during compaction wereassumed to decrease from 0.40 to 0.30 at intervals of 0.025. For
each of the resulting porosities, we then calculated the pore structure, the SWRC, and the intrinsic permeability.

6The
6
Discrete Element Method
General Principles

Discrete element models are being used increasingly in soil androck mechanics to explore deformation and shear strength properties of granular materials (Cundall and Strack, 1979; Scholtès
et al., 2009; Scholtès and Donzé, 2012) as well as for agriculturaland industrial applications such as soil tillage (Xu et al., 2002;
Shmulevich et al., 2007) and particle flows in silos (Coetzee andEls, 2009). In the DEM, grain particles are considered as spheres
defined by their radius, mass, momentum, and acceleration. A large
number of spheres are typically used to represent a particular soil
packing. The approach assumes that a spring-type behavior occurs
at the contact between spheres. Discrete element models for thisreason are often referred to as spring-dashpot models (Cundall andStrack, 1979). In our study, we used the Yade-DEM software tosimulate packings of particles for varying target porosities having
different particle size distributions.

VZJ | Advancing Critical Zone ScienceFig. 2. The effective diameter (D10) and mean diameter (D50) of thedifferent combinations of sand grains from sands S1 to S5.

We considered normal elasticity, tangential elasticity, and sliding
as the main processes taking place at the contact points. Following
a study by Labenne sand by Belheine et al. (2009), normal elasticity was calculated according to a linear contact law. When two
particles are in contact and pushed toward each other by the surrounding medium, they are flattened at their contact point. The
tendency of both particles to keep their initial shapes may cause
an elastic force fn:
f n = kn d n[1]

where kn is a coefficient describing the stiffness of the particles andd n denotes the normal displacement.
As a result of the shear at the contact point, a tangential force, ft,
may arise. When the tangential force reaches a certain thresholdvalue, the contact point will break, causing sliding of the two particles. This process is described by

f t £ f n tan f[2]

where f is the friction angle used to estimate the threshold value.
For each particle, a force balance is made by including forces due
to contacts with other particles, the boundary conditions, and theeffects of gravity. The resulting net force gives the acceleration of
the particle, using Newton’s second law. We refer to studies by
Cundall and Strack (1979) and Belheine et al. (2009) for details.

Packing Simulations
The DEM was used to generate a packing of spherical particles for acertain porosity and particle size distribution. Initially, we selected
5000 particles from a given particle size distribution. These wereinserted in a large box such that a cloud of particles was generated.
The grain-scale parameters were taken from Belheine et al. (2009),
p. 3 of 11

representing a dune sand. The values of particle stiffness, the friction angle, and Poisson’s ratio were 9.6 ´ 108 Pa, 30.0°, and 0.04,
respectively. All six boundaries of the domain were subsequently
allowed to confine the particles until a porosity value of 0.45 wasreached. After this point, only the top boundary was allowed tomove down to further confine the particle packing until a target
porosity was reached. From the resulting packing, the pore bodyand pore throat sizes were extracted using the pore unit method.
Moreover, the resulting pore spaces were used to calculate the
SWRC and the permeability. The final packing resulted in around9500 pore throats and 11,500 pore units.

Calculation of Pore Structure
A regular triangulation method was used to extract the geometryof the pore structure (Chareyre et al., 2012; Yuan et al., 2015). In
the application of this method, the pore space is divided into tetrahedrons, with each tetrahedron serving as a pore unit. As shownin Fig. 3, each tetrahedron is surrounding one pore body and has
four vertices that are located at the centers of the spherical particles. The sides, or facets, of a tetrahedron are the locations at
which two pore bodies are connected. The facet is the narrowest
transect between two connecting pore units and is referred to as
the pore throat. While some studies have suggested random poreconnectivities due to the natural structure of porous media (Raoofand Hassanizadeh, 2010; Vasilyev et al., 2012; Li et al., 2014), ourapproach assumed that four pore throats are connected to each pore
body. Each pore throat in this way is surrounded by three particles
to which an inscribed circle can be fitted (Fig. 3c). 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. The inscribed sphere is referred toas the Haines insphere (Gladkikh and Bryant, 2005; Unsal et al.,2009; Mousavi and Bryant, 2012; Prodanović et al., 2015), while
its corresponding radius is denoted as Ri. Compaction was incorporated in the pore-scale calculations by repeated repacking andrearrangement of the spherical particles to decrease porosities to theprescribed values, with concomitant changes in the pore structures.

Calculation of Hydraulic Properties
The permeability of the granular material was computed using
the pore finite volume method as implemented in the DEM byChareyre et al. (2012). This method allows solving for the pressure

distribution under fully water-saturated conditions. The approachassumes that one pressure value can be associated with each pore
unit. The flux, qij [L3 T−1], of water from the ith pore unit to thejth pore unit is given by (Vasilyev et al., 2012)
q ij = kijpi - p j
l ij

where pi and pj [M L−1 T−2] are the pressures in the ith and jth poreunits, respectively, lij [L] is the length between the centers of the
ith and jth pore units, and kij [M−1 L5 T] denotes the hydraulicconductivity of the facet (pore throat) between the ith and jthpores, given bykij =Aijf Rij22m

VZJ | Advancing Critical Zone Science

[4]

where Aijf [L 2] is the smallest transect of the pore throat, m[M L−1] is the dynamic viscosity, and Rij [L] is the hydraulic radiusof the facet:
Rij =

Q ij
w ij

[5]

where Qij [L3] is the pore throat volume and w ij [L2] the surface ofthe pore throat. See Chareyre et al. (2012) for more details.
Under steady-state conditions, and assuming incompressible flow,
the continuity equation provides that (Zhang et al., 2015)4

å qij = 0

[6]

j =1

Equations [3], [4], and [6] were solved implicitly to calculate pressure values (Raoof et al., 2010). The total flux through the medium
was subsequently used to compute the permeability of the medium.

Capillary Pressure–Saturation Calculations
To estimate the capillary pressure–saturation (soil water retention) curve, the bottom boundary of the simulated domain was
assumed to be a water reservoir at a fixed pressure, Pw, whereas the
top boundary was considered to be an air reservoir with pressure
Pair. Under quasi-static conditions, the capillary pressure, Pc, isthen equal to the pressure difference between the two reservoirs:

Pc = Pair - Pw

Fig. 3. Schematic of the (a) spherical packing with (b) a pore body withradius Ri and (c) a pore throat with radius Rij as used in this study.

[3]

[7]

For drainage (i.e., drying) simulations, the non-wetting air pressure
was increased in steps to allow air to invade increasingly smaller
pore throats. If the capillary pressure is larger than the entry pressure of a pore throat, air could invade a water-saturated pore unit.
The entry pressure (Pe) of a pore throat can be obtained using theYoung–Laplace equation:
p. 4 of 11

Pe =

2g
Rij

[8]

where g is the surface tension [M L−1 T−1] and Rij is the radius of
the inscribed circle in the pore throat (Fig. 3c). During imbibition,
the capillary pressure must be lower than the entry pressure before
water may invade an air-saturated pore unit (Zhang et al., 2013). At
each pressure step, the equilibrium positions of the air–water interfaces within the network are determined. From this information
we can calculate the average saturation of the network at any given
capillary pressure (e.g., Øren et al., 1998; Raoof and Hassanizadeh,
2012; Raoof et al., 2013).
To obtain the SWRC, the simulation procedure was as follows.
After selecting a capillary pressure, we determined if a water-saturated pore unit could be invaded by a neighboring air-saturated
pore unit. Both pore units had to be connected to their corresponding reservoir for invasion to occur. Disconnected water-saturated
pores were used to calculate the irreducible (i.e., residual) water
saturations. An algorithm was further implemented to keep trackof disconnected air and water regions to ensure that no displacement of the disconnected pores could occur. See Sweijen et al.
(2016) for more details.

66Effect of Compaction and

Mixing on Pore Structures

Pore Body Radius

We investigated changes in the pore body radius distribution due
to the compaction and sand mixing processes. Figure 4 showschanges in the cumulative frequency distributions of the pore bodyradii for Sands S2 and S5 and their mixtures. Results for the othersands and mixtures are shown in Supplemental Fig. S1. The plots
indicate that particle mixing significantly affects the pore-size distribution while the effects of compaction are far less substantial.
Compaction caused a decrease in the porosity and, compared withmixing, caused a more uniform change in the pore space (Fig. 4aand 4b). This is consistent with previous experimental and modeling studies on the effects of compaction on pore structure (e.g.,Matthews et al., 2010; Koch et al., 2011; Kuncoro et al., 2014).
Moreover, an increase in the degree of compaction reduced thepore body radius and consequently the volume of pores. This
causes a reduction in the volume of macropores as a fraction of
the total volume of the soil medium (Koch et al., 2011; Kuncoro
et al., 2014). Mixing Sand S2 with larger sand grains (i.e., Sands
S3, S4, and S5) created mixtures with much larger pore body radii(Fig. 4c; Supplemental Fig. S1).

Fig. 4. Cumulative pore body radii distributions as affected by soil
compaction (different porosities) of sands (a) S2 and (b) S5 and (c)as obtained by mixing different proportions of sands S2 and S5 for a
porosity f of 0.40.

Pore Throat Radiusfrequency distributions of the pore throat radius obtained with thecompaction and particle mixing simulations, again for Sands S2and S5 and several of their mixtures. These results, as well as thosefor the other sands shown in Supplemental Fig. S2, show that thepore throat radii increase with increasing grain sizes, as expected.
The simulations revealed a relatively uniform change in the porethroat radii of the sands during compaction, whereas the shapesof the pore throat radius distribution curves changed far less uniformly for the mixed samples. Because SWRCs are dominated bypore throat radii, at least in the wet range, one may expect uniformchanges in the SWRC of each sand during the compaction process
but more non-uniform changes as a consequence of mixing.

Fluid flow and capillary displacement processes in porous media
are strongly affected by pore throat sizes because they provide constrictions along the flow paths. Figure 5 depicts the cumulative

Mixing our main sand, S2, with larger grains created distributions with larger pore throats (Fig. 5c). Linear relationships were

VZJ | Advancing Critical Zone Science

p. 5 of 11

Fig. 6. Effect of the mean diameter (D50, expressed as the soil particlesize) of sands S1 to S5 and their mixtures on the average pore body andpore throat radii; AR is the average radius. For all mixtures, porosity
was assumed to be 0.40.

diameters of the sand mixtures and then estimating the slopes ofthe fitted lines, such as for D50 in Fig. 6. We did the same for othersoil particle D values (i.e., for D10, D30, and D60). Our calculations
showed that the slopes of the fitted lines were much higher for D10
than D50 and D60 (results not further shown here). This indicatesthat small soil grains are far more effective in changing the average pore body size, and hence that the average pore size is affected
mostly by changes in the D10 value of a soil. In other words, theradii of relatively large pore bodies and pore throats are controlled
mostly by the smaller soil grains.

Fig. 5. Cumulative pore throat radii distributions as affected by soil
compaction (different porosities) of sands (a) S2 and (b) S5 and (c)as obtained by mixing different proportions of sands S2 and S5 for a
porosity f of 0.40.

obtained between the average pore throat and pore body radii
during both compaction and mixing (Supplemental Fig. S3).

Effect of Grain Diameter on Pore Body and
Body Throat Radii
To investigate the effect of particle diameter on pore size during
mixing of the sands, the average pore body and pore throat radiiof all sands were calculated and plotted vs. D50 values (Fig. 6). The
overall linear correlation indicates that the average pore body andthroat size can be approximately estimated at a given porosity value.
These types of correlations can be used to advantage to generatepore network models of soils with known porosity values.
We also investigated the effect of particle diameter on the pore
body and pore throat radii. This was done by plotting the averagepore body and pore throat radii as a function of the various particle

VZJ | Advancing Critical Zone ScienceWe similarly plotted the average pore body and pore throat radiivs. D 50 of the different mixtures of S2 with the other sands(Supplemental Fig. S4). A comparison of four graphs showed thatthe relationship between the average pore body and/or throat
radius and D50 will be linear if the difference between the mean
grain sizes of two main sands is very small (e.g., for S1 and S2).
However, as shown in Supplemental Fig. S4, the relationships
became nonlinear when differences between the mean grain
sizes of the sands were more substantial, such as when mixing
S2 and S5.

66Evolution of Porosity andthe Intrinsic Permeability

We next present DEM simulation results showing the effects ofcompaction and particle mixing on the pore body and pore throat
radii and the intrinsic permeability. Figure 7 depicts changes in the
average pore body and pore throat radii as a function of porosity.
The average pore radii were found to increase linearly with porosity.
Changes in porosity had more effect on the radii for the coarser
sands, as indicated by the slightly steeper slopes in Fig. 7 especiallyfor Sands S4 and S5. The plots in Fig. 7 demonstrate that different
sand samples can have different average pore sizes at a given valueof porosity.
p. 6 of 11

Fig. 7. Effect of porosity on the average pore body radius (APBR) andthe average pore throat radius (APTR) of the different sand graincompactions (expressed as the porosity f) of sands S1 to S5. The R2
values of the regression were all between 0.99 and 1.00.

Results for the intrinsic permeability, k, are shown in Fig. 8. Soil
permeability is determined by both the bulk geometric properties of a sample (e.g., porosity and the pore-size distribution) andtopological properties such as pore space connectivity and tortuosity. The effect of the pore structure on the permeability wasstudied by mixing different fractions of the main sand (S2) with
two other sand samples (S1 and S3). The results in Fig. 8a showthat mixing the fine sand (S1) with S2 caused a steady reductionin the permeability, while mixing the coarser sand (S3) with S2
increased the permeability. This confirms that the permeability
is mostly strongly affected by the larger pores (Håkansson and
Lipiec, 2000; Assouline, 2006b). Mixing S2 with S3 hence causes
an increase in the fraction of smaller pores of the coarser sand (hereS3) and consequently a reduction in the permeability of Sand S3.
By contrast, the permeability of S2 increased due to an increasein the fraction of larger pores. Our results clearly show that thepermeability reductions of the coarser sand are dominated strongly
by the relative amounts of the finer sands.
Table 2 compares calculated values of the intrinsic permeability
with measured experimental values for Sands S1 and S2, as well
as for their mixtures. While the measurements and simulation
results both showed permeability reductions toward the finersand, the simulated values for intrinsic permeability were muchhigher than the experimental values. The differences may be dueto the fact that the pore finite volume method used in our study
solves the flow problem locally at a very small scale and not at the

VZJ | Advancing Critical Zone ScienceFig. 8. Effect of sand grain mixing and compaction on the intrinsic permeability, k: (a) effect of mixing on k for three sands (S1, S2,
and S3) and their mixtures for a porosity f of 0.40; (b) relationshipbetween pore body and throat radii and k for different sand grain mixtures; (c) effect of compaction (porosity) on k for three sands (S1, S2,
and S3); and (d) values of k vs. pore body and throat radii for different
porosities (porosities of the sands ranged from 0.3 on the left side ofeach curve to 0.4 on the right sides); AR is the average radius.

p. 7 of 11

macroscopic level. This has been noted also by Tong et al. (2012),
who, in a study of the effects of mixing on the permeability ofglass beads, found similar differences between their experimental
data and DEM simulations. Another possible reason may be theassumption in the DEM that the particles have a perfect spherical shape. Clearly, more studies on formulations of the DEM areneeded to resolve the discrepancy between calculated and measured permeabilities.
Figure 8b shows that the relation between permeability and poresize became slightly nonlinear for the various mixtures of S1 andS2, while a more linear relationship was obtained for the mixturesof S2 and S3. We attribute this to the fact that Sand S1 had amore bimodal grain size distribution (Fig. 1), while S3 had a very
regular (unimodal) grain size distribution (Fig. 1). The relationships between permeability and pore size were even more nonlinear
when the permeability for a given porosity was plotted as a function of either the pore body radius or pore throat radius, as shownby the results in Fig. 8d.

and pore throat size on the SWRC. Figure 9 shows the evolutionof the volumetric water content, q w, as a function of capillarypressure, Pc , for some of the compaction and mixing simulations. Figures 9a and 9b indicate how saturated water contents
decreased by diminishing the porosity due to the compaction.
Also, however, at a given water content away from saturation, thecapillary pressure was always higher for the more compacted samples due to the creation of smaller pore throats, which regulate
higher capillary entry pressures and control capillary displacement (Supplemental Fig. S5). These findings are consistent withprevious studies by Richard et al. (2001a), Assouline (2006a),and Matthews et al. (2010), among others, showing the effect ofsoil compaction on the SWRC.

We also evaluated the effects of compaction on calculated permeability values. For this, we used three sands (S1, S2, and S3) andcalculated the intrinsic permeability at different porosity values(i.e., f = 0.3, 0.325, 0.35, 0.375, and 0.4). Figure 8c indicates an
exponential decrease in permeability as a function of porosity
following compaction. The results showed that the changes in permeability are consistent with the changes in pore body and throat
radii during compaction (Fig. 7). As explained above, the fractionof larger pores decreased when the samples were compacted. Figure
7 clearly shows that the pore body and throat radii, and hence thepore volumes, diminished with increased compaction. This may
explain the reductions in permeability with compaction. Results
also indicate that, in terms of absolute values, the permeability
during compaction changed more for the coarser materials thatthe fine-textured sands (Fig. 8c).

66Soil Water Retention Curve

Finally, we present DEM simulation results showing the effectsof compaction and particle mixing on the SWRC. This curve
has a first-order effect on predictions of water flow in variably
saturated media. We investigated the effect of both pore body
Table 2. Comparison between experimental data and simulation results
for the intrinsic permeability.
Intrinsic permeability
Sand

Simulation

Experiment———————— ´ 10−11 m2 ————————
100% S1

24.5

6.4

40% S1, 60% S2

39.3

7.9

20% S1, 80% S2

50.8

8.6

100% S2

62.6

9.1

VZJ | Advancing Critical Zone ScienceFig. 9. Effect of compaction on the water retention curves of sands (a)S2 and (b) S5 and (c) of mixing different proportions of S2 and S5 on
the soil water retention curve for a porosity f of 0.40.

p. 8 of 11

The simulated soil water retention curves were analyzed in termsof the water retention equation of van Genuchten (1980):
q w ( pc ) = q rw +q sw -q rw

é1 + a p n ù 1-1/ n
ê (
c ) ú
ë
û

pc < 0[9]where q w (dimensionless) is the volumetric water content, q r w(dimensionless) is the residual water content, q sw (dimensionless) is
the saturated water content, a [M−1 L T2] and n (dimensionless)are shape parameters, and pc [M L−1 T−2] is the capillary pressure.
Equation [9] was fitted to the SWRC data using the RETC code(van Genuchten et al., 1991). Results show that a decrease in porosity due to compaction leads to lower values of a (Fig. 10a). Figure
10b shows that an essentially linear relationship between the van
Genuchten parameter a and the average pore throat radius (APTR)
was obtained. This is not surprising because the pore throat size
has a controlling effect on the entry pressure of pores. However, as
opposed to a clear effect of compaction on the average pore throat
radius, and hence on a, we did not find any noticeable effect of compaction on the van Genuchten parameter n. Supplemental Fig. S6shows that n had little or no correlation with porosity (Supplemental
Fig. S6). This should be expected because compaction, as modeledusing the DEM, does not materially change the particle-size distribution, while the relative width (or normalized variance) of thepore-size distribution will remain very similar except for a shift in
the distribution toward smaller pore sizes (and hence smaller valuesof the van Genuchten parameter a), at least for non-macroporous
soils as considered in this study. These conclusions are very muchconsistent with several earlier experimental studies using repacked
soils (e.g., Richard et al., 2001b; Stange and Horn, 2005).
For the mixing simulations, combining Sand S2 with a coarser
sand caused a decrease in the capillary pressure at a given volumetric water content, and hence larger values of the a parameter. As
shown in Fig. 10c, the van Genuchten parameter a increased nonlinearly with D50 expressed as the sand grain size. This was again
expected because coarser soils are known to have larger a values(e.g., Carsel and Parrish, 1988). Similar to the compaction simulations, however, no clear correlation between the van Genuchten
parameter n and sand grain size was found with the sand mixing
simulations (Supplemental Fig. S6).
Finally, Fig. 10d shows that the SWRCs of all combinations of
Sand S2 with the other sands changed non-uniformly with the
APTR, except possibly for the mixture of S2 and S3. Moreover,mixing different fractions of S1 with S2 had only a relatively small
effect on the SWRCs. This was due to minor changes in the porethroat radii of the mixtures compared with the main sand (S2).

66Conclusions

In this study, we used five different sands (S1–S5) to investigate
the effects of compaction as well as mixing of different sand types

VZJ | Advancing Critical Zone ScienceFig. 10. Effect of different degrees of compaction and mixing of sands
S1 to S5 in different combinations on the van Genuchten parametera. Calculations for the mixing scenarios assumed a porosity f of 0.40.
The APTR is the average pore throat radius and Dia is the characteristic grain diameter (effective diameter D10 or mean diameter D50).

p. 9 of 11

on pore structure, intrinsic permeability, and the SWRC using
coupled discrete element and pore-scale network modeling. We
found that the average pore body and pore throat radii increased
linearly with the median particle size (defined as D50). Results forthe various sand compaction and mixing scenarios showed that theeffects of sand mixing on the pore-size distribution was far more
important than compaction as such.
The results revealed that mixing with finer sands caused a reduction in the intrinsic permeability, while mixing with coarser sandincreased the permeability of the mixture. We found that therelationship between permeability and pore size is affected by thetype of grain size distribution—unimodal or bimodal. Mixing theunimodal sand, S2, with the bimodal sand, S1, caused a nonlinear
relationship between permeability and pore body and throat radii,
while mixing with other sands having a unimodal grain-size distribution produced linear relationships between permeability andpore body and throat radii.
An exponential increase in permeability with porosity was found,with permeability increasing much faster with porosity for thecoarser sands. Our results showed that permeability increased
nonlinearly with both pore body and throat radii. As expected,
the saturated water content decreased with diminishing porosity
under compaction. As opposed to the mixing simulations, compaction caused the SWRCs to change uniformly with porosity. The
decrease in porosity, and hence the increase in bulk density, due tocompaction caused smaller values of the van Genuchten parametera. Values of a further decreased nonlinearly with decreasing sandD50 expressed as the sand grain size. No clear relationship between
the van Genuchten parameter n and sand grain size was found.
The results of this study help to better understand the effects ofsand mixing and compaction on soil pore structure and macroscopic soil hydraulic properties such as porosity, permeability, andsoil water retention processes. This provides confidence that thepore-scale approach used here can be used ultimately to develop
correlation relationships for use in field-scale models that accountfor soil compaction and mixing.

Acknowledgments

We would like to thank Filcom Company (Papendrecht, the Netherlands; www.filcom.nl)
for supplying the sands and Hamed Aslannejad for taking microscopic images of sands. Comments by two anonymous referees and the VZJ Associate Editor further helped to improvethe manuscript. This work is part of the Veni Research Program with Project no. 016.151.047,
which is financed by the Netherlands Organization for Scientific Research (NWO).

References

Assouline, S. 2006a. 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. 2006b. Modeling the relationship between soil bulk density
and the hydraulic conductivity function. Vadose Zone J. 5:697–705.
doi:10.2136/vzj2005.0084
Baumhardt, R.I., 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

VZJ | Advancing Critical Zone ScienceBelheine, N., J.P. Plassiard, F.V. Donzé, F. Darve, and A. Seridi. 2009.
Numerical simulation of drained triaxial test using 3D discrete element
modeling. Comput. Geosci. 36:320–331.
Berli, M., F. Casini, W. Attinger, R. Schulin, S.M. Springman, and J.M. Kirby.
2015. Compressibility of undisturbed silt loam soil: Measurements andsimulations. Vadose Zone J. 14(8). doi:10.2136/vzj2014.10.0153
Bradford, J.M., J.E. Ferris, and P.A. Remley. 1987. Interrill soil erosion
processes: I. Effect of surface sealing on infiltration, runoff, andsoil splash detachment. Soil Sci. Soc. Am. J. 51:1566–1571.
doi:10.2136/sssaj1987.03615995005100060029x
Carsel, R.F., and R.S. Parrish. 1988. Developing joint probability distributionsof soil water retention characteristics. Water Resour. Res. 24:755–769.
doi:10.1029/WR024i005p00755
Chareyre, B., A. Cortis, E. Catalano, and E. Barthélemy. 2012. Pore-scalemodeling of viscous flow and induced forces in dense sphere packings.
Transp. Porous Media 94:595–615. doi:10.1007/s11242-012-0057-2
Cnudde, V., and M.N. Boone. 2013. High-resolution X-ray computed
tomography in geosciences: A review of the current technology andapplications. Earth Sci. Rev. 123:1–17. doi:10.1016/j.earscirev.2013.04.003
Coetzee, C.J., and D.N.J. Els. 2009. Calibration of discrete element
parameters and the modelling of silo discharge and bucket filling.
Comput. Electron. Agric. 65:198–212. doi:10.1016/j.compag.2008.10.002
Cui, Y., J. Wooster, P. Baker, S. Dusterhoff, L. Sklar, and W. Dietrich. 2008.
Theory of fine sediment infiltration into immobile gravel bed. J. Hydraul.
Eng. 134:1421–1429. doi:10.1061/(ASCE)0733-9429(2008)134:10(1421)
Cundall, P.A., and O.D.L. Strack. 1979. A discrete numerical model for granularassemblies. Geotechnique 29:47–65. doi:10.1680/geot.1979.29.1.47
Evans, E., and A.C. Wilcox. 2014. Fine sediment infiltration dynamics in agravel-bed river following a sediment pulse. River Res. Appl. 30:372–
384. doi:10.1002/rra.2647
Gee, G., and D. Or. 2002. Particle size analysis. In: J.H. Dane and G.C. Topp,
editors, Methods of soil analysis. Part 4. Physical methods. SSSA Book
Ser. 5. SSSA, Madison, WI. p. 255–294. doi:10.2136/sssabookser5.4.c12
Glab, T., and S. Kopec. 2009. Effect of soil compaction on root systemmorphology and yields of meadow fescue (Festuca pratensis). Pol. J.
Environ. Stud. 18:219–225.
Gladkikh, M., and S. Bryant. 2005. Prediction of imbibition in
unconsolidated granular materials. J. Colloid Interface Sci. 288:526–
539. doi:10.1016/j.jcis.2005.03.029
Gregorich, E.G., N.B. McLaughlin, D.R. Lapen, and P. Rochette. 2014.
Soil compaction, both an environmental and agronomic culprit:
Increased nitrous oxide emissions and reduced plant nitrogen uptake.
Soil Sci. Soc. Am. J. 78:1913–1923. doi:10.2136/sssaj2014.03.0117
Hadas, A. 2004. Seedbed preparation: The soil physical environment ofgermination seeds. In: R.L. Benech-Arnold and R.A. Sanchez, editors,
Handbook of seed physiology: Applications to agriculture. Hayworth
Press, New York. p. 3–49.
Haghnazari, F., H. Shahgholi, and M. Feizi. 2015. Review: Factors affecting
the infiltration of agricultural soils. Int. J. Agron. Agric. Res. 6(5):21–35.
Håkansson, I., and J. Lipiec. 2000. A review of the usefulness of relative
bulk density values in studies of soil structure and compaction. Soil
Tillage Res. 53:71–85. doi:10.1016/S0167-1987(99)00095-1
Keller, T., M. Lamande, S. Peth, M. Berli, J.-Y. Delenne, W. Baumgarten, et al.
2013. An interdisciplinary approach towards improved understandingof soil deformation during compaction. Soil Tillage Res. 128:61–80.
doi:10.1016/j.still.2012.10.004
Koch, K., A. Kemna, J. Irving, and K. Holliger. 2011. Impact of changes ingrain size and pore space on the hydraulic conductivity and spectral
induced polarization response of sand. Hydrol. Earth Syst. Sci. 15:1785–
1794. doi:10.5194/hess-15-1785-2011
Kuncoro, P.H., K. Koga, N. Satta, and Y. Muto. 2014. A study on the effectof compaction on transport properties of soil gas and water: I. Relativegas diffusivity, air permeability, and saturated hydraulic conductivity.
Soil Tillage Res. 143:172–179. doi:10.1016/j.still.2014.02.006
Li, S., A. Raoof, and R. Schotting. 2014. Solute dispersion under electricand pressure driven flows: Pore scale processes. J. Hydrol. 517:1107–
1113. doi:10.1016/j.jhydrol.2014.06.049
Mahmoodlu, M.G., S.M. Hassanizadeh, N. Hartog, and A. Raoof. 2014.
Oxidation of trichloroethylene, toluene, and ethanol vapors by apartially saturated permeable reactive barrier. J. Contam. Hydrol.
164:193–208. doi:10.1016/j.jconhyd.2014.05.013
Mahmoodlu, M.G., S.M. Hassanizadeh, N. Hartog, A. Raoof, and M.Th.
van Genuchten. 2015. Evaluation of a horizontal permeable reactive
barrier for preventing upward diffusion of volatile organic compounds

p. 10 of 11through the unsaturated zone. J. Environ. Manage. 163:204–213.
doi:10.1016/j.jenvman.2015.08.025
Masís-Meléndez, F., L.W. de Jonge, T.K.K. Chamindu Deepagoda, M.
Tuller, and P. Moldrup. 2015. Effects of soil bulk density on gas transportparameters and pore-network properties across a sandy field site.
Vadose Zone J. 14(7). doi:10.2136/vzj2014.09.0128
Matthews, G.P., G.M. Laudone, A.S. Gregory, N.R.A. Bird, A.G. Matthews,
and W.R. Whalley. 2010. Measurement and simulation of the effectof compaction on the pore structure and saturated hydraulicconductivity of grassland and arable soil. Water Resour. Res.
46:W05501. doi:10.1029/2009WR007720
Mehlhorn, T., S. Prohaska, U. Homberg, and V. Slowik. 2008. Modelling
and analysis of particle and pore structures in soils. In: K.J. Witt, editor,
Workshop Internal Erosion. 26–28 Nov. 2008. Schriftenreihe Geotechnik
21. Bauhaus-Universität, Weimar, Germany. p. 53–60.
Mousavi, M.A., and S.L. Bryant. 2012. Connectivity of pore space as acontrol on two-phase flow properties of tight-gas sandstones. Transp.
Porous Media 94:537–554. doi:10.1007/s11242-012-0017-x

Sela, S., T. Svoray, and S. Assouline. 2015. The effect of soil surface sealingon vegetation water uptake along a dry climatic gradient. Water
Resour. Res. 51:7452–7466. doi:10.1002/2015WR017109
Shmulevich, I., Z. Asaf, and D. Rubinstein. 2007. Interaction between soil
and a wide cutting blade using the discrete element method. Soil
Tillage Res. 97:37–50. doi:10.1016/j.still.2007.08.009
Šmilauer, V., E. Catalano, B. Chareyre, S. Dorofeenko, J. Duriez, A. Gladky,
et al. 2015. Yade documentation. 2nd ed. The Yade Project, Grenoble,
France. doi:10.5281/zenodo.34073
Stange, C.F., and R. Horn. 2005. Modeling the soil water retention
curve for conditions of variable porosity. Vadose Zone J. 4:602–613.
doi:10.2136/vzj2004.0150
Sweijen, T., E. Nikooee, S.M. Hassanizadeh, and B. Chareyre. 2016. The
effects of swelling and porosity change on capillarity: DEM coupledwith a pore unit assembly method. Transp. Porous Media 113:207–226.
doi:10.1007/s11242-016-0689-8

Øren, P., S. Bakke, and O.J. Arntzen. 1998. Extending predictive
capabilities to network models. Soc. Pet. Eng. J. 3:324–336.

Tong, A.T., E. Catalano, and B. Chareyre. 2012. Pore-scale flow
simulations: Model predictions compared with experiments on bidispersed granular assemblies. Oil Gas Sci. Technol. 67:743–752.
doi:10.2516/ogst/2012032

Peth, S., J. Nellesen, G. Fischer, and R. Horn. 2010. Non-invasive 3D
analysis of local soil deformation under mechanical and hydraulicstresses by mCT and digital image correlation. Soil Tillage Res. 111:3–
18. doi:10.1016/j.still.2010.02.007

Torskaya, T., V. Shabro, C. Torres-Verdín, R. Salazar-Tio, and A. Revil. 2014.
Grain shape effects on permeability, formation factor, and capillarypressure from pore-scale modeling. Transp. Porous Media 102:71–90.
doi:10.1007/s11242-013-0262-7

Prodanović, M., A. Mehmani, and A.P. Sheppard. 2015. Imaged-based
multiscale network modelling of microporosity in carbonates. Geol.
Soc. London Spec. Publ. 406:95–113. doi:10.1144/SP406.9

Unsal, E., G. Mason, N.R. Morrow, and D.W. Ruth. 2009. Bubble
snap-off and capillary-back pressure during counter-current
spontaneous imbibition into model pores. Langmuir 25:3387–3395.
doi:10.1021/la803568a

Rajkai, K., S. Kabos, M.Th. van Genuchten, and P.E. Jansson. 1996.
Estimation of water retention characteristics from bulk densityand particle distribution of Swedish soils. Soil Sci. 161:832–845.
doi:10.1097/00010694-199612000-00003

Vahyala I.E., A.I. Zata, G.A. Adediron, I.S. Dzarma. 2013. Impacts ofcompaction on selected physical properties of tropical soils at Yola,
Nigeria. Adv. Agric. Res. Dev. 1(1):9–15.

Raoof, A., and S.M. Hassanizadeh. 2010. A new method for generating
pore-network models of porous media. Transp. Porous Media 81:391–
407. doi:10.1007/s11242-009-9412-3

van Genuchten, M.Th. 1980. A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–
898. doi:10.2136/sssaj1980.03615995004400050002x

Raoof, A., and S.M. Hassanizadeh. 2012. A new formulation for porenetwork modeling of two-phase flow. Water Resour. Res. 48:W01514.
doi:10.1029/2010WR010180

van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC code forquantifying the hydraulic functions of unsaturated soils. EPA/600/291/065. USEPA, R.S. Kerr Environ. Res. Lab., Ada, OK.

Raoof, A., S.M. Hassanizadeh, and A. Leijnse. 2010. Upscaling transport ofadsorbing solutes in porous media: Pore-network modeling. Vadose
Zone J. 9:624–636. doi:10.2136/vzj2010.0026

Vasilyev, L., A. Raoof, and J.M. Nordbotten. 2012. Effect of mean network
coordination number on dispersivity characteristics. Transp. PorousMedia 95:447–463. doi:10.1007/s11242-012-0054-5

Raoof, A., H.M. Nick, S.M. Hassanizadeh, and C.J. Spiers. 2013. PoreFlow:
A complex pore-network model for simulation of reactive transportin variably saturated porous media. Comput. Geosci. 61:160–174.
doi:10.1016/j.cageo.2013.08.005

Vaz, C.M.P., M. Tuller, P.R.O. Lasso, and S. Crestana. 2014. New
perspectives for the application of high-resolution benchtop
X-ray microCT for quantifying void, solid and liquid phases insoils. 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. p. 261–281.
doi:10.1007/978-3-319-06013-2_12

Rawls, W.J., D.L. Brakensiek, and B. Soni. 1983. Agricultural management
effects on soil water processes: I. Soil water retention and Green
and Ampt infiltration parameters. Trans. ASAE 26:1747–1752.
doi:10.13031/2013.33837
Richard, G., I. Cousin, J.F. Sillon, A. Bruand, and J. Guérif. 2001a.
Effect of compaction on the porosity of a silty soil: Influence
on unsaturated hydraulic properties. Eur. J. Soil Sci. 52:49–58.
doi:10.1046/j.1365-2389.2001.00357.x
Richard, G., J.F. Sillon, and O. Marloie. 2001b. Comparison of inverse anddirect evaporation methods for estimating soil hydraulic properties
under different tillage practices. Soil Sci. Soc. Am. J. 65:215–224.
doi:10.2136/sssaj2001.651215x
Scholtès, L., and F.V. Donzé. 2012. Modelling progressive failure offractured rock masses using a 3D discrete element method. Int. J.
Rock Mech. Min. Sci. 52:18–30. doi:10.1016/j.ijrmms.2012.02.009
Scholtès, L., P.Y. Hicher, F. Nicot, B. Chareyre, and F. Darve. 2009. On thecapillary stress tensor in wet granular materials. Int. J. Numer. Anal.
Methods Geomech. 33:1289–1313. doi:10.1002/nag.767

VZJ | Advancing Critical Zone ScienceWall, A., and J. Heiskanen. 2009. Soil-water content and air-filled porosity
affect height growth of Scots pine in afforested arable land in Finland.
For. Ecol. Manage. 257:1751–1756. doi:10.1016/j.foreco.2009.01.029
Xu, Y., K.D. Kafui, C. Thornton, and G. Lian. 2002. Effects of material
properties on granular flow in a silo using DEM simulation. Part. Sci.
Technol. 20:109–124. doi:10.1080/02726350215338
Yuan, C., B. Chareyre, and F. Darve. 2015. Pore-scale simulationsof drainage in granular materials: Finite size effects and therepresentative elementary volume. Adv. Water Resour. (in press).
doi:10.1016/j.advwatres.2015.11.018
Zhang, Q., S.M. Hassanizadeh, N.K. Karadimitriou, A. Raoof, B. Liu, P.J.
Kleingeld, and A. Imhof. 2013. Retention and remobilization of colloids
during steady-state and transient two-phase flow. Water Resour. Res.
49:8005–8016. doi:10.1002/2013WR014345
Zhang, Q., A. Raoof, and S.M. Hassanizadeh. 2015. Pore-scale study offlow rate on colloid attachment and remobilization in a saturated
micromodel. J. Environ. Qual. 44:1376–1383. doi:10.2134/jeq2015.01.0058

p. 11 of 11