Relating Hydraulic Conductivity to Particle Size Using DEM

It is widely known that the packing characteristics of a granular material are governed by the size of the particles. It has been accepted for many years for example that the porosity (or density) of a granular soil depends on the distribution and range of particle sizes (e.g., Burmister 1938). The packing of a granular material strongly influences the macroscopic engineering properties; for example, in the absence of crushing, the density of a sand determines the dilatancy and peak strength when sheared (e.g., Bolton 1986), while the size and connectivity of the pores control the ease with which fluid can flow through the soil, that is, the hydraulic conductivity (e.g., Mitchell and Soga 2005).

There are many factors controlling the hydraulic conductivity of a soil–fluid system, which include both soil and fluid properties. Many predictive methods for example take into account the fluid viscosity and density, the degree of saturation, temperature, and the shape and tortuosity of the pores. For fully saturated conditions with a given fluid however, predictions of the hydraulic conductivity are based primarily on the particle size(s), which are assumed to control the size of the pores. Indeed, it has been accepted for over a century that the hydraulic conductivity of soil is largely governed by some characteristic smallest particle or pore size (e.g., Hazen 1892), although the number of published predictive equations reveal that this remains unsolved (e.g., Chapuis 2012). The focus here therefore is specifically on the governing role of pore size on the saturated hydraulic conductivity.

Recent work (de Bono and McDowell 2018b) has sought to establish what controls the porosity (the nonsolid fraction of volume) of granular materials with arbitrary particle size distributions and particle shapes. It was shown that for a confined granular material, there exists a group of “smallest particles” that govern the pore space and therefore the overall porosity or density. In this paper, it will be established that the same particles can be said to control the hydraulic conductivity, and this will be shown to be consistent with existing classical predictive methods.

Predicting the permeability or hydraulic conductivity of a granular (or any porous) material from known or estimated physical properties is a long-standing problem. The flow of groundwater or any other fluid through soils and other porous media is of major interest to geotechnical, geological, and petroleum engineers. The general established formula for relating hydraulic conductivity to particle (or pore) size is of the formwhere the conductivity k is in units of velocity and d is a characteristic particle or pore size in units of length. The coefficient c is typically empirical and may or may not include properties such as porosity or fluid viscosity, and so forth.

The conductivity k is used in Darcy's law for steady-state seepage through porous media according to v = ki = −k(dh/dl), where v is the superficial velocity, i is the hydraulic gradient, h is the potential, and l is the distance over which the potential drops. Probably the oldest and most well-known predictive method for hydraulic conductivity is that of Hazen (1892), which defines k aswhere d₁₀ is the particle diameter that 10% in terms of mass of the soil is finer than, measured from a graphical semi-logarithmic particle size distribution. Hazen's formula clearly states that the hydraulic conductivity of soil is governed by the d₁₀ size; most subsequent predictive methods for calculating k likewise use d₁₀ or similar (such as d₅). For an excellent overview of such methods, see Chapuis (2012), in which Table 3 provides a useful summary.

The usual starting point when deriving a prediction for hydraulic conductivity is steady laminar flow in a circular pipe, which, according to the Poiseuille formula, gives the volume flow rate Q (volume of fluid passing per unit time) as a function of the pipe diameter d and drop in piezometric pressure Δp as Q = −(πd⁴/128μl)Δp, where μ is the viscosity and l is the length over which the piezometric pressure drops. Dividing piezometric pressure by the unit weight of water γ_w gives the potential, and recasting this gives the velocity of the fluid as v_w = −(d²γ_w/32μ)(Δh/l). If this equation is to apply to an idealized set of pipes in soil through the voids, then one must account for the fact that v_w applies to fluid velocity and not superficial velocity, and that voids are tortuous. This can be simplified by writing the superficial velocity as v ∝ (d²γ_w/μ)(Δh/l), which is Darcy's law if the hydraulic conductivity is

To establish the permeability as a function of particle size, Kozeny and Carman (Carman 1956; Mitchell and Soga 2005) considered the pore space in a porous material as a series of pipes with tortuous paths. Taking into account the effect of surface friction on flow, they used the hydraulic radius, r_h,to derive a formula for the hydraulic conductivity, which for a given porosity, soil material, and fluid gave the conductivity as a function of r_hwhich is consistent with Eq. (3) and is therefore similar to the Hazen formula in using a characteristic linear size (but also dependent on the porosity). For a circular pipe, r_h = r/2 = d/4, with 4r_h representing the actual pipe diameter. For a cylindrical pipe with a uniform cross section, the hydraulic radius is given by

