Evolving pore orientation, shape and size in sheared granular assemblies

This paper presents new insights into the deformation response of sheared granular assemblies by characterising pore space properties from discrete element simulations of monodisperse particle assemblies in two-way cyclic shearing. Individual pores are characterized by a modified Delaunay tessellation, where tetrahedral Delaunay cells can be merged to form polyhedral cells. This leads to a natural partition of the pore space between individual pores with tetrahedral and polyhedral geometry. These are representative of small compact pores and larger well-connected pores, respectively. A scalar measure of pore orientation anisotropy during shearing is introduced. For triaxial shearing, larger pores align in the loading direction, while small pores are aligned perpendicular to the larger pores. Pore anisotropy mobilises at a slower rate than contact anisotropy or macroscopic stress state, and hence, is an important element to characterise in granular assemblies. Further, the distribution of pore volume remains isotropic. Pore shape was found to be a good micro-scale indicator of macroscopic density, with a strong relationship between averaged shape factor and macroscopic void ratio. Combining results for pore shape and orientation reveals an interesting interplay, where large elongated pores were aligned with the loading direction. These results highlight the importance of considering pore space characteristics in understanding the behaviour of granular materials.

tact networks have been extensively studied and this is aided by the fact that both can be regarded as discrete objects.
In contrast, the effect of pores (or more generally, pore space characteristics) on mechanical response of the medium has received relatively limited attention. This may be attributed to the complexity associated with the description of continuous, entirely interconnected pore space in granular materials. Prior studies have focused on how pore space characteristics affect fluid flow and transport properties, such as conductivity [3,33], water retention [4,8] and filtration [24,30].
However, features of the pore space can also provide insights into deformation characteristics in sheared granular assemblies. Oda et al. [23] observed the formation of elongated pores in bi-axial compression tests on photoelastic disks and quantified this pore anisotropy using a scan-line technique. Ghedia and O'Sullivan [6] extended this method to analyse digital images, while Muhunthan et al. [21] and Muhunthan and Chameau [20] employed a similar technique in exploring the concept of yielding and an ultimate state boundary surface.
An alternate approach is to tessellate the pore space, or identify individual pores forming the granular assembly. Common approaches include the classical Delaunay and Voronoi tessellations, which typically rely on the centroidal coordinates of particles. Other tessellation methods such as the approach proposed in Li and Li [16] explicitly consider the interparticle contact points. Space tessellation techniques have also been considered in the formulation of a micromechanical strain tensor for a granular assembly [2,14,15,27]. Li and Dafalias [18] explicitly considered a pore fabric measure in their model of anisotropic critical state based on the length of void vectors obtained from a slight modification of the space tessellation proposed by Li and Li [16]. They proposed that limiting states of pore space fabric would be attained at critical state conditions (in addition to the usual constraints of critical state soil mechanics). These examples suggests that the use of space tessellation is a useful means to gain fundamental understanding of deformation characteristics.
This study investigates how pore space characteristics affect the deformation response of sheared granular assemblies by employing the space tessellation method. Porescale data for sheared granular assemblies are extracted from numerical simulations employing the discrete element method, and individual pores are identified using the modified Delaunay tessellation [1]. While it is acknowledged that there is no unique way to tessellate the pore space, several studies [10,25,31] have suggested that the modified Delaunay tessellation provides a physically representative delineation of the pore space. Two particular properties of the pore space are investigated in detail: (i) pore orientation and the associated anisotropy of pores, and (ii) pore shape as an indicator for macroscopic density. The influence of pore size is also considered indirectly, while some elements of the contact networks are briefly discussed. Note that loose constant volume (CV) and simple shear (SS) samples were metastable and exhibited collapse mechanisms, where the samples were unable to sustain any shear stress. These metastable samples are not considered further in this study. Samples are generated by randomly inserting particles within a cuboidal domain with no interparticle contacts and subsequently isotropically compressed to 200kPa.
In the CP-L and CP-D simulations, the particles were bound by rigid frictionless walls, while periodic boundary conditions were imposed for the CV-D and SS-D simulations. The effects of gravity were not considered. All samples comprised monodisperse assemblies of particles with diameter d s = 1 mm with contacts characterised by the Hertz-history model (available in LIGGGHTS), with particle elastic modulus and Poisson ratio of E = 70 GPa and ν = 0.25. The local packing density in subsets of the particle assemblies were monitored to check if the packing density exceeded the random close packing limit. These checks found that there was no significant crystallisation in the simulations. Samples of varying density were generated by adjusting the interparticle friction angle during the sample generation process. The interparticle friction angle was set at μ = 0.25 for all samples during shearing. Each case simulated two-way cyclic shear test comprising initial loading, unloading and reloading phases in order to explore the influence of shearing direction (i.e. compression or extension), sudden load reversals and history effects on the evolving pore space features. All samples were sheared in a quasi-static manner. Further details on these numerical simulations are available in Sufian et al. [32]. The macroscopic response for all samples are reproduced in Figs. 1 and 2, which shows that the results agree well with typical, measured behaviour of granular materials. These simulations encompass a broad range of stress and volumetric conditions. The axisymmetric simulations include two different packing densities (i.e. loose and dense) and two different loading paths (i.e. drained and undrained). The direction of principal stress axes are fixed in the axisymmetric simulations, while rotation of principal axes is explored in the simple shear simulation. Further, the CP-L and CP-D simulations experience macroscopic volume change, while CV-D and SS-D maintain a constant volume constraint. Very few studies explore pore space characteristics in sheared granular assemblies for the broad range of stress and volumetric conditions considered here.

