Force chain and contact cycle evolution in a dense granular material under shallow penetration

The mechanical response of a dense granular material submitted to indentation by a rigid flat punch is examined. The resultant deformation is viewed as a process of self-organisation. Four aspects of the mechanical response (i.e. indentation resistance, failure, Reynolds’ dilatancy, the undeforming “dead zone”) are explored with respect to the linear and cyclic structural building blocks of granular media self-organisation: force chains and contact network cycles. Formation and breaking of 3-cycle contacts preferentially occur around and close to the punch uncovering a “dilation zone”. This zone encapsulates (i) most of the indentation resistance and is populated by force chains consisting of six or more particles, (ii) all buckling force chains, and (iii) a central, near-triangular, undeforming cluster of grains beneath the punch face. Force chain buckling is confined to the zone’s outer regions, beneath the corners and to the sides of the punch where surface material heave forms. Grain rearrangements here involve the creation of 6-, 7-, and 8-cycles – in contrast with Reynolds’ postulated cubic packing rearrangements (i.e. 3-cycles opening up to form 4-cycles). In between these intensely dilatant regions lies a compacted triangular grain cluster which moves in near-rigid body with the punch when jammed, but this dead zone unjams and deforms in the failure regimes when adjacent force chains buckle. The long force chains preferentially percolate from the punch face, through the dead zone, fanning downwards and outwards into the material.


