Unsupported browser

For a better experience, please update your browser to its latest version.

Your browser appears to have cookies disabled. For the best experience of this website, please enable cookies in your browser

We'll assume we have your consent to use cookies, for example so you won't need to log in each time you visit our site.
Learn more

Technical Paper: A particulate perspective on soil mechanics

Cracked earth

Catherine O’Sullivan, Imperial College. This paper was published in GE’s March 2016 edition.

Arup Tech Paper logo (NB)


This paper gives a summary overview of material presented in the 2015 Géotechnique Lecture.

The lecture showed how a particulate perspective can provide useful insight into soil mechanics. Use of both particulate discrete element modelling (DEM) and micro computed tomography (micro CT) were considered. Three research studies were considered to show how these techniques were used to revisit and supplement earlier published research in Géotechnique.

This synthesis of the lecture firstly outlines how DEM and micro CT are used in particle-scale soil mechanics, prior to giving a summary of each of the three research studies. The lecture can be viewed online by members of the British Geotechnical Association via the ICE website and selected slides and further details on the research is available on the author’s webpage.

Discrete element modelling

It is apt to consider DEM in a Géotechnique Lecture as the original paper documenting the DEM algorithm was published in Géotechnique in 1979 (Cundall and Strack, 1979). This is the most cited Géotechnique paper and has the most significant scientific impact among all papers that have been published in Géotechnique. The reason this paper has attracted such interest is that Cundall and Strack proposed a means to simulate granular materials at the particle scale. This technique has applications beyond soil mechanics; for example it can be applied to study particles of infant formula (powdered baby milk) and semi-solid metals.

The DEM algorithm is conceptually simple; the basic concept is that granular materials, including soil, can be idealised as rigid bodies with shapes that can be easily geometrically described. In 3D the simplest shape to consider is a sphere. The algorithm allows for overlap between the particles; this overlap is taken to represent the deformation that takes place at the contacts between real physical particles.

The forces between particles are calculated using relatively simple force:deformation models. In the contact tangential direction the shear force is capped at m times the normal force, where m is the coefficient of friction. In each time step, the DEM software will determine which particles are contacting and calculate the contact forces, calculate the resultant force acting on each particle by consideration of the contact forces, gravity forces, and boundary forces. The particle accelerations are calculated from the resultant forces according to Newton’s second law. The particle velocities and displacements are determined from numerical integration of the velocity values.

Typically an explicit, conditionally stable time-stepping approach is used and so very small time increments are considered. The small time step and the need to search to find contacts mean that DEM simulations are highly time-intensive; however, the algorithm is also well suited to exploit high performance computing.

In principle, DEM simulations are useful to provide insight into soil deformation mechanisms in situ; for example Bym et al (2013) used DEM to look at tunnel induced settlements. However, the bulk of geomechanics related research has used DEM in fundamental studies of soil behaviour.

Using DEM, a virtual laboratory can be created. DEM simulations can explore stress states that cannot easily be attained in the laboratory (as discussed further below). DEM simulations can provide insight into the particle-scale mechanics underlying the response observed in conventional laboratory tests. Detail on the number of particle-particle contacts, the magnitude of the contact forces and the local stresses can be extracted.

Figure 1: DEM simulation of direct shear test

Figure 1: DEM simulation of direct shear test

For example figure 1(a) illustrates the response observed in a simulation of a direct shear test on spherical particles, as discussed by Cui and O’Sullivan (2006). The particle displacement vectors are shown in figure 1(b). The displacement field is clearly highly non-uniform. Figure 1(c) illustrates the manner in which forces are transmitted in a soil sample in a direct shear test. The cylindrical tubes in the diagram represent the contact forces between spherical particles in the simulation. The colour and tube diameter are proportional to the contact force. It is clear that the force transmitting network of contacts is highly heterogeneous.