Pore space tessellation
The pore space is tessellated using the modified Delaunay tessellation proposed by Al-Raoush et al. [1]. This method initially considers the classical Delaunay tessellation, which produces a space-filling set of tetrahedral cells, where the vertices of the tetrahedron are the centres of the four particles forming the tetrahedron. The classical Delaunay tessellation over-segments the pore space, and thereby does not capture the wide variation in pore sizes (that is, the existence of small and large pores) and the variable connectivity of the pores (as tetrahedrons have exactly four four faces). Merging of these tetrahedral Delaunay cells using the modified Delaunay tessellation addresses both these concerns to produce a more physically representative tessellation of the pore space.
Each Delaunay cell has an inscribed sphere contained entirely within the pore space. Part of the inscribed sphere may be located outside the Delaunay tetrahedron and intersect neighbouring inscribed spheres. Individual Delaunay cells are then merged with adjacent cells if their respective inscribed spheres overlap. This results in a space-filling set of polyhedral pore units which are representative of individual pores. Figure 3 provides several examples of pores identified using the modified Delaunay tessellation, and further details can be found in Sufian et al. [31].
This modified Delaunay tessellation is used to track pore units for each of the granular assemblies through the two-way cycle of shearing described above. Pore unit properties are obtained at specified strain increments, with δε a = 0.083% for CP-L and CP-D, δε a = 0.016% for CV-D and δγ = 0.042% for SS-D, allowing for a detailed investigation of pore space characteristics.

Properties of individual pores
For each individual pore, there are several geometric properties of interest in this study: (i) pore volume, (ii) pore orientation, and (iii) pore shape.
Pore volume, V p , is calculated by subtracting the solid volumes from the total volume of the unit cells, and can be determined analytically with the assumption of nonoverlapping particles (which is reasonable for the numerical simulations considered in this study).
The orientation of individual pores can be determined using the pore orientation tensor introduced by Sufian et al. [31]: where, n α and a α are the unit normal vector and corresponding area of the face of the unit cell, and the summation is over all N f faces forming the pore unit. The major principal orientation of the pore, p, is determined from an eigendecomposition, where p is the eigenvector corresponding to the minimum eigenvalue of P i j . The eigen-decomposition of P i j has a physical basis, where the eigenvalues are the projected areas in the principal directions specified by the eigenvectors, as detailed in Sufian et al. [31]. The pore orientation tensor in Eq. 1 is of a similar form to the Minkowski interface tensor W 0,2 1 = S n ⊗ ndS [29] and has been considered by Schröder-Turk et al. [28] and Kraynik et al. [13] to explore cell shape anisotropy. In particular, Schröder-Turk et al. [28] considered the ratio ξ min ξ max , where ξ min and ξ max are the minimum and maximum eigenvalues of the interface tensor, to quantify shape anisotropy of Voronoi cells in disordered sphere packing. In a similar manner, this study employs the shape factor: Numerical stress-strain and stress path response for the CV-D and SS-D simulations [32] (a) (b) (c) (d) Fig. 3 a Tetrahedral Delaunay cell, or pore unit with no merging (termed unmerged in this paper); b pore unit formed by merging two Delaunay cells; c pore unit formed by merging four Delaunay cells [33] where ξ min and ξ max are obtained from the eigenvalues of the pore orientation tensor (Eq. 1). The shape factor β ∈ [0, 1], where β = 0 represents a pore with regular geometry (e.g. cube, regular tetrahedron or sphere), while β = 1 implies that the pore is a line or plane (depending on the intermediate eigenvalue). The slight difference in the definition of the shape factor ensures consistency with the following pore anisotropy analysis, where a value of zero implies isotropy, and a non-zero value reflects some degree of anisotropy.