For a granular material with a uniform porosity n, the pore space can be considered as consisting of a series of pipes or as a single pipe with a highly complicated cross section. Using Eq. (6), the hydraulic radius can therefore be expressed aswhere V_V = volume of voids (i.e., the pore space) and S_T = total surface area of the solid particles. This can be rewritten in terms of the porosity n and void ratio e aswhere S = surface area per unit volume of the bulk material and S₀ = specific surface area of the particles (the surface area per unit volume of solid), and the void ratio e is the ratio of the volume of voids to the volume of solids. The well-known form of the Kozeny–Carman equation used in geotechnical engineering uses (e/S₀) instead of r_h.

The Hazen and Kozeny–Carman methods are the two most well-known methods for predicting soil hydraulic conductivity, with the latter regarded as being more accurate. The values of r_h and d₁₀ used in the above methods will now be reconciled with the recent work on the smallest particle size, which has been shown to control the overall pore space and density (de Bono and McDowell 2018b).

Recent work based on fractal theory and statistics of soil particle strengths has quantified the rate at which the volume of soil reduces when subjected to monotonic increasing stress (de Bono and McDowell 2018a; McDowell 2005). This was based on the fact that the porosity or void ratio of a granular soil must be controlled by the smallest particle size, d_sm. The total volume of voids (pores) was considered to be proportional to the cumulative volume of the smallest particles. The theory was shown to be valid for the compression of real sands (de Bono and McDowell 2018a), for which fractal size distributions emerge and self-similarity enables the changing size of the smallest particles to be estimated. However, in subsequent simulations investigating the packing of spheres, it proved difficult to quantify the smallest particles in an expanding fractal distribution of sizes due to there always only being one or two particles with the absolute smallest dimension out of many thousands or millions. It was then suggested (de Bono and McDowell 2018b) that particles that are mainly surrounded by larger particles should be considered the smallest and therefore representative of the surrounding pore space. Specifically, those particles for which ≥50% of the contacts are with larger neighboring particles are the smallest particles, of size d_sm. This novel approach took into account both the particle size distribution and the geometrical packing. The use of interparticle contacts to classify the smallest particles also meant it was possible to account for local variations in the particle size distribution. Using discrete element method (DEM) simulations, it was shown that the cumulative volume, V_sm, of these particles demonstrated direct proportionality with the volume of voids, V_V, over a range of particle shapes, independent of stress path and particle size distribution (including nonfractal distributions). This was also confirmed to be the case in idealized Apollonian sphere packings (de Bono and McDowell 2018b), for which it has been shown analytically that the overall pore volume should also be proportional to the volume of the smallest particles (Anishchik and Medvedev 1995).

Considering that these smallest particles have been shown to govern the void space, and that it is accepted that hydraulic conductivity is related to the squared value of some characteristic pore diameter (e.g., Carrier 2003; Chapuis 2012), it follows that these so-called smallest particles might also determine the hydraulic conductivity. That is, the characteristic size d_sm of the smallest particles defined by de Bono and McDowell (2018b) may accurately predict the hydraulic conductivity and/or provide a physical justification of traditional parameters such as d₁₀ or r_h.

Real fractal distributions of sizes are finite, and from the definition of a fractal (e.g., Palmer and Sanderson 1991), it can be deduced that the number of particles N of a given size d (or size range) is $N(d)\sim d^{-D}$, where D is the fractal dimension. For a granular material subjected to confined crushing, the fractal dimension invariably tends to D = 2.5 (Turcotte 1986), which is the same value found in Apollonian sphere packings (Anishchik and Medvedev 1995). In assuming that the volume of voids V_V was proportional to the volume of the smallest particles, the volume of voids was given in McDowell (2005) asand for such a fractal distribution, in McDowell and Bolton (1998) the total surface area of all particles S_T was given bywhich means that the hydraulic radius parameter in Eq. (7) and used in the Kozeny–Carman theory is given by

This analysis will work for any fractal dimension D, and means that, according to the theory of McDowell (2005) and the Kozeny–Carman theory (Carman 1956), the smallest particle size d_sm is the correct characteristic size to use in predicting the permeability or hydraulic conductivity of soil, whatever that smallest particle size d_sm can be reasonably deemed to be.