Subnetworks of contacts carrying above average forces form; these are called force chains and observation of these contact forces over the entire test duration indicates that during deformation there is an ongoing process of buckling of force chains followed by formation of new force chains. Considering either the displacements or the force transmission, the response in this test is markedly more complicated than the schematic diagrams of the direct shear test comprising rigid bodies moving relative to each other typically provided in undergraduate soil mechanics textbooks. DEM can therefore challenge long-held perceptions of soil behaviour.

DEM simulations are computationally expensive. Simulations with more than 100,000 particles are rare, yet there will be over 150,000 200 micro sand grains in a cubic volume 10mm x 10mm x 10mm. DEM is, however, well suited to parallel processing, or high performance computing. The DEM simulations presented here used a modified version of the LAMMPS code (Plimpton, 1995) and were run on high performance computers.

Micro computed tomography

Computed tomography was originally developed for medical applications. This technology generates 3D images of the internal structure of solid bodies by combining the information provided in a series of 2D radiographs (X-ray images).

Figure 2(a): Resin impregnation in triaxial apparatus prior to sub coring for micro CT scanning

Figure 2(a): Resin impregnation in triaxial apparatus prior to sub coring for micro CT scanning

The combination process is called reconstruction and the resultant 3D image shows the attenuation of X-rays by the material in the scanned sample; different colours (grayscales) in the image indicate different levels of X-ray attenuation. The amount of attenuation depends on the material density and atomic number and so the image can be segmented to identify different material phases by identifying a threshold colour intensities marking the transition from one material to another. Individual soil particles can be differentiated from each other using the watershed algorithms originally developed to identify catchments associated with particular rivers.

In a materials micro CT scanner suited to examine the internal structure of soil, a small, typically cylindrical, soil sample is placed between an X-ray source and a detector. The sample can be rotated in a controlled manner and radiographs are recorded at small angular intervals.

Figure 2(b): Representative micro CT image of Reigate sand

Figure 2(b): Representative micro CT image of Reigate sand

For the research discussed here samples were prepared in a conventional triaxial apparatus (figure 2(a)). Once the sample has been created and taken to the stress or deformation state of interest it is impregnated with an epoxy resin. Once the resin is hardened a sub-sample is cored from the triaxial specimen as the resolution of the micro CT image depends on sample size.

In the research considered in the lecture the sub-sample diameters ranged from 3mm to 10mm in diameter. These small sample sizes gave a small voxel (3D pixel) size relative to sand particle sizes, thus enabling the soil particle morphology to be accurately described and particle contacts to be identified. Figure 2(b) illustrates a typical output from a micro CT of a sand sample (Fonseca et al, 2012).

DEM analysis of the state parameter

The response of soil during shearing depends on the sample density (loose sands contract, dense sands dilate) and on the stress level (increasing the stress level will suppress dilatancy).

Figure 3(a,b,c): Comparison of DEM simulations

Figure 3(a,b,c): Comparison of DEM simulations

In a Géotechnique paper, Been and Jefferies (1985) proposed that these two effects are captured by considering the difference between the current void ratio of a sand and the void ratio on the critical state line at the same mean effective stress; they termed the difference in these void ratios the state parameter (ψ).

Been and Jefferies (1985) and Jefferies and Been (2006) presented correlations between various aspects of soil behaviour and the state parameter based on empirical data. These empirical data were developed in laboratory triaxial tests. In situ soil stress states are often not triaxial and it is very difficult to establish a critical state line in the laboratory under general stress conditions.

A virtual laboratory can be created using DEM and true-triaxial tests can be simulated to large strain levels to attain a critical state line. Huang et al (2014) documented a series of true-triaxial test simulations on a model sand sample comprising spherical particles with a grading similar to that of Toyoura sand.