Pore orientation and induced anisotropy of the pore space
The directional distribution of a collection of pores at any instance during shearing can be quantified by the probability density function of pore orientation, E k (p), given by the following second-order polynomial expression: where D k i j is termed the deviatoric directional probability density tensor [17]. While this formulation is typically used to describe the directional distribution of contact normal vectors, it can also be applied to describing the distribution of individual pores.
In Eq. 3, the superscript k refers to a particular subset of pores. This concept of partitioned subsets was previously applied to the interparticle contact network [32], where dominant contributors to shear strength were identified. Sufian et al. [32] found a unique linear relation between the mobilised internal friction and the anisotropy of a partitioned subnetwork. The current paper considers whether a partitioned subset of pores can be identified as the dominant contributors to pore space anisotropy. In addition to considering the complete collection of pore units, this study will also investigate subsets of unmerged and merged pores, as delineated by the modified Delaunay tessellation and denoted k ∈ {u, m}, respectively. An absence of any value for k implies a quantity calculated for the complete collection of pores. Note that unmerged pores are representative of smaller pores in the assembly, while merged pores represent larger pores. Hence, the {u, m} partition also explores the influence of pore size on the anisotropy of the pore space.
Following Kanatani [9], D k i j in Eq. 3 is given by: where is the moment tensor of pore orientation or alternatively a pore space fabric tensor. N k p is the number of pores within the k th subset. A scalar measure of the degree of anisotropy can be defined by the second invariant of D k i j : Note that a signed measure of anisotropy is presented, where for axisymmetric simulations (Fig. 4a-c), positive values imply preferential alignment in the vertical direction, while negative values imply preferential orientation in the horizontal direction (Fig. 5a). For the simple shear simulation (Fig. 4d), negative values imply preferential orientation of pores tilted to the right (Fig. 5b), while positive values imply preferential orientation of pores tilted to the left (when observing the x-z plane).