Introduction
The contact between a granular material and a solid body lies at the core of essentially all human activity on the Earth's surface.Structures, vehicles and machines operate on geological materials which are granular in form, e.g.soil, rock, ice rubble, snow (Paavilainen and Tuhkuri, 2013;Estep and Dufek, 2013;Murthy et al., 2012;Calvetti et al., 2004).Laboratory tests which mimic, and are aimed at understanding the fundamentals of, deformation and dewatering processes in geological formations involve mechanical configurations where the geomaterial is loaded in fully confined or partially confined solid boundaries (e.g.Carson and Berglund, 1986).In most of these aforementioned cases, the structure or mechanical system not only imposes forces on the surface of, but also penetrate into, the granular material.A prevailing message in these and other studies of geological hazards (e.g.rock debris De Blasio and Crosta, 2014 and dense volcanic flows Estep and Dufek, 2013, ice rubbling Paavilainen and Tuhkuri, 2013 etc.), all aided by DEM simulations, is that force chains, regardless of breakage, provide the primary mechanism for load transfer, and their failure by buckling governs failure and the attendant energy release at higher spatial scales (De Blasio and Crosta, 2014).
Granular-solid interaction systems, under shallow and deep penetration conditions, are especially widespread in certain sectors of industry, including land contaminant remediation, mining, land transportation, infrastructure, agriculture and urban development.The particular type of granularsolid interaction system of interest here is the classical contact mechanics problem of indentation (shallow penetration) of a material by a rigid punch (Johnson, 1985;Tordesillas and Shi, 1999): see, for example, the experiments and physical models of foundations and embankments by Biarez (1962) (Fig. 1).This problem is of core importance not only in the geophysical and geotechnical operations mentioned above, but also in the manufacturing and processing of powders and grains in the ceramics, chemical, pharmaceutical and agricultural industries (e.g.Sanchez-Arevalo et al., 2013;Kadiri and Michrafy, 2013).
Granular-solid interaction systems give rise to a plethora of complex and emergent phenomena (e.g.Duran, 1999;Sandnes et al., 2011).Among the most studied is the formation of self-organised load-bearing particle groups called "force chains" which vary in size from a few to many grains that can percolate through the whole sample (e.g.Geng et al., 2003;Oda and Kazama, 1998;Oda et al., 2004;Tordesillas, 2007).From the standpoint of material characterisation and modelling, these emergent structures pose great scientific challenges as their analytical formulation embodies strong nonlinearities and multiscale dynamics in space and time.The root cause of this complexity are the vast number of internal degrees of freedom and the many-body interactions that govern their evolution (Tordesillas and Muthuswamy, 2009).In this context, robust model predictions of granularsolid interaction hinge crucially on a rigorous characterisation of these interactions, both between constituent grains in a single force chain and between force chains in the bulk material.In recent earthquake centrifuge tests, which investigated the response of earth structures during earthquakes, these interactions manifested in fabric evolution -with associated force chain evolution implicated as the key mechanism behind the different liquefaction resistances observed for the various sand samples (Yu et al., 2013).Yu and co-workers highlight the need to take into account the micromechanics of force chains in future research into liquefaction and soilstructure interaction under seismic loading conditions.Similarly, experimental studies of microstructural fabric evolution during dewatering in saturated granular materials under compression suggest that shear failure zones "become both loci of water escape and zones of differential movement" (Carson and Berglund, 1986).Recent complex network analysis of interstitial void evolution in a dry granular material corroborates this observation, viz.that the formation of shear zones, governed by the dynamics of force chains, places grains in a configuration more prone to the efficient transport of material through its pore space, as evidenced by the location of percolating shortest path routes (Walker et al., 2013).The complementary but inter-related nature between the solid-grain matrix and interstitial pore space has also been recognised in studies of vibrated powders (Barker and Mehta, 1992).
In this study, we aim to unveil new insights into the internal (mesoscopic) processes responsible for the mechanical response of a dense granular material to a symmetric indentation by a rigid flat punch.Our strategy is to treat the resultant granular deformation as a process of self-organisation.To achieve this we primarily represent the granular system as a complex network.We have already established the structural building blocks of granular media self-organisation (Tordesillas et al., 2010).The structural building blocks comprise linear and cyclic mesoscopic structures, i.e. the major load-bearing force chains and their surrounding contact cycles, respectively.We also uncovered cooperative behaviour between these emergent structures (Tordesillas et al., 2010(Tordesillas et al., , 2011)).The nature of this cooperative interplay is reminiscent of rules employed in the design of trusses and columns in architectural high-rise constructions (e.g.communication towers, power towers, offshore platforms, bridges, cranes, etc.).In essence, the material self-organises to form "truss-lacedcolumns" with high compressive load-bearing capacity (Ali and Moon, 2007;Keii et al., 2010).Triangular trusses, i.e. the 3-cycles, brace the force chain columns both to increase rigidity where needed (i.e.prop-up column misalignments) and resist relative grain rotations, the precursory mechanism for force chain failure by buckling.The combination of the linear building blocks of force chains and the cyclic building blocks of contact n-cycles suggests parallels with the concept of linear and complex bridge, or arch, structures found in granular materials (Mehta et al., 2004;Carlevaro and Pugnaloni, 2012;Cao et al., 2013).To date, only qualitative connections have been suggested, with much work required to realise a precise characterisation of the relationship between bridges and the force chain and n-cycle building blocks we examine in detail here.
Broadly, the question we seek to answer in this study is: How does the load-indentation curve -i.e. the principal measurement in tests of the macroscopic response of materials under this type of boundary condition -relate to the evolution of structural building blocks inside the material?Specifically: what aspects of this evolution govern the following key aspects of mechanical response: (i) penetration resistance, (ii) dilatancy, (iii) loss of load-carrying capacity, and (iv) the development of a "dead zone" of undeforming material?Indeed, the development of dilatancy in this system is synonymous with the circumstance surrounding Reynolds' discovery of this phenomenon while walking along the beach: the wet sand, indented by his foot, dilated causing the adjacent surface grains to "whiten" as their wet gloss receded.Reynolds (1885) attributed this illusion of the sand drying up around his foot to specific grain rearrangements that created large open pores, into which the surface water drained under gravity.We aim to quantify these dilatant grain rearrangements with respect to the evolution of structural building blocks.
The rest of this paper is arranged as follows.A brief summary of the salient aspects of the simulation is given in Sect.2, followed by a delineation of our strategy for material characterisation using complex networks in Sect.3. We present our results in Sect.4, and then conclude with a summary of key findings and open questions for future research in Sect. 5.