As illustrated in figure 3(a), the DEM simulations qualitatively captured the relationship between the peak angle of shearing resistance (ϕ’p) and the initial state parameter (ψ0); the difference in the ϕ’p values can be attributed to the spherical shape of the DEM particles. Interestingly, the DEM simulations quantitatively captured the relationship between both the difference between ϕ’p and the critical state angle of shearing resistance (ϕ’cv) and ψ0, see figure 3(b), and the relationship between the dilatancy at peak (Dp) and ψ0, see figure 3(c).

These results are important as many constitutive models for sand incorporate the state parameter concept. These data show that DEM simulations of soil behaviour under complex stress states that cannot easily be attained in the laboratory can be used to supplement triaxial tests data to develop constitutive models. Huang et al (2014) successfully calibrated the Norsand constitutive model against their DEM dataset.

Particle-scale analysis of a locked sand

Cresswell and Powrie (2004) documented a series of triaxial tests on samples of Reigate sand. Reigate sand is a locked sand. Analysis of geological thin sections revealed that the contacts between the sand grains include the concavo-convex contacts and long straight contact indicative of this type of sand.

Cresswell and Powrie carried out triaxial tests and compared the response of intact (undisturbed) samples with the response of reconstituted samples with similar void ratios. They found that for these similar void ratios, samples tested at the same confining pressure gave differing response characteristics with the intact samples mobilising higher peak stress ratios; they were also significantly stiffer and more dilatant. This difference in behaviour falls outside of conventional understanding of soil response and violates the observations of Jefferies and Been discussed above.

Fonseca et al (2013) used micro CT to provide quantitative explanation of the differences in response. They obtained sand samples from the same location as that used by Cresswell and Powrie. They followed the approach of Cresswell and Powrie and succeeded in carefully carving out undisturbed samples and using the trimmings to create reconstituted samples with the same void ratio.

Data from triaxial tests on these samples was in agreement with the earlier observations. A number of triaxial tests were repeated at the same confining pressure, but stopped at different strain levels, when they were impregnated with an epoxy resin to hold the fabric in place. Four stages were considered: (1) prior to isotropic compression, (2) at the onset of dilation, (3) at the start of a visible shear band, and (4) approaching a critical state. Once the samples had hardened, sub cores were extracted for micro CT scanning.

Figure 4(a,b,c): Particle-scale analysis

Figure 4(a,b,c): Particle-scale analysis

The micro CT data provided a rational explanation for the differences in the responses of the intact and reconstituted materials. The reconstituted material was found to have a slightly different particle size distribution to the intact material; this was likely caused by particles fragmenting along pre-existing cracks in the sand grains, see figure 4(a). More significant quantitative differences were observed when the fabrics of the soils were compared. The coordination number – number of contacts per particle – was lower in the reconstituted material than in the intact material.

Figure 4(b) gives the cumulative distributions of coordination number at each of the stages described above for both sample types. As shearing progressed there was a reduction in the coordination number for both materials, with the intact sand’s coordination number approaching the reconstituted sand’s coordination number at large strains. The contact index (proportion of the particles’ surfaces that is occupied by contacts) was a better discriminator of the differences in fabric between the intact and reconstituted materials.

Figure 4(c) gives the cumulative distributions of contact index for both materials at each of the four stages considered and contact index of the intact material was much higher than the reconstituted contact index and the contact index reduced during shearing. This study illustrates how micro CT can unravel perplexing soil behaviour.

Internal instability of dam filters

As outlined by Skempton and Brogan (1994), seepage-induced internal instability occurs in cohesionless gap graded materials with fines contents less than 25%. When the finer material is eroded from the void space that exists between the coarser particles that form the primary, stress transmitting fabric of the material. This phenomenon is known to present a hazard to dams and flood embankments.

Figure 5: Variation in stress-reduction factor with fines content

Figure 5: Variation in stress-reduction factor with fines content

Skempton and Brogan showed that this can occur at hydraulic gradients that are significantly lower than the critical hydraulic gradient for base shear (one) which is typically used in engineering design. They hypothesised that the finer particles in this type of material carry only a small proportion of the applied stress and so the effective stress in these finer particles becomes zero at hydraulic gradients smaller than one. Skempton and Brogan proposed that in materials with fines contents exceeding 35% the finer fraction separates the coarse particles and is not susceptible to preferential erosion.