Orthogonal preferential alignment of unmerged and merged pores
When considering the complete collection of pores, there is a preferential alignment in the direction of expansion. This implies that in the axisymmetric simulations, pores are orientated horizontally in the loading phase, before becoming vertical in the unloading phase and returning to a horizontal preference in the reloading phase. This is consistent with the observations in the simple shear simulation, where pores align to the right in the loading phase (which is the direction of tilt and the macroscopic direction of expansion), while they are orientated to the left in the unloading phase and align to the right in the reloading phase.
A more detailed perspective is provided by considering the anisotropy in the {u, m} partition. Unmerged pores are the dominant contributor to pore space anisotropy, with a u p > a m p and this is particularly noticeable in the CP-L (Fig. 4a) and CP-D (Fig. 4b) simulations. In contrast, the magnitude of pore anisotropy is significantly smaller in the constant volume shear simulations (CV-D in Fig. 4c and SS-D in Fig. 4d). These cases show a near linear trend in pore anisotropy due to the (relatively) small strain range considered in these simulations. Similar linear response in pore anisotropy was found in the CP-L and CP-D samples for shearing to |ε a | = 2% (not presented here).
In all simulations, unmerged pores exhibit strong preferential orientation in the expansion direction (perpendicular to the compression loading direction), while merged pores have a slight orientational preference parallel to the compression direction. This is schematically shown in Fig. 5 and was also observed in numerical simulations of bi-axial compression tests and direct shear tests [10,11].
This orthogonal preferential orientation of merged and unmerged pores can be explained by considering the interparticle contacts along the edges of a pore unit. Recall that a pore unit formed by the modified Delaunay tessellation has polyhedral geometry (Fig. 3). Some edges of the polyhedral unit cell may contain interparticle contacts, while other edges will have a separation distance between non-interacting particles. Suppose κ c = N c N e is the ratio of active contacts, where N c is the number of contacting edges and N e is the number of edges of the pore unit. Also of importance is the ratio of sliding contacts, κ s = N s N c , where N s is the number of contacts which are sliding (i.e., at the Coulomb friction limit). Consider also the magnitude of the interparticle normal force, f n of the active contacts. The average force magnitude of contacts forming a pore unit is defined as κ f = f n / f n , where the force is normalised by the mean normal force, f n , and the average is taken only over active contacts. The parameters κ c , κ s and κ f effectively capture the influence of the contact network on individual pores. Figure 6a shows that the ratio of active contacts (κ c ) and sliding contacts (κ s ) fluctuate around a constant value for pores of all size, but the average interparticle normal force increases with pore volume (Fig. 6b). Note that an approximate delineation between unmerged and merged pores is shown in Fig. 6. Even though merged pores contain the same proportion of active and sliding contacts as unmerged pores, it is clear that merged pores can only be sustained by interparticle contacts with large normal force (Fig. 6b). This is an intuitive result, as large pores would collapse if the surrounding interparticle contacts had relatively small normal force. It is well-known that contacts with relatively large normal force magnitudes tend to align parallel to the loading direction. Hence, if merged pores predominantly comprise strong contacts, then it would be anticipated that these pore units would also exhibit a preferential orientation in the direction of compressive loading, as indicated in Fig. 4. Based on biaxial compression tests on photoelastic disks, Oda et al. [23] had previously indicated that elongated pores were formed parallel to the major principal stress direction due to arching structures of strong force chains. This study suggests that the observation by Oda et al. [23] is only applicable for merged pores. Figure 7 shows that steady state pore anisotropy is attained at large strains (and hence only observed in the CP-L and CP-D simulations). This steady state pore anisotropy is attained even though the sample is macroscopically undergoing volume change, suggesting that there are additional factors contributing to the deformation response. Comparisons of pore anisotropy in the loose and dense assemblies suggests the attainment of a "critical state" for the a p parameter (i.e. for the complete collection of pores), although shearing to higher strains is required to confirm this observation. Additional analysis are required to establish uniqueness of this pore space parameter for shearing to large strain critical state conditions. This critical state characteristic is not found for the subsets of unmerged or merged pores (Fig. 4a, b). Li and Dafalias [18] proposed that there is a unique pore fabric parameter that characterises critical state conditions, along with the usual constraints of stress state and void ratio. This suggests that the current pore anisotropy parameter, a p , may be incorporated into such an approach, provided it is expressed in a thermodynamically consistent manner as per the suggestion of Li and Dafalias [19]. Figure 7 also indicates that the steady state pore anisotropy mobilised in the unloading phase (i.e. triaxial extension) is not equal and opposite to that observed in the loading and reloading phase (i.e. triaxial compression). The steady state value in compression is a p ss, comp ≈ − 0.27, while in extension it is a p ss, ext ≈ 0.46. Note that there is a natural asymmetry in the maximum attainable pore anisotropy in the (a) (b) Fig. 5 Schematic diagram showing the sign convention (shown by plus and minus symbols) used for the determination of pore space anisotropy in axisymmetric and simple shear simulations. The diagram illustrates the orthogonal preferential orientation smaller unmerged and larger merged pores, with unmerged pores aligning the expansion direction (as discussed in Sect. 3.1). It also illustrates that smaller unmerged pores are more rounded in shape compared to the elongated larger pores (as discussed in Sect. 4) radial and axial direction, due to the fact that alignment in the axial direction allows for only perfectly vertical pores, while alignment in the radial direction can have an infinite number of orientations in the horizontal plane. This asymmetry is also evident in the stress-strain response in Fig. 1 during triaxial shear, as there is one major principal axis in the axial direction during compression, but two major principal axes in the radial direction during extension. However, this significant asymmetry was not noted for the contact network anisotropy parameters [32]. For example, Fig. 8 shows that the steady state contact anisotropy is (a c ) ss, comp = 0.74 in compression and (a c ) ss, ext = − 0.76 in extension for the CP-D simulation. Contact anisotropy, a c , is also obtained by using the contact normal vectors in Eqs. 3-6 in place of the pore orientation vectors, thus enabling comparison with the pore anisotropy parameter, a p . This asymmetry in pore anisotropy, and the lack of it in the particle contacts, makes it difficult to correlate these two properties. Fu and Dafalias [5] used numerical simulations of (2D) bi-axial compression tests to suggest that pore space and contact network fabrics can be related (noting that they employ different definitions of pore (a) (b) Fig. 6 a Average percentage of active and sliding contacts for individual pores as a function of normalised free pore volume, where V p is the pore volume, V p,min is the minimum pore volume assuming no particle overlaps and d s is the particle diameter. This normalisation enables comparison between samples of different densities and particle diameter. Note that the averaging is conducted over all data in the initial loading phase for all simulations. b Average interparticle normal force as a function of pore volume. In both figures, the green line represents the approximate delineation between unmerged and merged pores and contact fabric). However, a direct correlation between a c and a p has not been found in this study. Fig. 8 Comparison of stress ratio, q p , contact network anisotropy, a c , for the complete contact network, and pore anisotropy, a p , for the complete collection of pores in the CP-D simulation. Note that the x-axis for two-way shearing has been unfolded into a single linear strain axis, and that a negative value of pore anisotropy is shown to allow for ease of comparison with the other parameters In addition to this asymmetry, comparison of the pore and contact anisotropy indicates that pore space anisotropy is substantially smaller than observed for the contact network with a p < 0.5, while |a c | < 1.0 (similar trends are noted for the partitioned subsets of pores and contacts). This was also noted in Oda et al. [23] (albeit with different measures of contact and pore space anisotropy) and reflects the fact that a highly anisotropic pore space (e.g. all pores aligned in a particular direction) would not form a mechanically stable assembly when subject to shear. Figure 8 also shows that the stress ratio response has a slight asymmetry ( q p = 0.71, 0.59 for critical state compression and extension conditions, respectively). However, a direct correlation with the stress ratio is hindered by the fact that pore space anisotropy a p evolves at a much slower rate than the contact network anisotropy (a c ) and the macroscopic stress response. This is most evident by considering the unloading and reloading phases in Fig. 8. In the unloading phase, the strain required to change loading direction (or anisotropy direction) is |δε a | = 0.8% for the stress ratio, while it is |δε a | = 2.5% for contact anisotropy and |δε a | = 4.1% for pore anisotropy. In the reloading phase, sign reversal requires |δε a | = 2.0, 7.1, 14.3% for stress ratio, contact anisotropy and pore anisotropy, respectively. This apparent lag in the response of the pore space reflects the fact that re-alignment of pores requires a substantial displacement of particles, which in a "dense" or slow-flowing granular assembly (with low inertia number) takes substantial amount of time. In contrast, the stress state can switch directions quickly by adjusting the magnitude of contact forces initially and then allowing the contacts to re-align as observed by Sufian et al. [32].