It will now be examined whether there is any agreement between the value of d_sm and the value of d₁₀ that Hazen or r_h that Kozeny and Carman would have used. A simple simulation of crushable spheres was performed using the DEM and the strength parameters specified in McDowell and de Bono (2013). Details are the same as used previously (de Bono and McDowell 2018b): to summarize, the sample of spheres was subjected to increasing isotropic normal stress, and when the average octahedral shear stress within any sphere reached its strength, the sphere was replaced by two smaller spheres, obeying conservation of mass. It should be noted that although this simulation uses spheres, the fundamental concept has been shown to be valid for nonspherical particles (de Bono and McDowell 2018b). In terms of hydraulic conductivity, particle shape will affect the tortuosity and therefore also the hydraulic conductivity, but the focus here is solely on the governing effect of the pore size. The cylindrical sample initially contained approximately 6,000 spheres of d = 2 mm and contained around 155,000 spheres at the end of compression. The macroscopic results in terms of void ratio and stress are given in Fig. 1. If the smallest particles are defined as any particle for which at least 50% of neighboring particles are larger, then the cumulative volume of these particles, V_sm, is directly proportional to the void space. This is shown and explained in more detail in de Bono and McDowell (2018b), along with additional simulations using different particle shapes.

There are many different ways of calculating a representative size of these smallest particles (or the pores, which will be discussed later). A straightforward method of obtaining an effective size of these smallest particles, d_sm, is using the diameter of the volumetric average (V_sm divided by the total number of smallest particles). The values of d_sm are shown on particle size distributions in Fig. 2, showing the evolution due to crushing. It can be seen that d_sm approaches d₁₀, particularly at high stresses as the grading evolves to a fractal, lending new credibility to the Hazen approach.

The volume of voids and the total particle surface area can be calculated directly in DEM for use in Eq. (7) to calculate the hydraulic radius. The values of d_sm are plotted against the hydraulic radius r_h in Fig. 3 for the same simulation. The values are directly proportional at medium-to-high stresses, showing a clear relationship. It is generally accepted that the Kozeny–Carman approach to calculating hydraulic conductivity is one of the most accurate predictive methods (e.g., Carrier 2003), and the fact that the hydraulic radius is entirely consistent with d_sm highlights the significance of both d_sm and the method of defining the smallest particles. Carman (1956) stated that 4r_h can be regarded as some mean pipe diameter of significant hydrodynamical importance, so it is interesting to note from the gradient in Fig. 3 that 4r_h ≈ 0.57d_sm.

For granular materials, predictions of hydraulic conductivity tend to be based on particle size(s), due to the fact that the particle sizes govern the pores. For nongranular porous materials (e.g., sedimentary rock), predictions of k tend to be based directly on estimations of the pore or throat size(s). Pore-size distributions can be estimated using, for example, mercury intrusion porisometry, which shows the cumulative pore volume against (decreasing) pore diameter. It is possible to obtain analogous information about the pore network in DEM from the size distribution of the smallest particles. Given that the overall pore space is proportional to the cumulative volume of these smallest particles, it can be assumed then that each of these particles is proportional to its surrounding pore space. That is, each small particle has an associated volume of pore space, the sum of which is equal to the overall pore space. This assumption means that the volumetric pore-size distribution should be equivalent to that of the smallest particles. Fig. 4 shows the evolution of the cumulative volumetric size distribution for the smallest particles from the DEM simulation. The constant of proportionality between the total volume of voids and the volume of the smallest particles is 1.7, that is, V_V = 1.7 × V_sm (de Bono and McDowell 2018b). So for any small particle of size d, the diameter of the surrounding pore, d_p, might be assumed to be d_p = d × 1.7^(1/3); this pore-size scale is also plotted in Fig. 4.

Two well-known methods for predicting permeability using a characteristic pore size are by Winland (e.g., Nelson 2009) and Katz and Thompson (1986). The Winland method uses the pore size for which 35% of the pore space is in larger pores, as the size in Eq. (1), which, similar to the Hazen approach, uses a single size and does not include information about the distribution in its entirety. The Katz and Thompson method suggests using the pore size corresponding to the inflection point on the cumulative pore-size distribution (e.g., Fig. 4). This value, d_pi, represents a modal pore size in terms of volume and was interpreted as being the smallest throat size forming part of a continuous pore network through the material. The suggestion that d_pi represents the minimum pore size of significance implies that this is more useful than using an arbitrary fixed value from some distribution (such as the Hazen or Winland methods); however, it does not necessarily account for the range of the pore sizes. From the curves in Fig. 4, inflection points have been simply estimated, and the corresponding values of d_pi compared with the hydraulic radius and d_sm in Fig. 3. The different measures were all consistent with each other, which supports the assumptions used here in estimating the pore-size distribution.