Shire et al (2014) revisited Skempton and Brogan’s ideas using DEM. Samples with different fines contents and densities were considered. Each sample was taken to an isotropic stress state of 50kPa. The contact force data available from the DEM analyses were used to determine the stress in the finer grains. Figure 5 is a plot of the proportion of the applied effective stress carried by the finer material (a) against fines content. This figure largely confirms Skempton and Brogan’s hypotheses.

Figure 6: Variation in stress in finer fraction with density

Figure 6: Variation in stress in finer fraction with density

At fines contents less than 24%, the finer material carries only a small proportion of the applied stress; at fines contents exceeding 35% the finer material does not carry a particularly low stress. At fines contents between 24% and 35%, the stress transmitted by the finer material is density-dependant. Figure 6 compares the particle stresses in equivalent loose and dense samples with 30% fines. For the loose sample, almost all of the fine particles are shaded in blue and have a values of zero indicating they are not stressed, while many of the fine particles are shaded in red in the dense sample indicating they are fully stressed, transmitting a stress equal to the applied confining pressure. This study demonstrated the ability of DEM to confirm particle scale hypotheses that have an impact on design guidelines.


Engineers have long understood that the complexity of soil behaviour can be at least partially attributed to the inherent particulate nature of the material. Hitherto engineers have proposed hypothetical particle scale mechanisms to explain the overall material response. The advent of DEM and micro CT enables a more scientific assessment the particle-scale mechanics.


This article has presented work of four former PhD students: Liang Cui (now based at the University of Surrey), Xin Huang (now based at Tongji Univeristy, China), Joana Fonseca (now based at City University in London), and Tom Shire (now based at Imperial College).

Support was provided by other colleagues and former colleagues at Imperial College, particularly Kevin Hanley. Funding for the research was provided by the Irish Research Council for Science Engineering and Technology and the Engineering and Physical Sciences Research Council.


Been, K and Jefferies, MG (1985). A state parameter for sands. Géotechnique 35, 99–112

Bym, T, Marketos, G, Burland, JB and O’Sullivan, C (2013). Use of a two-dimensional discrete element model to gain insight into tunnelling-induced deformations, Géotechnique 63, 791-795

Cresswell, A and Powrie, W (2004). Triaxial tests on an unbonded locked sand. Géotechnique 54, 107-115.

Cui, L and O’Sullivan, C (2006). Exploring the macro- and micro-scale response characteristics of an idealised granular material in the direct shear apparatus, Géotechnique 56, 455-468

Cundall, PA and Strack, ODL (1979). A discrete numerical model for granular assemblies, Géotechnique. 29, 47–65

Fonseca J, O’Sullivan C, Coop MR, and Lee PD (2012). Non-invasive characterization of particle morphology of natural sands, Soils and Foundations, 52, 712-722

Fonseca, J, O’Sullivan, C, Coop, MR, and Lee, PD (2013). Quantifying the Evolution of Soil Fabric during Shearing using Scalar Parameters, Géotechnique 63, 818-829

Fonseca, J, Sim, WW, Shire, T, and O’Sullivan, C (2014). Micro-structural analysis of sands with varying degrees of internal stability, Geotechnique 64, 405-411

Huang, X, O’Sullivan, C, Hanley, KJ, and Kwok, CY (2014). DEM Analysis of the State Parameter, Géotechnique 64, 954-965

Jefferies, MG and Been, K (2006). Soil liquefaction: a critical state approach. Taylor and Francis.

Plimpton, S. (1995). Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117, 1–19.

Have your say

You must sign in to make a comment

Please remember that the submission of any material is governed by our Terms and Conditions and by submitting material you confirm your agreement to these Terms and Conditions. Links may be included in your comments but HTML is not permitted.