Isotropic directional distribution of pore volume in sheared assemblies
The observations in Sect. 3.1 for pore anisotropy suggests that parallel to the loading direction there are larger pores but fewer of them, while perpendicular to the loading direction, there are smaller pores but in greater number. This is schematically illustrated in Fig. 5. This raises an interesting hypothesis as to whether these two competing effects cancel each other out and result in an isotropic distribution of pore volume. To further investigate this hypothesis, consider the following distribution function: which represents the average pore volume in a given direction. This analysis is only conducted for the complete collection of pores and hence, the superscript k denoting partitioned subsets is omitted. T i j is given by: where is a moment tensor and also of the same form as the global orientation tensor introduced in Sufian et al. [31].
Note that M i j incorporates both competing effects of distribution by number and distribution by volume, and hence, differs from the representative value analysis presented in Li and Yu [17] which is typically used to explore directional distribution of normal and tangential force magnitude independent of any contact normal anisotropy [32]. Within the context of this paper, it is desirable to explore pore distribution by number and volume together.
It can be shown that V 0 = M kk = V p − V p,min , that is, the mean free pore volume, and in a similar manner to above, a scalar measure of anisotropy can be defined by: If the distribution in Eq. 7 is isotropic, then a t = 0 and the global orientation tensor M i j must be an isotropic tensor with diagonal elements given by . Sufian et al. [31] showed that the global orientation tensor was isotropic for a collection of static granular assemblies obtained experimentally and numerically. Figure 9a shows that the anisotropy parameter a t ≈ 0 in sheared granular assemblies, with some slight deviations away from perfect isotropy (as shown in the inset to Fig. 9a), but this deviation is two orders of magnitude less than the observed values for a p . Fig. 9b shows that the global orientation tensor is approximately isotropic with diagonal elements defined by the mean free pore volume. Note that this is not observed when considering the unmerged or merged pores separately. Therefore, deforming granular assemblies maintain isotropic pore volume distribution for all the simulations considered in this study. This suggest that when subject to external mechanical loading, particles adjust in a manner to maintain isotropic volume distribution, despite the induced anisotropy in pore orientation.