It is not currently feasible to use the contact distribution in a real soil to identify the smallest particles, which would enable V_sm or d_sm to be correlated to the volume of voids or permeability, confirming the theories proposed here. Nor is it possible to determine the permeability of any numerical DEM samples with confidence. This limiting factor in the DEM simulations is generally the particle size: the existence of increasingly fine particles causes the time step to become very small and the simulations to become unwieldy. To simulate the flow of fluid, the mesh size would need to be several orders of magnitude smaller than the smallest particle, meaning that this is not currently computationally within our means (not to mention the inherent assumptions involved). Having said that, it remains a priority to physically validate the theory here; one option that may be possible in the future is using X-ray tomography to obtain particle and contact distributions prior to permeability testing.

Although the lack of validation is an issue, the volume of the smallest particles has been rigorously shown to correlate directly with the pore space, for any particle shape (de Bono and McDowell 2018b), and it is logical that these particles also govern the hydraulic conductivity–for which it has been accepted for more than a century is controlled by the squared value of some characteristic diameter. In addition to providing physical justification for classic measures (e.g., the hydraulic radius), it is hoped that this theory will prove useful in future applications and lead to further microscopic insight, when more advanced investigative techniques become available.

The method of categorizing the smallest particles using contacts (as described in de Bono and McDowell 2018b) has the advantage of taking into account both the particle size distribution and the geometrical packing of the bulk material, and the value of d_sm correlated very well with the hydraulic radius used in many of the most reliable permeability predictions, providing strong evidence that these are the particles of most significance.

Interestingly, changing the contact selection criterion for the smallest particles, for example from 50% (or more) to 35% or 60%, has no effect on the proportionality between the overall pore space and the volume of the smallest particles, or d_sm, which is consistent with the self-similarity expected from a fractal distribution. This observation was also found to be consistent with the analysis of Apollonian packings. Using an extreme value however causes the direct proportionality observed between the cumulative volume of small particles and the overall pore space to be lost. For example, if 100% is used—that is, only particles entirely surrounded by larger particles are deemed the smallest—then this ignores many pairs of two small particles surrounded by larger ones. Similarly, using a low value such as 10% results in the vast majority of particles being classed as the smallest, and does not correlate with the volume of voids.

At the very least, this work provides a very simple and easy-to-implement method of quantifying or discretizing the void space and deducing pore-size distributions in numerical DEM samples for which the continuous pore space is not easily dividable. From these pore-size distributions, it should be possible to identify additional or more useful micro properties, for example a more appropriate characteristic size for use in explaining permeability than simply the diameter of the volumetric average (d_sm). It may also be possible to further refine this work and relate the pore sizes to the actual throat sizes between the pores, or as part of a network model (e.g., Bryant and Blunt 1992).

This paper has sought to establish how the smallest particles in a granular material control the (saturated) hydraulic conductivity, a concept that has been accepted for a great many number of years. It has been shown here that the smallest particles—defined as those with fewer (or equal) smaller neighbors than larger ones—should be those that govern the permeability. These particles have previously been shown to govern the volume of pore space, and a simple fractal theory presented here confirms that the size of the smallest particles, d_sm, is the correct characteristic size when used to estimate the hydraulic conductivity in traditional equations. DEM results have shown a strong correlation between a simple characteristic measure of size of these particles, d_sm, and the hydraulic radius used by Kozeny–Carman as well as d₁₀ used traditionally by the geotechnical community. In addition to providing a new analytical basis for predicting hydraulic conductivity, this work lends new micro mechanical credibility to classical predictive formulas. Volumetric pore-size distributions were also deduced from the smallest particles, from which it was possible to perform further pore-size analytics.

Some or all data, models, or code generated or used during the study are available from the corresponding author by request (i.e., raw results), while others may only be provided with restrictions (i.e., DEM software and code).

This work was supported by the Engineering and Physical Sciences Research Council (Grant number EP/L019779/1).

The following symbols are used in this paper: