Percolating contacts network and force chains during interface shear in granular media

The concept of force chains transmitting stress through granular materials is well established; however identification of individual force chains and the associated quantitative analysis is non-trivial. This paper proposes two algorithms to (1) find the network of percolating contacts that control the response of loaded granular media, and (2) decompose this network into the individual force chains that comprise it. The new framework is demonstrated considering data from discrete element method simulations of a ribbed interface moving against a granular sample. The subset of contacts in the material that transfers load across the sample, namely the percolating contact network (Gperc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G_{perc}$$\end{document}), is found using the maximum flow algorithm. The resulting network is fully-connected and its maximum flow value corresponds to the force percolating the system in the direction normal to the ribbed wall. Gperc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G_{perc}$$\end{document} re-orientates in response to the ribbed interface movement and transmits 85–95%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$95\%$$\end{document} of the stress, with only 40–65%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$65\%$$\end{document} of the contacts in the sample. Then, Gperc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G_{perc}$$\end{document} is split into individual force chains using a novel implementation of the widest path problem. Results show that denser materials with increased force-chain centrality promote a higher density of force chains, which results in a higher macro-scale strength during interface shearing. The contribution of force chains in the network is revealed to be highly centralized, composed by a small set of strong and long-lived force chains, plus a large set of weak and short-lived force chains.