Pore shape as an indicator for the macroscopic density
Further insights into deformation characteristics can be gained by exploring pore shape properties. Figure 4a, b show that CP-L and CP-D specimens show similar pore anisotropy in loading, although the loose specimen undergoes contraction during shearing and the dense dilates (Fig. 1b). This difference in their volumetric behaviour can be captured by exploring pore shape. Pore shape is defined using the shape factor, β, in Eq. 2. Figure 10 plots the average shape factor throughout the twoway loading cycle. Both the complete collection of pores, and the {u, m} partition are considered.
When comparing the CP simulations, Fig. 10a shows that the loose sample (CP-L) initially comprises more elongated pores (i.e. higher β ), while the dense sample (CP-D) contains more rounded pores (i.e. lower β ), when averaging over the complete collection of pores. This was noted in Schröder-Turk et al. [28] with regards to Voronoi cell shape in static assemblies of spheres, where it was suggested that elongated Voronoi cells contain more free volume and hence, represents a loose assembly. In the {u, m} partition, Fig. 10b shows that smaller unmerged pores have more reg- Moreover, Fig. 10 shows that β (i.e. for the complete collection of pores) closely correlates with the void ratio response for the CP-L and CP-D simulation (Fig. 1b), suggesting that it is a good indicator of macroscopic density. This is confirmed in the results for the CV-D and SS-D simulations in Fig. 11, where β maintains a constant value, reflecting the constant volume constraint in these simulations. The constant shape factor is also observed in the {u, m} partition for the CV-D and SS-D simulations. This strong correlation between shape factor and void ratio is illustrated in Fig. 12 and suggests that (averaged) pore shape characteristics represent a micro-scale indicator of the tendency of the assembly to contract or dilate. Nevertheless, it is important to recognise that at the scale of individual pores, there is a distribution in pore shape as shown in the inset to Fig. 12. Relationship between average shape factor, β , and the macroscopic void ratio for all simulations. Note that a similarly strong correlation is noted for unmerged pores, β u , indicating the dominant contribution of unmerged pores. Inset: Distribution of pore shape factors in the CP-D simulation at a = 5% in the initial loading phase for the complete collection of pores (solid line), unmerged pores (dashed line) and merged pores (dotted line) The average shape factor β can be viewed as a measure of shape anisotropy, while a p provided a measure of orientational anisotropy of the granular assembly. Note that these two forms of anisotropy have opposing trends. The loose assembly (CP-L) exhibits a lower degree of orientational anisotropy (Fig. 4) but a higher degree of shape anisotropy compared to the dense assembly (CP-D). This suggests that deformation characteristics may be represented as a competition between shape and orientational anisotropy.
There is an interesting interplay between shape and orientation of pores when combining the results of the orientational analysis in Sect. 3 and the above shape analysis. Smaller unmerged pores are more rounded and align in the expansion direction, while larger merged pores are more elongated and align in the compression direction (schematically shown in Fig. 5). This is reflected in the directional averaged β in Fig. 13.
In the initial loading phase for the CP-D simulation (Fig. 13a), sub-horizontal pores (θ = 0 • − 30 • , perpendicular to the loading direction) are more rounded, while sub-vertical pores (θ = 60 • − 90 • ) are more elongated, and pores with θ = 30 • − 60 • fall in between. As the loading direction is changed in the unloading phase, sub-horizontal pores become more elongated (as they are now aligned parallel to the loading direction), while sub-vertical pores become more rounded (and so on for the reloading phase). This demonstrates that pore shape properties are closely tied to the orientational anisotropy of the pore space, which is confirmed in the CV-D simulation (Fig. 13b). Although the average shape factor was constant in the CV-D simulation, there is a clear directional dependence, and in fact the trends noted for the CP-D simulation are equally applicable for the CV-D simulation.

Conclusion
This paper provides insights into the deformation characteristics in sheared granular assemblies from consideration of the orientation and shape properties of the pore space. This was explored using discrete element simulations for common geotechnical elemental shear tests including drained and undrained axisymmetric and simple shear modes. Individual pores were identified using the modified Delaunay tessellation [1], which provides a physically intuitive representation of the pore space. The modified Delaunay tessellation was