Discrete element simulations
We confine this discussion to a summary of the salient features of the simulations.Complete details of the penetration model, constructed using the discrete element method (DEM), are provided in Muthuswamy and Tordesillas (2006).This DEM model has been employed to examine the constitutive response of granular assemblies subject to a range of monotonic and cyclic compression and penetration tests (e.g.Muthuswamy and Tordesillas, 2006;Tordesillas, 2007;Tordesillas and Muthuswamy, 2009).Moreover, this model has been subjected to a comprehensive comparative analysis of multiple aspects of granular deformation and failure (e.g.dynamics of force chain evolution, particle kinematics and contact evolution) -against various experimental tests on assemblies of photoelastic disks and on sand (e.g.Tordesillas et al., 2013Tordesillas et al., , 2011)).
Resistance to relative motion at the particle-particle and particle-wall (punch and confining box) contacts is governed by combinations of a linear spring, a dashpot and a friction slider.The model is designed to mimic the response of assemblies of noncircular particles.It incorporates a moment transfer to account for rolling resistance, in accordance with Oda and Iwashita (2000).This modification to the classical DEM model of Cundall and Strack (1979) which allows the relative rotations of particles at contacts to be controlled, has been found to be essential in achieving more realistic rotations and stress predictions (e.g.Gutierrez and Mohamed, 2010).
A combination of Hooke's law and Coulomb's law defines the interaction between contacting particles, namely where f n and f t are the normal and tangential components of the contact force, k n and k t are spring stiffness coefficients, b n and b t are the viscous damping coefficients, and µ is the Coulomb friction coefficient.Following the modified DEM in Oda and Iwashita (2000), the rolling resistance or contact moment, defined in an analogous fashion to Coulomb's law, is expressed as where R min denotes the smaller of the radii of the two interacting particles, k r is the spring stiffness coefficient, and µ r is the friction coefficient.The remaining quantities in the above relations are the relative normal and tangential displacements and relative rotations denoted respectively by u n , u t , and α, as well as the relative normal and tangential translational and rotational velocities: v n , v t , and α.Table 1 provides a summary of the simulation and material parameters used.The stiffness coefficients -rolling resistance and contact moment -are selected according to the assumption that subject to equilibrium conditions, contact moments due to rolling resistance are comparable to contact moments due to tangential forces (Iwashita and Oda, 1998) coefficients for normal and tangential forces are given by b n,t = 0.1 , where m min is the mass of the smaller contacting particle.Analogously, b r = 0.1R min √ m min k r .The granular bed in its undeformed configuration is shown in Fig. 2. The dimensions of the granular bed and ratio of punch width to box width (1 : 10) were chosen in order to reduce boundary effects (see Muthuswamy and Tordesillas, 2006).The punch is pushed down into the granular bed at a constant increment of vertical displacement (i.e.0.2 % of box height).A total of 64 displacement increments is applied, with the granular material in mechanical equilibrium at each stage of the quasi-static loading program.Since we are considering shallow penetration only, the punch indents to 12.4 % of the box height.At each stage, i.e. a fixed indentation depth, data acquired for each grain i in the material in equilibrium comprise centroid position, kinematics, contacts, and contact forces.The macroscopic mechanical response is shown in Fig. 3.The loading history comprises six regimes, alternating between periods of increasing resistance and failure, as evident in the rises and relatively brief falls in the total force exerted by the granular material on the punch.Region 2 is the first failure regime from the beginning of loading.The conditions here mimic those encountered in geotechnical engineering studies of incipient failure beneath a foundation.In such studies, this failure point determines the ultimate bearing capacity of the underlying soil, i.e. the soil's capacity to withstand the load of the structure.We confirmed the trends reported here are robust for other punch systems (i.e. a range of packing fractions from 0.894 to 0.899; number of particles from 11 743 to 12 364).

Material characterisation using complex networks
Our strategy for material characterisation proceeds in three phases.In Phase 1, we map the data relating to the physical or simulated sample to a set of complex networks, each representing the sample at a particular stage of the loading history.In Phase 2, we analyse the properties of these complex networks.In Phase 3, we map the results of the complex network analysis back to the domain of the sample, and extract new physical insights into material behaviour with respect to the question posed.
The question under consideration is centred around (i) penetration resistance, (ii) dilatancy, (iii) loss of loadcarrying capacity, and (iv) the development of a "dead zone" of undeforming material.Although these aspects of mechanical response are intertwined, (i) and (iii) deal directly with force transmission, while (ii) and (iv) can be quantified through grain arrangement and rearrangements.As indentation proceeds, the particles rearrange: some contacts break while new contacts form.The material's deformation can thus be entirely tracked through its evolving network of contacts.Also, forces are transmitted through these contacts.Thus the complex network to be constructed from data in Phase 1 is the contact network: here the nodes represent the grains, and the links represent the contacts between grains.
In Phase 2, we are to identify the structural building blocks of self-organisation and characterise their evolution.The cyclic structural building blocks, i.e. the n-cycles, are found from the minimal cycle basis of the contact network (Horton, 1987;Mehlhorn and Michail, 2006).The identification of the linear structural building blocks of force chains require additional information, namely, the contact forces as well as their orientation.The methods we employ for identifying force chains and n-cycles are explained in detail in Appendix A. Figure 4 shows an example of these conjoined mesoscopic structures.In this study, we also quantify the failure of structural building blocks and details of the methods used to identify them are also given in Appendix A.
In Phase 3, we tie the evolution of these structures to the four aforementioned aspects of mechanical response.Given the immediate physical significance of the complex contact network, the network patterns and structures obtained in Phase 2 lend themselves well to physical interpretation.

Results
In what follows, we first focus on the temporal evolution of structural building blocks in conjunction with the four aspects of penetration resistance, failure, Reynolds' dilatancy and the existence of a near-undeforming cluster of grains (i.e. the dead zone).Next, the spatial distributions of building blocks are explored, in particular, how the different size building blocks are organised and assembled together to form one global structure that resists the indentation.Finally, we tie these spatial trends to the earlier results on temporal evolution to understand how the material evolves its structure.What building blocks does it create more versus less under increasing indentation?How are the different building blocks adapted and reorganised to support the punch as it penetrates deeper into the bed?