Introduction
The response of granular media to applied shear deformation or shear stress is controlled by a interconnected network of load-bearing contacts (usually referred as the strong contacts network) [1]. Based on this observation, [1,2] proposed a partition of the contact network into strong/weak subnetworks, based on whether the magnitude of the contact force is above/below the mean force in the system. Developing upon this simple approach to partitioning authors have argued that the overall contact network is formed by 'primary' load-carrying contacts, and 'secondary' contacts that mostly provide lateral support [1,[3][4][5][6][7]. Various studies have found that the distribution of contact forces follows a power law with few contacts carrying high forces, and most contacts carrying relatively low forces [1,6,[8][9][10][11], and this distribution is independent of pressure, sample size and particle size distribution [1,12]. However, [13] found that this approach to partitioning is not mechanically robust, and using the mean force to differentiate strong and weak contacts is not always appropriate.
Recent studies have used graph theory, representing particles and contacts as nodes and edges [14], to propose partitioning methods based on the topology of the contacts. For instance, [15][16][17] used persistent homology to study tapped granular media, the compression of soft granular matter, and the stick-slip events in granular media (respectively). Community detection and clustering have been used to study granular materials in the context of sound transmission [18], deformation localisation [19,20], kinematics [21] among many others.
Flow networks (i.e. maximum flow and minimum-cut algorithms) have also been used to study the mechanical behaviour of granular media. [22,23] applied the MFA to biaxially/triaxially loaded granular media, and found that the set of edges forming the bottleneck for force transmission in the system is located in the region of strain localisation (i.e., shear bands). Similarly, [24] developed a family of network flow models and used the MFA to find the percolating sub-network of contacts that transmits the highest units of force at least dissipated energy. These studies used expressions based on the local topology (3-cycle) and/or the magnitude of the relative displacement of the particles involved in the contact (edge) to set the capacity and/or cost of flow through each edge; effectively using the maximum flow as a relative measure of the force that could be transmitted through the contact network. Here we propose an algorithm that extends this capability, using the MFA on a directed network where the capacity of the edges is equal to the magnitude of the contact force in the direction of load transfer.
Other authors have found that the mechanical behavior of granular systems is controlled by the collective contribution of force chains, i.e. sets of particles that form 'column-like' structures that carry load [3,5]. Methods used to identify force chains include 'geometric-mechanic' methods, which calculate the major principal stress orientation of the set of particles which is under above-average stress, and then defines force chains as the groups of particles which align with the principal stress orientation [25,26]. Other 'similarity' methods use community detection, a clustering technique capable of finding chain-like structures in graphs [27]. Having identified individual force chains, prior researchers have studied force chain buckling [10], quantified their contribution to deviatoric stress [6], and tracked their temporal evolution [28][29][30][31][32].We propose an algorithm based on the widest-path problem to identify individual force chains from G perc that collectively capture the force transmission across the granular media, and ultimately its mechanical response.
In this contribution we apply our developed algorithms to study the behaviour of granular media in contact with a rigid ribbed interface that applies a shear deformation in the direction of the initial major principal stress. This geometry and loading scenario has been used to measure soil properties in the laboratory [33] and for in-situ characterisation (using textured cone penetrometer sleeves) [34], and to study the seismic response of granular media inside gouges [35,36]. These studies have found that the response of the interaction is controlled by the stress conditions, the density of the sample, and the relationship between particle size and the size of the ribs in the rigid body. Based on these observations, the dataset was developed using DEM to test the influence of different stress conditions and with different grain sizes.

DEM model
The dataset considered in this study was developed by simulating a ribbed interface moving against a granular material at a constant velocity. The texture morphology illustrated in Fig. 1d was developed following work by [37][38][39][40], In the current context, the role of the data generated are to illustrate the new methods developed herein.
The DEM models were built using the molecular dynamics code LAMMPS [41]. The simplified Hertz-Mindlin contact model was used [42], the particle shear modulus G = 29 GPa , the Poisson's ratio = 0.12 , and the particle density = 2600 kg m −3 , are consistent with the properties of quartz. The use of spherical particles enabled large systems to be studied and minimize boundary effects. We acknowledge that excessive rotation of spheres does lead to a lower shear strength in comparison with sands. However, simulations using spheres capture key aspects of soil behaviour, e.g. the state dependency of strength and dilation, phase transformation, etc., as discussed in [43,44].
To showcase the developed algorithms a parametric study was designed to consider a wide range of stress conditions and particle sizes to assess their influence in the response of the system. A total of 27 ( 3 3 ) simulations were performed, combining three particle size distributions, three vertical stresses ( � v ∶ 50, 100, 150 kPa ), and three particle friction coefficients during preparation ( p ∶ 0.01, 0.1, 0.25 ). The preparation particle friction coefficient ( p ) determined the density and initial anisotropy of the samples. The particle size distributions (PSD) were scaled (scaling factor k = 1, 2, 4 ) from the PSD of Toyoura sand, a fine rounded sand with a mean particle size d 50 = 0.22 mm , coefficient of curvature C c = 0.96 and coefficient of uniformity C u = 1.39 [45]. The resulting PSD's, referred herein as PSD std , PSD mid and PSD lrg , have median particle diameters ( d 50 ) of 0.22, 0.44, and 0.88 mm respectively.
The sample creation process is illustrated in Fig. 1. The width of the model ( W m ) was set to: W m = h r + 50 ⋅ d 50 , while the depth was set to D m = 15 ⋅ d 50 , where d 50 is the median particle diameter in the samples. The sample size ( H m × W m × D m ) was achieved with copies of a random cubic sample of side 15 × d 50 (see Fig. 1c). Random sphere packings yield relatively loose samples, therefore, an extra densification step was included to reduce the void ratio of the initial cubic sample, and with it, reduce the computational cost required to reach the desired confinement stress. In the first step, the spherical particles were placed randomly to create a cubic sample with a void ratio e = 1.0 (see Fig. 1a). In this initial stage, the diameter of each particle was sampled from the target PSD, and reduced by a factor of 1.2 prior to its placement. In the second step, a simulation was set up to incrementally increase the diameter of the particles in the cubic sample until they equalled the original diameters sampled from the PSD (see Fig. 1b). After this step, the cubic sample had a new void ratio e = 0.67 . In the third step, copies of the cubic sample were tiled to reach the desired sample dimensions (see Fig. 1c). Lastly, the sample and the ribbed wall were merged, and any particles overlapping or penetrating the ribbed interface (described below) were removed (see Fig. 1d). A summary of the dimensions of the DEM models are shown in Table 1 where the values of H m correspond to the densest scenarios ( p = 0.01 , � v = 150 kPa). Planar periodic boundaries were inserted at the top and bottom of the sample (i.e. at the minimum and maximum z values) and at the front and back of the sample (i.e. at the minimum and maximum y values). The left boundary wall (parallel to the yz plane at x = 0 ) was a planar rigid wall. A rigid, textured interface was inserted at the right hand side of the sample (i.e. the maximum values of x). The ribbed interface comprised wall particles with a diameter d w = d min ∕2 , where d min is the minimum particle diameter the simulated sample. These particles were placed on a square grid with a centre-to-centre spacing of d w and there were no bonds or interactions between them. The heights of the trapezoidal ribs ( h r = 5 mm ) correspond to about 23, 11 and 6 times d 50 for each PSD tested; these values are sufficient to generate texture clogging and induce particleto-particle shear response [46]. The ribbed geometry follows [34,46], comprising alternating trapezoidal ribs with a height h r = 5 mm and flat sections, with a length of 9 mm each. The wall includes two rib-flat sections, for a height of 36 mm , plus an additional flat region (added to account for sample compression during confinement) of 4 mm on top for a total sample height H m = 40 mm . The number of wall and sample particles for the different PSDs tested are shown in Table 1.
After sample preparation, the models were allowed to reach equilibrium to remove any initial overlap between particles, after which residual stresses were negligible. Then, a servocontrol mechanism was used to confine the sample particles by adjusting the position of the top periodic boundary (parallel to the z-direction) until the target vertical effective stress value ( ′ v ) was reached. The dimensions of the samples along the  x-and y-directions remained unchanged during the simulation-therefore, the horizontal (x-axis) stress in the samples is consistent with K 0 stress conditions (see Fig. 1d). Gravity was not considered in the simulations. The displacement of the vertical boundary did not affect the wall particles that make up the textured rigid wall. The vertical strains induced during this stress-controlled confinement, zz , ranged between 4.53% (for PSD std , p = 0.25 , � v = 50 kPa ), and 13.48% (for PSD mid , . Once the target ′ v was reached, the particle friction coefficient was changed to s = 0.25 , and kept constant during wall movement. To simulate interface shear, the ribbed interface was displaced vertically (i.e. in the z-direction, see Fig. 3a) at a velocity V z = 25 mms −1 . During the interface shear the wall particles behaved as a rigid body moving at constant velocity. The chosen velocity is similar to the penetration speed of cone penetrometers [37]. Following [47], the inertial number (I) was calculated as shown in Eq. 1.
where ̇ is the strain rate, d 50 is the median particle diameter of the material, is the particle density ( 2600 kg m −3 ), and P is the stress level in the material, taken as ′ v . The strain rate is calculated according to [48,49] as: ̇= V z ∕L p ; where L p is the width of the shear zone developed adjacent to the textured wall. The shear zone thickness was calculated using the bi-linear relationship between L p (normalized with by d 50 ) and the d 50 defined by [50]. Following that relationship together with observations of the extent of the shear zone in the simulations, conservative values of L p = 9, 7.5, 4 × d 50 were used for PSD std , PSD mid and PSD lrg respectively. Obtained values of I range between 3.7 × 10 −4 (for PSD std and � v = 150 kPa ) and 1.4 × 10 −3 (for PSD lrg and � v = 50 kPa ). The total displacement was 36 mm , equal to the total height of the ribbed wall.

Stress-deformation response
The average stress in the granular assembly during shear deformation is calculated from the contact forces and particle velocities according to [51] as: where V is the sample volume (excluding the volume of the textured wall), N c , N p are the total number of contacts and particles (respectively), f i is the ith component of the contact force vector, l j is the jth component of the branch vector that joints the centers of the particles forming contact c, and v i is the ith component of the particle velocity. According to CPT interface friction simulations from [49], quasi-static conditions are expected or I < 1 × 10 −2 , as was the case here, however the dynamic component of the stress tensor (Eq. 2) was included in the calculations as a conservative measure. The magnitude of the dynamic component of stress remained < 5% of the static component. The stress anisotropy in the samples is quantified with the lateral earth pressure coefficient where ′ kk is the normal component of stress in the direction of the k axis, calculated from Eq. 2. The inter-particle friction coefficient during initial compression ( p ) determined both the void ratio and initial stress anisotropy of the samples, with higher values of p promoting higher void ratios (lower packing density)-ranging between 0.35 and 0.40, and increased stress anisotropy following the initial sample consolidation, with K 0 values between 0.65 and 0.91. The average values of e and K 0 for different values of p are summarised in Table 2.
The overall stress response during shear deformation, shown in Fig. 2, is quantified in terms of the mobilised stress ratio q eq ∕p � , where the mean effective stress p ′ is calculated . The particle friction coefficient during preparation p , is the variable with the biggest influence in the stress mobilisation in the system, as it controls the density of the material (in terms of e), and the stress in the x-direction ′ xx (and hence K 0 ). Higher mobilised stress ratios ( q eq ∕p � ) were obtained for denser samples and higher values of K 0 . Moreover, these dense, anisotropic samples reach a steady state at significantly larger wall displacements; arguably even longer than the 36 mm of displacement tested in the simulations. These results are in agreement with the sleeve friction during CPT testing [52][53][54][55], and with results from interface friction with ribbed walls of similar geometry to this study [38][39][40].
The small increase in q eq ∕p � (and its fluctuations) with PSD, is a phenomenon that has been attributed to particle size effects, and to changes in the relative roughness ( R r = h r ∕d 50 ) of the ribbed wall. Butlanska et al. [56] modeled CPT tests with scaled PSD and found that the increase in the magnitude of the fluctuations of the soil response is related to the decrease in the number of particles in contact with the rigid surface. While [38] studied ribbed walls with the same geometry, and suggested that for higher R r values the relatively small particles clog and mobilise the particle-particle shearing resistance,

Percolating contacts network
In contrast to direct shear laboratory tests where the major principal stress is normal to the shearing direction, for the dataset considered here the shearing direction (z-axis) is parallel to the major principal stress orientation. Previous studies considering similar loading conditions [34,37,46] have shown that the response of the system is controlled by the force chains normal to the shearing direction, i.e. the x-direction in the current study. In this section, we propose a new method to partition the contact network into two complimentary sub-networks: (1) a principal network that transmits the force percolating through the sample in the x-direction, namely the percolating network G perc , and (2) a secondary contact network formed by the remaining particles and contacts that support the structure of the material Fig. 2 Mobilised stress ratio q eq ∕p � as a function of wall displacement for the three particle size distributions considered. The stress data were calculated considering the sample particles (i.e. excluding the wall particles) Fig. 3 Example model, percolating network and strongest force chains for a test with PSD mid , p = 0.25 and � v = 50kPa . a Complete sample, b particles that belong to the percolating network G perc , and c five strongest force chains in G perc but do not contribute directly to the force transmission in the x-direction, namely the supporting network G supp .
The partition of the contact network is achieved using network analysis. Particles and contacts in the samples are represented as nodes and edges (respectively) in a directed graph G, as described by [14]. The node (particle) information includes the particle diameter, position (x,y,z coordinates), coordination number (equivalent to the node degree in graph theory) and type (i.e wall or sample particle). The edge (contact) information includes the IDs of the particles forming the contact, and the force components (normal and tangential). The edge direction is set such that the contacts are oriented from the ribbed wall towards the planar rigid wall at x = 0 (i.e. the left wall see Fig. 3a which shows the complete sample) as: where x i is the x-coordinate of node i, and the contact c(i, j) is oriented from node i to node j.
The percolating network G perc is found using the maximum flow algorithm (MFA). Finding the maximum flow between a source s and a target t node in a graph with edges of finite capacity is a classical problem in optimisation theory [57]. The solution involves finding the set of paths connecting s, t which form the bottleneck of flow in the network. Adopting this approach ensures that: (1)  The response of the system is controlled by the force transmission in the x-direction, therefore, the capacity C of edge in the algorithm was set to the x-component of its contact force, namely C =| f n x + f t x | , where f n x and f n x are the x-components of the normal and tangential forces. In this way, the flow through corresponds to the percolating force transmitted by that given contact, and is ≤ C . Virtual nodes were added to represent the ribbed interface (source s) and the opposite left wall at x = 0 (target t). Next, artificial edges of infinite capacity were added from s, to every wall particle, and from t to every particle in contact with the left wall. In this way, the algorithm captures all the paths of force transition in the x-direction between the left wall and the ribbed interface. Then, MATLAB's [58] implementation of the Edmonds-Karp's maximum flow algorithm [57] was used to find the maximum flow network and its maximum flow value.
The maximum flow value corresponds to the total percolating force in the sample in the x-direction, while the subset of edges (contacts) with non-zero flow (i.e. contacts that contribute to the percolating force) form the percolating network G perc . The resulting G perc is fully connected, meaning that for every pair of nodes (particles) in G perc , a continuous path with non-zero percolating force exists. The mass-conservation principle of the algorithm ensures that the sum of the in-coming and out-going percolating forces on each particle (node) in G perc are equal.

Characteristics of particles and contacts in the sub-networks
The characteristics of the sub-networks ( G perc and G supp ) are analysed for the interval X p = [0, 36] mm of wall displacement, with only minor changes to the characteristics of the sub-networks during the interval. In the following, representative characteristics of the system are shown for X p = 36 mm . Figure 3b shows the particles in the percolating network G perc for a representative sample. Figure 4 summarises the proportion of particles, N p , and contacts, N c , in G perc as a function of PSD and p . G perc contains between 40 and 60% of the total number of contacts in the network, and involves between 60 and 85% of the particles in the samples. The proportion of particles and contacts in G perc increases with sample density (i.e. decreasing p ), and with increasing particle size. This observation suggests that denser materials, with larger particles create more pathways to percolate force in the system, this observation is verified later in the study of the force chains in the system (see Fig. 11). Figure 5 shows the cumulative distribution functions (CDF) of the particle sizes, the contact forces f c , and (the ratio between the tangential and normal contact force) in G perc and G supp . The particle sizes and contact forces are normalised as percentiles for easier comparison between simulations, while is shown in the range [0, s ] (where s = 0.25 is the inter-particle friction coefficient during shearing, equal for all simulations). Lines in Fig. 5 correspond to the mean CDF for all the simulations, while the interquartile range-IQR (i.e. the range between the 25th and 75th percentiles of the data) is shown as shaded regions. Due to the small differences between the CDF's in Fig. 5, the region enclosed by IQR is virtually indistinguishable from the mean lines.
The results show that G perc is formed by larger diameter particles and contact forces than G supp . However, both subnetworks contain particles of every size and contacts forces of every strength, rather than exclusively small/big particles or weak/strong contact forces, as suggested by a-priori partition techniques based on contact force [9,22,59]. The distribution of shows that, on average, contacts in G supp are closer to sliding than those in G perc , with 38% sliding contacts in G supp and 26% sliding contacts in G perc ; similar results have been observed by [1,59].

Stress transmission in the sub-networks
The stress at steady state is measured from the stress tensor of the contact network at 36 mm of wall displacement, as before. The stress tensor of each sub-network, i.e. ′ perc and ′ supp is calculated according to Eq. 2, including only the particles and contacts in G perc and G supp respectively. For each tensor, three components of stress: the 2-norm of the stress tensor ‖ ′ ‖ 2 (equivalent to the major principal stress ′ 1 ), the isotropic stress ( p ′ ), and the deviatoric stress ( q eq ) components are calculated. Figure 6 shows the distribution of each component of stress in G perc relative to G. Results show that G perc carries most (80-95% ) of the stress in the material, while having only 40-60% of the contacts of the network. Moreover, the fraction of q eq transmitted by G perc appears to be slightly higher than the fraction of p ′ this partition transmits, although the difference is not significant. Previous studies on element tests based on strong and weak sub-networks [1,5,7], have also found (expectedly) that strong contacts (above the mean contact force in the system) carry most stress in the material.
Representative rose diagrams are presented on Fig. 7 for at simulation at 36 mm of wall displacement (similar rose diagrams were obtained for all simulations). The rose diagrams show significant differences in the preferential contact orientation amongst the sub-networks. G perc has the strongest contacts in the network, and exhibits a preferred orientation between 35 • and 55 • (measured from the horizontal), which indicates a re-orientation of the contacts during wall movement. Conversely, contacts G supp are weaker in comparison, and their preferred orientation is closer to 90 • , i.e. close to the the vertical direction, consistent with the direction of the application of ′ v . The projection in the XZ-plane of the orientation of the major principal stress ( 1 ), is calculated as 1 = tan −1 ( z ∕ x ) , where x and z are the x and z components of , the principal direction (eigenvector) associated with the major principal stress ′ 1 . 1 is calculated for G, G perc and G supp , and presented in Fig. 8, and as arrows in Fig. 7. The results indicate that the displacement of the ribbed wall induces a rotation of the principal stresses in the material from its initial K o state. This reorientation is reflected in the major principal stress orientation in G perc , while G supp remains oriented closer to the vertical direction, consistent with the rose diagrams shown on Fig. 7. Comparable results were encountered by [24,59] during uniaxial loading, where the strong network re orientates in the direction of loading, while the weak network remains mostly unchanged.

Contact and particle importance
The relative importance of the particles and contacts in G perc is measured in terms of node/edge centrality. The contact/edge centrality ( C cen ) is the percentage of the total percolating force transmitted by each contact in G perc , while the particle/node centrality ( P cen ) is the percentage of the total percolating force transmitted through each particle. Centrality is grouped by percentiles, from the most (1st percentile) to the least (100th percentile) important groups of particles/contacts. Figure 9 shows the distributions of C cen (the mean contact C cen ), and (the mean at the contacts). The C cen follows a exponential distribution, with few 'highly important' contacts. The values in Fig. 9a show that in dense samples ( p = 0.01 ), each contact in G perc transmits a higher percentage of the total percolating force when compared to less dense samples. The distribution of in in Fig. 9b indicates that contacts with higher centrality are, on average, further away from sliding.
Similarly, the distribution of the mean P cen , mean normalised particle size ( d∕d 50 ) and mean coordination number C N of the particles in G perc are presented on Fig. 10a, b and c respectively. The results show no significant differences between samples, with consistent exponential distributions of P cen , i.e. few 'highly important' particles. These highly important particles are in average 10% larger than d 50 , and have higher coordination numbers (between 6.5 and 7.0). The relation between particle importance with d∕d 50 and C N follows a sigmoid, logistic function.

Micro-scale: force chains
The percolating network G perc is the complete set of paths transmitting force across boundaries of the model, and thus can be decomposed into individual force chains, each of which transmits a unique value of percolating force. Each force chain is a simple path (i.e. a path without loops) of particles (and contacts between them) between the sink/ source, which transmits a given force between the ribbed wall, and the opposite boundary, i.e. the force chain's percolation force. Due to the orientation of the contacts in G (and subsequently in G perc ), particles in each force chain have an ordered set of x coordinates (see Eq. 3). Figure 3c shows an example of five force chains extracted from G perc .
Force chains are found using our implementation of the widest path problem based on Dijkstra's algorithm, which is traditionally used to find shortest paths in graphs. The widest path, also known as 'maximin' path, is defined as the path between a sink and source in the graph, which maximises the weight of the minimum-weight among the edges in such path [6]. In our context, the widest path is the force chain that transfers most percolating force between the ribbed wall and the opposite wall.
The pseudo-code of the proposed algorithm is shown in Alg. 1. the computation when only 10% of the initial percolating force remains in the percolating network. The pseudo-code of the algorithm that splits the percolating network into individual force chains is shown in Alg. 2. Fig. 10 a Mean particle centrality P cen , b mean normalised particle size D∕D 50 , and c mean coordination number C N of the particles in G supp , as a function of percentiles of particle centrality. Relative particle importance calculated as percentiles of particle centrality. Lines show the mean value across simulations and shaded regions show the inter-quartile range (IQR), namely the range between the 25th and 75th percentiles of the data from the 27 simulations In order to find all the force chains in the system, Alg. 1 is applied to G perc sequentially. Once a force chain is found, it is recorded, and it's percolating force removed from the contacts along its path, updating the graph. Due to the significant computational cost of the algorithm, we stop The distribution of the percolating force transmitted by the force-chains is highly centralized, with the top 10% and 20% of the force chains carrying 35% and 50% of the total force percolating the sample. The overall distribution of the percolating force follows an exponential trend (see Fig. 11a), similar to the distributions of particle importance (in Fig. 10a) and contact importance (in Fig. 9a).
The overall system centrality (also called the concentration or inequality) is measured using the unitary Gini coefficient g c . A g c value of zero corresponds to a system where all

Distribution of force chains
The characterisation of the force chains in the system is done for the interval between [33,36] mm of wall displacement (a window of 0.12 s ). The interval is sampled every 1d − 4 s , for a total of 1200 sampling points. The characteristics of the force chains were constant over time, with negligible differences between different points within the interval. The subsequent analysis on this section is based on the mean value of the quantities over the sampling interval. the force chains carry the same force (i.e. homogeneous system), and the value grows as the difference between forces increase (i.e. increased system centrality). Figure 11b summarizes the distribution of g c as a function of PSD, variable that showed the largest statistical influence. Results indicate that samples with larger particles 'depend' more heavily on a relatively few number of strong force chains, whereas similar samples with smaller particles distribute the load among force chains more evenly.
The 'density' of force chains D fc in the system is calculated as: where N pw is the average number of particles in contact with the ribbed wall, and N fc is the number of force chains that account for 90% of the total percolating force in the system (i.e., the maximum flow value transmitted by G perc ). Figure 11c shows the distribution of D fc for the 27 simulations as a function of PSD and p . Figure 11d shows a direct relationship between the density of force chains per unit area D fc and the average coordination number C N . This relationship is expected, as denser arrays have more contacts between particles, facilitating the creation of more pathways for force transmission. However, the range of C N does not vary with particle size, indicating that the increase in D fc with increasing PSD is not due to higher coordination numbers in the system, but instead, due a higher fraction of particles/contacts in G perc , in agreement with Fig. 4. (4) These results indicate that sample density, and K o increase the density of force chains in the system, while making it more centralised (in terms of g c ). These findings agree with the trends observed in Fig. 2, linking the increase in the mobilised stress in the material with an increased force chain density and system centrality.

Life-span of force chains
The longevity or duration of force chains during the induced shear deformation, a matter also addressed by [28], is considered next. We measure the duration of each force chain in the system within the same interval: [33,36] mm of wall displacement ( 0.12 s ), sampled every 1d − 4 s , for a total of 1200 points.
The duration of a force chain in the contact network is equal to the sample interval ( 1d − 4 s ) times the number of consecutive intervals where the same force chain can be found. Each force chain fc m is identified by as the set of contacts that form it, namely fc m ∶ c i,j , c j,k , ... , where c i,j is the contact between particles i and j. The similarity index, a modified version of the overlap coefficient [60], is used as the criteria to compare force chains at different times: The similarity index of force chains fc m and fc n , namely Sim (m,n) , is calculated as the number of contacts common to both chains, divided by the average number of contacts in fc m and fc n , as shown in Eq. 5. If two force chains in consecutive time steps (say fc m and fc n ) have Sim (m,n) ≤ 75% , we consider them to be the same force chain, and increase its duration. In the vicinity of the start/end of the time window studied the algorithm may underestimate the duration of the force chains.
Once the duration of the force chains in the system was calculated, we studied the distribution of the force-chains' percolating force for three groups of different duration: the baseline group contains all the force chains in the system regardless of their duration, namely 'All'; the second group includes only the force-chains with a duration above 2% of the analysed time interval, namely t fc > 2%t T ; similarly, the third group includes the force chains with a duration above 5% of the analysed time window, namely t fc > 5%t T .
The results in Fig. 12 show that force chains with longer life-spans, i.e. long-lived force chains, are (in average) stronger than short-lived chains. These results are in agreement with Fig. 9, as strong, 'important' contacts in the longlived force chains have lower values of (less prone to sliding), explaining their increased duration.
Moreover, the comparison between different PSDs in Fig. 12 shows that in samples with PSD lrg , long-lived chains are comparatively stronger that force chains of the same duration with smaller PSD. This observation is supported by the increased centrality of systems with larger particles shown in Fig. 11b.

Conclusions
This contribution has proposed new algorithms to: (1) extract the fully-connected set of percolating contacts ( G perc ) that controls the stress/deformation response of granular systems, and to (2) decompose G perc into the underlying set of force chains that form it. These algorithms were then applied to a particular scenario, the response of granular materials to shear deformation induced by a ribbed interface moving at a constant speed and the steady-state response of the material was considered here. The data show that granular materials behave as highly centralised systems, where the mechanical response to shear deformation is ultimately controlled by a small set of strong and long-lived force chains.
In Sect. 4, we proposed a partitioning method that splits the contact network into a stress-carrying percolating network G perc and a supporting sub-network G supp . This approach is not the first to use the maximum flow algorithm, but (to the best of our knowledge) is the first to use directed edges with a physically based capacity (transmitted force). By doing so, G perc captures the force percolating between boundaries, and all the paths that transmit it. Thus, this partitioning method has a stronger fundamental basis than the commonly used approach that splits the force network by considering the mean force.
Results show that G perc contains 35-65% of the contacts (see Fig. 4), but carries most of the stress in the sample ( > 80% , see Fig. 6). Interestingly, it was found that G perc has a larger proportion of large particles, but both sub-networks include particles across the entire PSD. Similarly, the contact-force distribution shows that G perc has stronger contactforces in average, but both sub-networks contain weak to strong contacts. The major principal stress orientation shows that G perc re-orientates at an angle between 35 • -55 • as a response to the ribbed wall movement, while G supp carries little stress and remains oriented along the vertical direction.
In Sect. 5, we proposed a method to split G perc into individual force chains based on a new implementation of the widest path problem. With this method, each resulting force chain is an aligned set of particles/contacts in G perc that transmits a specific value of percolating force.
The force chains in the system form a highly centralised system that follows an exponential distribution, with a small number of strong and long-lived force chains, and a large number of weaker, ephemeral chains. These observations have been related to the stress-deformation response of the system, showing that systems with a more centralised distribution of contacts and higher density of force chains are able to mobilise higher mobilised stress ratios during shear deformation, providing micro-mechanical insight to previous laboratory and in-situ testing observations.
The methods discussed here may be applicable to a wide range of applications in geomechanics, including boundary problems, monotonic/cyclic loading and even pore network modeling. Future studies may also focus on the role of polydispersity and particle shape in the centrality and density of the force chains in the material. Understanding the role of the particle-scale properties in the behavior of force chains may allow to design granular materials that maximise the performance or resiliency of the system, i.e. the capacity to accommodate changes in loading into the existent network of force chains, under specific loading cases. Data availability The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.