Temporal evolution of structural building blocks
We present in Fig. 5 the evolution of the structural building blocks of (a) m-force chains and (b) n-cycles, classified according to the number of particles (m or n) that form the structure.We observe a general decline in the population of force chains in each class with increasing indentation.By contrast, the evolution of the minimal contact cycles differ according to size n: the smallest members (i.e.3-cycles) decrease in population while larger cycles (n ≥ 6) increase in number.This trend suggests an overall increase in the porosity of the material (i.e.Reynolds' dilatancy) as loading proceeds: the small 3-cycles cleave to form large cycles.

Indentation resistance
The evolution of the material's resistance to the punch indentation, shown in Fig. 3, comprises interlacing periods of increasing resistance and failure, as evident in the rises and relatively brief falls in the total force exerted by the material on the punch.This resistance is mostly borne by the particles in force chains, the principal load-bearing structures in a deforming granular material (Peters et al., 2005).The major (λ + ) and minor (λ − ) principal stresses acting on a given particle can be computed as described in Appendix A. The magnitude of the "particle load vector" for particle i is the magnitude of the most compressive principal stress.This particle property, averaged over different subgroups of particles in the bed, provides a measure of that subgroup's resistance to the punch indentation.On average, we observe a distinct difference in distributions of this property for short mforce chains, m < 6 and long m-force chains with m ≥ 6 (see, Fig. 6).These longer force chains provide the greatest resistance to the increasing punch load.This load allocation is far from obvious given the previously established cooperative evolution between force chains and 3-cycles (Tordesillas et al., 2010).As shown in Fig. 7, except for the first few stages of loading, the short force chains reside in those parts of the contact network with the highest average degree (i.e.supporting neighbours) and highest truss-like 3-cycles.Initially, the long force chains have higher 3-cycle supports compared to other non-force chain particles, but this trend is reversed in the latter stages of loading, counter to all past observations of mechanical response in biaxial and triaxial compression tests (Tordesillas et al., 2010(Tordesillas et al., , 2011(Tordesillas et al., , 2012)).In such tests, 3-cycle membership of force chain particles is consistently higher or at least in parity with those of non-force chains at all stages of loading, including during failure.To unravel the reasons why, we examine the spatial distributions of these structures later in the next section.

Loss of load-carrying capacity
Force chains, like architectural columns, fail by buckling.
The degradation of the lateral supporting structures of 3cycles, evident in Figs. 5 and 7, render force chains particularly prone to collapse by buckling.The results shown in Fig. 8 provide unambiguous evidence of the underpinning mechanism for bulk failure.Observe here that the sharp increase in the population of mesoscopic failure events, i.e. the failure of the structural building blocks of 3-cycles, precipitates the failure of force chains by buckling and the sudden loss of load-carrying capacity of the sample at the bulk macroscopic level, i.e. drops in P over regimes 2, 4,  and 6.Also note that in regime 1, i.e. the regime prior to incipient failure (regime 2), there are large fluctuations in the number of 3-cycles before the sharp increase at failure (Fig. 8b), while the force chain columns maintain their integrity throughout (Fig. 8a): this pattern holds true in subsequent periods of increasing resistance and failure (regimes 3-4 and 5-6).

Reynolds' dilatancy
Densely packed granular materials expand in volume when compressed under shear.The first reported study of this phenomenon was by Osborne Reynolds (Reynolds, 1885).Reynolds' illustration of the grain rearrangements that caused the dilatancy beneath, and to the sides of, his foot is reproduced in Fig. 9, along with its analog in planar deformation.Rowe (1962) also presented a similar view in relation to Coulomb's slip plane (shear band) for uniformly sized spheres.In these cases, the expansion in the interstitial void space translates to 3-cycles opening up to form 4-cycles.However, dilatancy in the present system is due to the creation of even larger cycles, particularly n = 6 (Fig. 5d), and this occurs concomitant with the progressive loss of 3-cycles both throughout the system (Fig. 5c).and around force chains in Fig. 7.The population of 4-cycles and 5-cycles remain almost constant throughout loading.
To unravel further details on the mechanisms giving rise to dilatancy, as well as penetration resistance, failure and the dead zone, we need to examine the spatial distributions of structural building blocks.

Spatial evolution of structural building blocks
In all past studies of the structural building blocks of m-force chains and 3-cycles, we have shown that these structures exhibit cooperative behaviour, irrespective of loading condition, material properties and dimensionality (Tordesillas et al., 2010(Tordesillas et al., , 2011(Tordesillas et al., , 2012)).Here we uncover new insights into the emergence and cooperative interplay among these structures -in different regions of the sample.Consider a spatial region extending radially inwards to the granular bed from the centre of the punch face to a boundary radius r.For a given r we can identify particles within this region and outside of this region.For all particles inside the region we can find the number of 3-cycles they are a constituent member.We define the average 3-cycle membership in this radial band region to be the average of the 3cycle membership counts per particle in the band.Examining Fig. 10 we see that for non-zero indentation time steps and radial bands extending to five punch widths (i.e.r = 5w), all curves have a common point of inflection.The radius (i.e.r = 2.5w) at which this inflection point occurs defines a region we hereafter refer to as the "dilation zone".We observe and will show a distinct difference in the behaviour of m-force chains and 3-cycles inside and outside of .For example, the percentages of particles with respect to the whole assembly residing inside versus outside this dilation zone remains almost invariant throughout loading, however, changes to the average 3-cycle membership per particle is essentially confined to the dilation zone (see, Fig. 11).The width of this dilation zone is comparable to the lateral extent of the deformed region identified in experimental studies (Murthy et al., 2012).The distance at which this inflection point occurs provides an indication of the diminishing influence of the punch.

Indentation resistance
Typical spatial distributions of m-force chains are given in Fig. 12.The long m-force chains, m ≥ 6, increasingly concentrate directly below the punch, fanning downward and outward from the punch face, as indentation proceeds.The spatial distribution of these chains is similar in shape to that observed in the experiments of Geng et al. (2003), in particular, the mean field response (i.e.averaged over many realisations) to a concentrated force at the boundary of an assembly of bidisperse photoelastic disks.In contrast, the short m-force chains, m < 6, are more evenly spread throughout  the material with less starkly apparent preferential areas of concentration at any stage of the loading history.Inside the dilation zone, the percentage of particles in force chains start from 34 % and progressively rises to 60 % at the last stages of loading (data not shown).Most of these chains are long force chains: short force chains in the dilation zone remain almost constant at around 15 %, while long force chains reach as high as 45 % in the last stages of loading.On average, the long force chains in the dilation zone provide the greatest resistance to punch indentation prior to incipient failure (regime 2); however, they bear almost the same, albeit marginally higher, loads compared to the short force chains for the rest of the loading history (Fig. 13).

Failure and dilatancy
All failure of force chains by buckling occurs inside the dilation zone .More precisely, these failure events are concentrated near the corners of the punch, downward and outward from the punch corner and sides, as shown in the accumulated spatial distributions of buckling force chains for all regimes of loading (Fig. 14a).Unlike force chains, which exclusively fail by buckling inside , there are a small number of 3-cycle deaths occurring outside .However, as can be seen in Fig. 14b, the great majority of deaths in 3-cycles throughout the entire loading process occur inside the dilatation zone.
It is interesting to note that the locations of force chains failing by buckling coincides with the dilation regions identified by Reynolds, i.e. adjacent to the indenting body.Recall from Figs. 11 and 14b that the degradation of 3-cycles is essentially confined to the dilation zone.Even more revealing of the material evolution with indentation is the spatial distribution of a special subset of the 3-cycle population, i.e. the most stable, so-called "persistent 3-cycles" (Fig. 15).Persistent 3-cycles are those which have lived from the beginning of loading to the stage in question.In previous studies of biaxial and triaxial compression tests, the destruction of the persistent 3-cycles in the stages leading up to the onset of shear banding was confined to that region where the shear band ultimately develops (Tordesillas et al., 2014).Furthermore, these persistent structures become essentially depleted in the shear band once the band is fully developed.Here the progressive degeneration of persistent 3-cycles is confined to the dilation zone, with the population reaching negligible values by the final failure regime of the loading history, i.e. regime 6 (Fig. 15a).The loss of these stable, persistent 3-cycle structures has implications for particle rotations and, in turn, force chain buckling.As previously mentioned, relative rotations are precursors to buckling and truss-like 3-cycles frustrate these at contacts.Accordingly, we expect the dilation zone to be the region where particle rotations concentrate.As shown in Fig. 15b, the development of rotations is similarly progressive in, as well as being mainly confined to, the dilation zone.Moreover, the instantaneous rotations (i.e.rotations developed over a single indentation increment) are temporally correlated with the failure of 3-cycles and force chains by buckling.Note the pronounced peaks in instantaneous rotations during the three failure regimes of 2, 4, and 6.This corroborates the findings in Fig. 8: the cleavage of 3-cycles result in an increase in the intensity of rotations and, in turn, force chain buckling.Based on the above findings, we now offer a new perspective on Reynolds' dilatancy in surface indentation (shallow penetration) systems.The phenomenon of dilatancy arises from the combined failure of 3-cycles and force chain buckling.The rearrangements associated with buckling are, however, distinct from the symmetric cubical topology proposed by Reynolds' in his original manuscript (recall Fig. 9).Force chain buckling is a highly asymmetric rearrangement event which typically creates larger n cycles.

Dead zone
Recall from Fig. 3 that the mechanical response comprises a series of rises in resistance, followed by relatively brief stages of failure in which a drop in load-carrying capacity can be observed.During periods of increasing resistance to indentation, grain motions relative to the indenting punch remain small: this is evident in the low values of the magnitude of relative displacement vector of grain i with respect to the punch, |v i − v p |.For particles with 3-cycle membership values ranging from zero to six, the width of these distributions during the stable regimes 1, 3, and 5 is almost constant (see, for example, regime 5 in Fig. 16).By contrast, during the failure regimes, the distribution of |v i − v p | widens, suggesting that significant deviation in grain motion relative to the punch develops for those grains with low 3-cycle membership (see regime 6 in Fig. 16).Those with high 3-cycle membership, i.e. at least four, mostly remain in sync with the punch motion; |v i −v p | remains low for most of these grains.Zooming in to the region beneath the punch, we find a cluster of these grains: see [inset, left] of Fig. 16.Note that there are other grains which exhibit very low |v i −v p | and high 3-cycle membership, but these predominantly lie outside the dilation zone.During regimes 1, 3, and 5, this cluster of grains manifest the near-triangular shaped region of undeforming material, i.e. the dead zone (recall Fig. 1b).However, during failure regimes 2, 4, and 6, connectivity among these grains degrades (Fig. 16,inset,right).This structure, however, is neither persistent nor possesses a well-defined shape.
In summary, the overall mechanical response to the punch penetration is captured in Fig. 17 Close-up of spatial distribution of grains having at least four 3cycles near the punch: regime 5 (inset, left); regime 6 (inset, right).The 3-cycle membership is the number of 3-cycles that grain i is part of; |v i − v p | is the magnitude of v i − v p , where v i is the displacement of grain i and v p is the displacement of the punch.
punch, transmitting relatively low contact forces.Between these intensely dilatant regions, directly below the punch face, is a region that compacts and embodies the dead zone.
Rearrangements here mostly involve the reduction of interstitial voids through an increase in the number of 3-cycles.
The long force chains percolate from the punch face, through the dead zone that is fortified with truss-like 3-cycles, fanning downwards and outwards into the material: see Figs. 12 and 17.

Conclusions
The mechanical response of a dense granular material, subject to indentation by a rigid flat punch, was examined as a process of self-organisation.In this framework, the deformation of the material evolved through the continual formation and collapse of linear and cyclic structural building blocks of self-organisation: force chains (determined by structural mechanics) and contact cycles (delivered by complex network analysis).In essence, the material forms trusslaced columns that are robust to axial compression and shear.Four aspects of mechanical response were explored with respect to the columnar force chains and truss-like contact cycles, i.e. indentation resistance, failure, Reynolds' dilatancy, and the undeforming "dead zone".With increasing indentation, we observed a degeneration of particles in 3cycles concomitant with a decrease in the number of particles participating in the force chain network.Despite the overall decline in 3-cycle population, there is considerable birth and death of 3-cycles occurring preferentially around and close to the punch.We uncovered a "dilation zone", where the most stable, persistent 3-cycles exclusively and progressively degrade, until essentially depleted.Recent studies of the evolution of localised failure in biaxial and triaxial compression tests showed that the shear band ultimately develops in this zone.Here, this zone is centred at, and extends radially into the material to ≈ 2.5 times the punch width from, the initial punch centreline.This zone encapsulates (i) most of the indentation resistance and is populated by long force chains, (ii) all buckling force chains, and (iii) a near-triangular cluster of grains, directly beneath the punch.Buckling is confined to the outer regions, beneath the corners and to the sides of the punch where surface material heave forms.Reminiscent of Reynolds' discovery of dilatancy from his foot's indentation of wet sand, the topological changes that give way to this surface heave are characterised by the creation of large open pores.However, our findings show the creation of not only larger but also highly asymmetric cycles than imagined by Reynolds.Grain rearrangements to the side of the punch involve the cleavage of 3-and 4-cycles to as large as 6-, 7-, and 8-cycles -in contrast with Reynolds' postulated cubic packing rearrangements (i.e.3-cycles opening up to form 4-cycles).In between these intensely dilatant regions of high porosity, lies a compacted triangular grain cluster which moves in near-rigid body with the punch when jammed (i.e.high levels of 3-cycles).The long force chains percolate from the punch face, through the dead zone, fanning downwards and outwards into the material.During stages of failure, when force chains buckle, the energy released unjams and www.nonlin-processes-geophys.net/21/505/2014/ Nonlin.Processes Geophys., 21, 505-519, 2014 deforms the dead zone.There are other structures that emerge in penetration systems, besides those studied here.For example, vortices were recently observed in Murthy et al. (2012); Hossain and Fourie (2013) in sand experiments.These merit future study, given the recent findings on the role that they play in the initiation of localised failure Tordesillas et al. (2014).
Complex network techniques have shown great promise in unravelling multiscale and nonlinear problems in geophysics.Recent complex network analysis of the interstitial void space (Walker et al., 2013) and the solid grain phase (e.g.Smart and Ottino, 2008;Tordesillas et al., 2010Tordesillas et al., , 2012Tordesillas et al., , 2013) ) provide compelling impetus for its future application in solving fundamental problems germane to a broad range of geophysical phenomena, including comminution, dewatering, hydraulic fracturing, fault gouge dynamics, and localised and diffuse mode failure in geomaterials.

Fig. 1 .
Fig. 1.Experiments of Biarez (1962): development of deformation field beneath an indenting flat punch relative to a fixed reference frame (a); and close-up of undeforming so-called "dead zone" with image taken relative to the moving punch (b).

Fig. 2 .Fig. 3 .
Fig. 2. The punch system at the initial state.Particles are represented by a point at their centre.

Fig. 4 .
Fig. 4. (Colour online)A close-up of particles in a 5-force chain (gray filled circles) beneath the punch and its conjoined n-cycle membership (filled polygons), taken from one of the stages within regime 2 of loading.Counterclockwise from the top left side of the force chain are three 3-cycles (red), a 5-cycle (green), a 3-cycle (red), a 4-cycle (blue), a 6-cycle (yellow), and an 8-cycle (cyan).Confining neighbour particles of the force chain are drawn as open circles overlaid with the contact network.

Fig. 8 .
Fig. 8. (Colour online) Punch load P and failure of structural building blocks.Population of force chains that have buckled by at least an angle θ * where θ * = 15 • , 10 • , 15 • (a), and 3-cycles that cleaved (b).Population of failure events are computed over a uniform sequence of indentation increments of size 1 and plotted at the end of the indentation increment.

Fig. 9 .
Fig. 9. (Colour online) Classical view of grain rearrangements associated with dilatancy.Reynolds (1885) illustration of dilatancy for monodisperse rigid spheres: from a face-centred cubic packing on the left to a simple cubic packing on the right (cycles are added) (top left).Corresponding grain rearrangement in planar deformation with two 3-cycles opening up to a 4-cycle (bottom left).Rowe (1962) illustration of grain rearrangement for shear band formation in a regular packing (right).

Fig. 10 .
Fig. 10.(Colour online) Average 3-cycle membership of particles versus strip outer radius for different indentations.Inset: the dilation zone (w is the width of the punch).

Fig. 11 .
Fig. 11.(Colour online) Average 3-cycle membership of particles versus indentation: inside and outside of the dilation zone.Percentage of particles relative to total number of particles in system versus indentation: inside and outside of the dilation zone.

Fig. 16 .
Fig.16.Summary statistics, mean (solid) and standard deviation (dotted), of representative distributions of |v i − v p | versus 3-cycle membership of grain i: at end of regime 5 and at end of regime 6. Close-up of spatial distribution of grains having at least four 3cycles near the punch: regime 5 (inset, left); regime 6 (inset, right).The 3-cycle membership is the number of 3-cycles that grain i is part of; |v i − v p | is the magnitude of v i − v p , where v i is the displacement of grain i and v p is the displacement of the punch.

Fig. 17 .
Fig. 17. (Colour online)Close-up of dilatant regions to the sides of the punch at the end of regime 5, mediated by a central region of densely connected, dead zone that is rich in truss-like 3-cycles.Links in the contact network are coloured red (black) if part of 3cycles (not part of 3-cycles).Line thickness reflects the magnitude of the (depth-normalised) contact force; thicker lines are contacts transmitting higher forces.