# Spanning the scales of granular materials through microscopic force imaging

Nicolas Brodu, Joshua Dijksman, Robert Behringer

► To cite this version:

Nicolas Brodu, Joshua Dijksman, Robert Behringer. Spanning the scales of granular materials through microscopic force imaging. Nature Communications, 2015, 6, pp.6361. 10.1038/ncomms7361. hal-01128576

HAL Id: hal-01128576

https://hal.science/hal-01128576

Submitted on 10 Mar 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE

Received 4 Nov 2014 | Accepted 23 Jan 2015 | Published 5 Mar 2015

DOI: 10.1038/ncomms7361

OPEN

# Spanning the scales of granular materials through microscopic force imaging

Nicolas Brodu<sup>1,2,*</sup>, Joshua A. Dijksman<sup>1,3,*</sup> & Robert P. Behringer<sup>1</sup>

If you walk on sand, it supports your weight. How do the disordered forces between particles in sand organize, to keep you from sinking? This simple question is surprisingly difficult to answer experimentally: measuring forces in three dimensions, between deeply buried grains, is challenging. Here we describe experiments in which we have succeeded in measuring forces inside a granular packing subject to controlled deformations. We connect the measured micro-scale forces to the macro-scale packing force response with an averaging, mean field calculation. This calculation explains how the combination of packing structure and contact deformations produce the observed nontrivial mechanical response of the packing, revealing a surprising microscopic particle deformation enhancement mechanism.

<footnote n="1">¹Department of Physics, Duke University, Physics Building, Science Drive, Box 90305, Durham, North Carolina 27708, USA. ²Institut National de Recherche en Informatique et en Automatique, Bordeaux Sud-Ouest, 200 avenue de la Vieille Tour, 33405 Talence, France. ³Laboratory of Physical Chemistry and Colloid Science, Wageningen University, PO Box 8038, 6700EK Wageningen, The Netherlands. * These authors contributed equally to this work. Correspondence and requests for materials should be addressed to R.P.B. (email: bob@phy.duke.edu).</footnote>
The way that a disordered packing of macroscopic particles responds to mechanical deformations, such as compression, was investigated as early as the seventeenth century, when Stephen Hales studied the expansion of dried peas submerged in water¹. The motivation for such studies is clear. Packings of particles surround us: coffee beans and rice; soils, embankments, many industrial processes and geophysical processes from earthquakes to landslides involve granular materials², hence, the gamut of recent investigations³⁻⁵. Part of the relevant physics also applies to emulsion and foams, which consist of highly deformable frictionless particles, and to colloidal systems⁶,⁷, where thermal agitation is relevant. A major challenge for granular materials is relating microscopic contact and force details to global system response, such as applied stresses or strains. The link from micro- to macro-scale behaviour informs the mathematical link, known as a constitutive relation, between macroscopic strain or stress and macroscopic response. These relations appear routinely for conventional materials such as Newtonian fluids (shear stress proportional to shear strain rate) or linear elastic solids (stress proportional to strain). For granular materials, there are continuum models containing empirical constitutive relations⁸, which are widely used in the soil mechanics communities. However, they suffer from problems: the connection between grain-scale and macro-scale behaviour is still an open question, and these models are known to contain complex mathematical instabilities⁹. The lack of a well-established constitutive relation for granular materials is highly limiting; it is essential in applications where a continuum description is the only computationally feasible approach.

To develop models such as constitutive relations, it is essential to have direct access to microscopic experimental information that shows how granular materials respond to applied stresses or strains. The challenge is to have access to all microscopic

quantities inside a granular packing, including at every contact between grains, not just for a single packing under fixed boundary conditions, but during the entire evolution of a packing under realistic deformations. Experimental probes of granular microscopic properties that provide structure and contact forces in three-dimensional (3D) packings are limited. Photoelastic methods provide this type of information in two dimensions¹⁰, but are hard to implement in 3D. Pioneering measurements by Brujić et al.¹¹ and by Zhou et al.¹² for emulsions yielded the distribution of forces at contacts, $P(f)$, relative to the average force and information on contact networks. Tomographic techniques (such as X-ray microscopic computed tomography (microCT)¹³,¹⁴ and refractive index matching¹⁵,¹⁶) have been successfully implemented to measure structure and particle contacts. Extracting forces was done in microCT data¹⁷,¹⁸, by imposing a global force balance constraint. This constraint infers local contact forces from (measured or simulated) average particle stresses and boundary conditions, even when contact forces cannot be easily extracted. An alternative approach would be to measure forces directly, at each individual contact. Both the methods could be combined to improve accuracy. Importantly, all existing methods are limited in their ability to yield microscopic structure, particularly contact forces, over many closely spaced macroscopic state changes and with a precise control of applied strains.

In this work, we use refractive index matching tomography to provide full access to microstructure, for packings of deformable hydrogel particles that are compressed and decompressed uniaxially. The crucial aspect of our work is the unique combined ability to measure individual contact forces in vectorial detail, while straining the sample in small increments, enabling us to track the system-scale stress tensor over many small strain steps. This feature gives access to the complete micro-macro range of mechanical details of the packing, including ingredients for constitutive modelling.

## Results

**Description of the experiment.** In a typical experiment, sketched in Fig. 1a, 514 hydrogel particles that have been saturated with fluorescent dye are contained in a Plexiglas box. The particles are roughly, but not perfectly, spherical, with a mean diameter of 2.1 cm. The box is filled with a solution of water and polyvinylpyrrolidone, such that the index of refraction of the particles is well matched to the solution. This allows optical access to a vertical laser sheet that is scanned horizontally. A camera, whose image plane is parallel to the laser sheet, records a series of images as the laser sheet is swept. To reconstruct particles and contacts, we have developed dedicated tomography algorithms that yield the force (including direction) for each individual particle contact. We discuss this method further, along with additional experimental information, in the Methods section. The particles are contained in a box with five rigid transparent walls. The sixth and topmost wall is a porous piston, which moves in the vertical direction, providing compression/decompression. A force gauge, which is in-line with the piston, measures the vertical force acting on the top layer of particles, and hence the pressure at the top boundary. The piston location and the force gauge give us the macroscopic strain and stress imposed on the system. For each scanned sequence of images, we infer the particle shapes, locations and all the contact vector forces, which give us the full microstructure (Fig. 1b). For the present system, friction is negligible and forces are nearly normal to contacts. However, the technique used here can be generalized to contact forces with friction. The density of the particles is also nearly matched to the fluid, so that the particles experience an effective gravity of about 0.01 g (g is the acceleration of gravity).

<figure><img src="image_2.png" /><figcaption>**Figure 1 |** View of the experiment setup, grain tomography and forces between particles. (a) Schematic representation of the setup: hydrogel grains saturated with a fluorescent dye are immersed in an index-matched bath, which is subject to a uniaxial strain. A moving laser sheet illuminates cross-sections, which are imaged by a digital camera; a force sensor measures the compressive stress. (b) Superimposed image data showing part of a single slice fluorescence image (black with grey particle cross-sections), the reconstructed particles (blue) and force networks (coloured rods). The blue semitransparent spherical shapes show the reconstructed particle surfaces from one complete set of slices. Rods show forces between grains, as determined from the reconstructed grain surfaces and contact deformations. Rod thickness (thin-thick) and colour (blue-red) represent strength of the force (low-high) at the contact, scaled with the cumulative distribution function. The scale bar represents 1 cm at this position in the perspective.</figcaption></figure>
The experimental protocol consists of a series of 20 uniaxial
compressions, each followed by a corresponding expansion to the
original boundary configuration. Each cycle consists of a
compression phase imposing a strain up to 13.4% of the initial
height, followed by a decompression phase returning the top plate
to its original position. A full cycle is carried out in 60 quasi-static
steps of 1 mm each, but the top plate does not touch the grains in
the first five and last five steps. After each quasi-static step, we
carry out a complete volume scan.

**Microstructure and 3D force information.** To obtain the contact forces from the reconstructed scanned images, we implement the following steps, which we discuss further in the Methods section. First, we identify the boundary of each particle, which is not perfectly spherical, including regions of contact with other particles. Then, from the areas of contact, we use linear elasticity, and the independently measured Young modulus of the particles ($E \approx 23$ kPa), to infer the contact normal forces. The particles are nearly incompressible (volume change from uncompressed to fully compressed is <1%) so we assume that they have a Poisson ratio of $\approx 0.5$. Our best estimate of the hydrogel friction coefficient is $\mu \approx 0.03$. Hence, for these experiments, the tangential (frictional) forces are at most 3% of the corresponding normal forces and lie below experimental resolution (the minimum average contact force is $\langle f \rangle = 10^{-2}$ N, see Methods). We then compute the coarse-grained continuum stress tensor$^{19}$. We integrate the normal contribution of the stress tensor over the upper boundary to obtain the force on the top plate. We emphasize that this measurement uses only information from the micro-scale and the independently measured particle Young's modulus. The force data resulting from the stress integration and from the independent in-line force gauge are nearly identical, as we show in Fig. 2, for the last 15 of the 20 uniaxial compression/decompression cycles. Here red and blue distinguish, respectively, the macroscopic force measurement of the gauge from the microscopically derived measurement. Although there is a small offset between the two different force measurements, the overall agreement is very good, and the results are very reproducible from cycle to cycle. The difference in the two types of measurements may be due to friction, which we neglected, as well as small differences in the properties of the hydrogels between compression and decompression$^{20}$. These may also be responsible for part of the hysteresis. In the remainder of this paper, we use the force values derived from the images, for consistency with other microstructure measures.

<figure><img src="image_3.png" /><figcaption>**Figure 2 | Accuracy of the image processing.** Comparison of the compressive force *F* measured by the in-line force gauge (red) with the force inferred from the tomographic reconstruction (blue), for each of the 900 full scans taken during 15 loading/unloading cycles of the uniaxial compression experiment (the first cycle is blown-up). The top plate does not touch the grains in the first and last five scans in each cycle, but in between, the scan numbers can be viewed as an 'unfolded' measure of strain expressed in mm.</figcaption></figure>

With large data sets of particle-scale data, it is possible to
obtain reliable statistics on all quantities of interest. We measured
data from the 20 compression cycles, but we avoid the first five
cycles in the analysis as they show weak transient effects. We
compute statistics (microscopic averages noted by ⟨·⟩) only on
grains that do not touch the walls, and on contacts between such
grains, to reduce any boundary effects. We also separate the
compression and the decompression motion at the same
compression level to display the full mechanical trajectory. The
measures presented in Fig. 3 are typically averaged over 10,400
samples for contacts and 4,100 samples for grains.

These data allow us to link microscopic quantities to similar
global properties of the packing, which we discuss in the context
of Fig. 3. In particular, we provide insight into the nominally
power-law relation between the macro-scale force *F* and the total
strain *Δ*: *F* ~ *Δ*<sup>*b*</sup>. We contrast *F* and *Δ* to their microscopic
equivalents, ⟨*f*⟩ the average force at contacts, and ⟨*δ*⟩ the average
contact deformation. These last two variables are related by Hertz' law, together with the radius of curvature, *r*, at contacts. A key
observation is that *F*(*Δ*) and ⟨*f*⟩⟨(*δ*)⟩ follow substantially different
functional relations<sup>21</sup>.

**Linking microstructure and global packing response.** Figure 3a shows the uncompressed and fully compressed states corresponding to one of the cycles. In Fig. 3b, we show data for the global force, *F*, versus the global compression *Δ*. A fit of these data to a power law yields an exponent *β* = 2.2 ± 0.2, which is significantly larger than the microscopic (Hertzian) force law exponent of 1.5 measured for intergrain forces. This non-Hertzian packing mechanics has even been observed in much larger packings<sup>21</sup>, so it is not a finite size effect—the key point is that *F* depends on ⟨*f*⟩ as well as on additional microscopic properties. This is already partly revealed through the small hysteresis visible in the compression loop of Fig. 3b. The pressure, exerted on the top plate, can be expressed as a combination of the averaged microscopic quantities. Using a mean field argument detailed in the Methods section, we obtain:

$$
F \propto \langle f \rangle \langle Z \rangle \langle \phi \rangle \langle b \rangle \quad (1)
$$

where ⟨Z⟩ is the average number of contacts per particle, ⟨b⟩ is
the average distance between grain centres and ⟨ϕ⟩ is the average
packing fraction of grains within the packing (not touching the
boundaries). In Fig. 3c, we show this mean field relation is well
satisfied. This result gives a simple explanation of why the global
force response with strain has a larger effective exponent than the
particle-scale force law.

Figure 3d shows data for ⟨Z⟩ versus ⟨ϕ⟩. Near jamming, for
⟨ϕ⟩ near 0.64, we expect ⟨Z⟩ ≈ 6, whereas we measure Z between
4 and 5. This error is due to experimental uncertainty in
distinguishing between two particles that are very close, but not in
contact, versus actually in contact. By contrast, near maximal
compression, the experimental error in ⟨Z⟩ is much lower. For the
largest packing fractions, ⟨ϕ⟩ > 0.74, our particles approximately
organize into a crystalline lattice (see Methods), with a
coexistence of hexagonal compact and cubic face-centred
organization modes that implies crystal defects. Note that
⟨ϕ⟩ = 0.74 corresponds to the packing fraction for an ordered
close-packed lattice of hard spheres (FCC or HCP), although our
particles are deformable. Using an overestimate for errors in all
relevant quantities, including ⟨Z⟩, we obtain the error bars
indicated by crosses for the scaling relation of Fig. 3c and
equation 1. Note that the scaling relation of equation 1 is not
strongly affected in absolute terms by errors in ⟨Z⟩, since
the region where ⟨Z⟩ is most uncertain corresponds to small
values of ⟨f⟩.
<figure><img src="image_4.png" /></figure>

**Figure 3** | Relations between microstructure and global quantities. (a) Packing reconstruction of uncompressed and fully compressed particles, indicating the global compression. (b) Empirical power law fit of the compressive force, F, versus Δ; the average best fit gives an exponent $β = 2.2 ± 0.2$. The legend applies to all panels. (c) Scaling relation, based on a mean field pressure argument, derived in the Methods section. Error estimates are overly conservative. (d) The number of contacts per grain, $\langle Z \rangle$, varies with the volume fraction within the packing. Errors in $\langle Z \rangle$ are largest near jamming, where we underestimate this quantity. The random close packing (RCP) and hexagonal close-packed (or face-centred cubic) limits (HCP/FCC) are indicated with dashed lines.

<figure><img src="image_5.png" /></figure>

**Figure 4** | Approximate power-law relations. (a) Between $\langle \delta \rangle$ and $\Delta$, (b) between $\langle f \rangle$ and $\Delta$. The legend applies to both panels. The fits are power laws with exponents as indicated.

Now that we have established the microscopic ingredients for the macroscopic mechanical response via equation 1, we can further elucidate the micro-macro link of the constitutive relation $P(\Delta)$ or, equivalently, $F \propto \Delta^β$. Importantly, we note that $F \propto \Delta^β$ is approximately, but not strictly a power law. For example, $\langle \phi \rangle$, in equation 1, varies with $\Delta$ as $\langle \phi \rangle \propto 1/(c-\Delta)$, where the constant $c$ depends on the packing geometry. We further see that $F$ depends on $\Delta$ indirectly through $\langle Z \rangle$, so the strain-evolving topology of the microstructure matters as well—indeed, the micro-macro link is very complex. However, the fact that $F$ is roughly a power law in $\Delta$ suggests that there may be other effective power-law-like relations between microscopic and macroscopic properties. We indeed find such relations. For instance, we find that $\langle \delta \rangle$ is well approximated by $\langle \delta \rangle \propto \Delta^\gamma$, with $\alpha \approx 1.4 \pm 0.1$, providing we shift $\langle \delta \rangle$ by a small offset, $\langle \delta \rangle_{min}$, which is below the experimental resolution (Fig. 4a), and that that $\langle f \rangle \propto \Delta^\gamma$, with an empirical exponent $\gamma \approx 1.7 \pm 0.2$ (Fig. 4b). It is

important to note that we cannot directly link all the empirical exponents $\alpha$, $\beta$ and $\gamma$, as $\langle \delta \rangle$ appears both in Hertz' law and as part of $\langle b \rangle$ in equation 1 (see the Methods section).

These observations yield the notable finding that, because $\alpha > 1$, there is an effective mechanism through which the global strain $\Delta$ produces nonlinearly amplified deformations $\langle \delta \rangle$ at the grain level. This makes the nonlinear deformation amplification highly relevant, while also showing the deep entanglement of force and structure in the constitutive modelling of disordered packings of particles.

Discussion

We present a novel experimental technique that yields all microscopic contact and force vector information from a 3D packing of particles subject to controlled deformation. With these microscopic data, we have been able to verify a quantitative relation between the macroscopic mechanical response of the
packing and microscopic structural metrics of the packing. Furthermore, our experiment has revealed an important non-linear enhancement of contact deformation in response to a global packing deformation, which significantly impacts the understanding and modelling of granular materials. Our new experimental approach has great potential for yielding understanding of other granular systems, including the granular response to shear, the jamming transition, particle diffusion, effects of particle shape and so on. It supplements other 3D imaging methods, such as microCT, NMR and confocal imaging, opening up a wide range of opportunities to shed new light on the poorly understood mechanics of disordered materials.

## Methods

**Derivation of the contact forces from geometric properties.** We compute the deformation, $\delta$, at a given contact from the geometric properties of the grains involved in that contact. All grains have the same Young's modulus, so $\delta$ is split equally between the two grains$^{22}$. The radius of curvature, $r$, of the undeformed

<figure><img src="image_6.png" /></figure>

**Figure 5 | Contact and grain geometries, sketched in a 2D plane orthogonal to the contact area.** The surface of each grain (blue/green) is shown as solid lines, the hypothetical 'undeformed' surface of the grains is shown as dashed lines. The deformation, $\delta$, for this contact is the unknown quantity that we infer from linear elasticity. We measure the area of contact $A$, the position of the contact centroid $C$ and its distance $d$ to the center of mass $G$ of each grain. Note that $d$ (hence $r$) can vary for different contacts on the same particle, as particles are weakly non-spherical. The distance between grain centres $b$ depends mainly on the grain polydispersity and weakly on their non-sphericity so, on average, $\langle b \rangle \approx 2\langle d \rangle$. Our force inference algorithm uses the observed surfaces, without assuming sphericity, and provides experimental access to all the relevant quantities in order to reconstruct $\delta$.

surface (see Fig. 5), which is different for each grain and each contact, is therefore $r = d + \frac{1}{2}\delta$, with $d$ the distance between the centre of the contact area and the grain centre. We observe that the force is given by $F = E_e r_e^{1/2} \delta^{3/2}$ and $F = E_e \delta a$, with $E_e$ the effective Young's modulus, $1/r_e = 1/r_1 + 1/r_2$ the effective contact radius and $a = \sqrt{(A/\pi)}$ the radius of the area of contact, $A$. Equating the two expressions for the contact force we obtain a cubic relation, which gives us $\delta$. We have independently measured the Young's modulus for our particles: $E \approx 23$ kPa, see Fig. 6.

**Derivation of the mean field scaling relation.** The stress tensor can be computed$^{23}$ with a relation of the form $\sigma = 1/V \sum_{c \in V} b_c \otimes f_c$, with $V$ an averaging volume and $c \in V$ the contacts in that volume ($V$ is actually replaced by a smoothing kernel$^{19}$). For each contact $c$, $b_c$ is the vector between the centres of the grains and $f_c$ is the force vector along the contact normal. Neglecting friction and non-sphericity, $b_c$ and $f_c$ are nearly aligned; hence, the trace tr $(b_c \otimes f_c) \approx b_c \cdot f_c$. The number of terms in the sum depends on the density of contacts, which is about $1/2 Z \phi$, with $Z$ the number of contacts per grain and $\phi$ the grain volume fraction within the packing (for grains not touching the boundaries). The grains are nearly density matched with the surrounding fluid (density difference of $\approx 10 \text{ kg m}^{-3}$) so we can ignore the hydrostatic pressure. With roughly spherical grains, the isotropic pressure is very nearly proportional to tr $b_c \otimes f_c$. Hence, $p \propto \langle b f Z \phi \rangle$. If we write any one of the quantities, $q$, on the right as $q = \langle q \rangle + e$, the contact force, $f$, has qualitatively different statistical properties compared with $b$, $Z$ and $\phi$. The contact forces are distributed broadly over $0 \le f \le f_{\text{max}}$, while the others have relatively narrow fluctuations around a non-zero mean value. For instance, $b$ is typically twice the roughly constant particle radius minus a typical small deformation $\delta$. Since $f$ is the only quantity that can change by orders of magnitude in the experiment, we use the mean values of $b$, $Z$ and $\phi$ in the expression for $p$ and assume that deviations from the mean are uncorrelated, which means that $p \propto \langle b \rangle \langle Z \rangle \langle \phi \rangle \cdot \langle f \rangle$. The integral of $p$ over the top plate yields the top plate force, $F$, which reflects the mean field relation for $p$. When the top plate barely touches

<figure><img src="image_7.png" /></figure>

**Figure 7 | Three-dimensional realization of the distribution of directions for contacts at the maximal compression level, $\Delta = 25$ mm.** Contacts directions show a median plane with an hexagonal symmetry, typical of closed packed lattices. The directions above and below that plane correspond to either face-centred cubic or hexagonal close-packed variants.

<figure><img src="image_8.png" /></figure>

**Figure 6 | Hertz law checks.** (a) Data averaged over all contacts for the 15 compression cycles described in the main text. Here, $f$ is the force at a contact, and $\delta$ is the deformation of one of the particles at the contact. Hertz law, $f \propto r^{0.5}\delta^{1.5}$ for an individual contact, also holds for averaged quantities. That is, by a similar argument as described in the 'Methods' section, $\langle f \rangle \propto \langle r \rangle^{0.5}\langle \delta \rangle^{1.5}$ can be used in place of $\langle f \rangle \propto \langle r^{0.5}\delta^{1.5} \rangle$. This figure also gives a conservative estimate on our experimental force resolution of $\langle f \rangle = 10^{-2}$ N in the least compressed case. (b) Force measured when compressing a single hydrogel particle in a microstrain analyzer. We applied a cycle of compression/decompression up to 7.2%. During the test, the grain was immersed in a bath with temperature $21 \pm 1^\circ\text{C}$. Temperature variations do not affect the power law scaling, only the amplitude. This test verifies the Hertz law, $f \propto \delta^{3/2}$, and yields the Young's modulus $E = 23 \pm 1$ kPa for this temperature.
the grains. *F* should be null and, due to the density matching, so should *f*.
We subtract the small minimal measured values from (*f*) and *F*, which are present
at the instant the plate touches the grains and which we attribute to measurement
inaccuracies. This guarantees that *F*=0 when (*f*)=0. The coefficient of
proportionality is the averaging volume, *V*, divided by the top plate area. As shown
in Fig. 3c, this relation is remarkably well respected.

**Crystalline structure of the packing.** We analyse this crystalline structure by looking at the frequencies of the preferred directions of contacts, in each angular sector around a grain's center of mass. This gives a spherical probability distribution, the amplitude of which in each angular sector can be represented as a surface (see Fig. 7).

**Inference of 3D geometrical properties from images.** For a given state of our system, we measure images of the packing cross-sections. From these we infer the geometric properties of the grains and their contacts. We use standard image denoising, local grey level renormalization, and image rescaling to compensate varying optical path lengths. We combine all the slice images into a 3D representation in terms of voxels. We detect border voxels by using a combination of statistical tests and signal processing techniques to account for physical properties of the hydrogels (for example, homogenous grey levels within grains, compensation of artifacts induced by index mismatch at grain boundaries due to their surface properties and so on). These methods are detailed in the publicly available source code. They ensure that very few false-border outliers remain. We then uniquely attribute border voxels to grains: a local tangent plane is fit at each voxel, such that the plane normal determines the direction towards a rough estimate of the grain centre at an average radius distance. These protocentres are considered to be distributed around a unique grain identifier, which we compute with 3D kernel density estimation. This process yields a cloud of border voxels uniquely attributed to each grain. We then fit an analytic shape to these borders such that for each direction **u** (unit vector), the surface of the grain from the center of mass in that direction is given by a function $s(\mathbf{u})$ that best interpolates through all border voxels. This function $s$ is expressed in terms of a basis of spline functions on $s^2=\{\mathbf{u}\}$, the unit sphere. Unlike their Cartesian counterpart (for example, non-uniform rational B-splines (NURBS)), the triangular spherical b-splines²⁴ do not introduce any singularities on the sphere, and they are maximally isotropic. We use the convexity of the hydrogel grains as a regularizing condition, and we then perform a Levenberg-Marquardt least-square optimization to find the coefficients of $s$ in the spline basis, so as to best fit border voxels, discarding any remaining outliers. Once we have analytic descriptions of the grain surfaces, it is easy to compute their centres of mass and any other geometric quantity, including the areas of contact.

**Experimental methods.** The granular material used here consists of 514 hydrogel beads¹⁶,²⁰. The beads are approximately spherical and roughly monodisperse with a typical diameter of 2.1 ± 0.1 cm. They are immersed in a water-polyvinylpyrrolidone 360,000 MW(polyvinylpyrrolidone) solution to match the particle index of refraction to the surrounding fluid. This matching allows interior optical access, which is necessary for the refractive index-matched tomography technique used here¹⁵. Since the particles are almost entirely composed of water, they are also nearly density matched with the fluid; the particle density is <10 kg m⁻³ greater than the fluid density. By itself, the fluid-particle system is completely transparent. To obtain contrast, we dye the particles with a hydrophobic fluorescent dye (Nile Blue 690), which can be excited with a laser sheet. This laser sheet (Lasiris SNF 635 nm, 25 mW) moves on a fast linear stage, which sequentially illuminates the sample slice-by-slice. A fast camera (Basler ava1000-120) equipped with lens and long pass filter records the fluorescent image of each slice. This produces a sequence of particle cross-sections. The camera is mounted on the same stage as the laser, as in the design of ref. 25. A complete scan typically consists of 360 slices. We confine the granular system in a rectangular box with base size 16.5 × 16.5 cm and carry out cyclical uniaxial compression/decompression using a stage-controlled piston as described in the main text. The top piston is made from a 6 mm thick perforated sheet, which allows flow of index-matching fluid into and out of the packing. Except for gravitational gradients, the pore pressure inside the fluid is uniform and matched to the ambient atmosphere. The piston is driven by a linear stage (Newport MTM250) and controller (Newport XS4), with a step resolution of 1 µm. The force on the piston plate is measured with a force sensor (Loadstar RSB4-005M-A). The compression speed during compression/decompression is 0.1 mm s⁻¹ to reduce fluid-induced shear stresses on the particles and to drive the granular system completely quasi-statically. The spacing between the piston plate and the walls is small enough that particles cannot escape confinement, so the number of particles during an experiment is constant. The container height depends on the compression level, but a typical uncompressed height is about 15 cm. More experimental details can be found elsewhere²⁰.

**References**

1. Hales, S. *Vegetable Staticks* 53–54 (Oldbourne, 1727) as referenced in Aste, T. and Weaire, D. *The Pursuit of Perfect Packing* (CRC, 2007).

2. Duran, J. *Sands, Powders, and Grains* (Springer, 2000).

3. Forterre, Y. & Pouliquen, O. Flows of dense granular media. *Rev. Fluid Mech.* **40**, 1–24 (2008).

4. Jaeger, H. M., Nagel, S. R. & Behringer, R. P. Granular solids, liquids and gases. *Rev. Modern Phys.* **68**, 1259–1273 (1996).

5. van Hecke, M. Jamming of soft particles: geometry, mechanics, scaling and isostaticity. *J. Phys. Condens. Matter* **22**, 033101 (2010).

6. Dinsmore, A. D., Weeks, E. R., Prasad, V., Levitt, A. C. & Weitz, D. A. Three-dimensional confocal microscopy of colloids. *App. Optics* **40**, 4152–4159 (2001).

7. Pusey, P. N. & van Megen, W. Phase behaviour of concentrated suspensions of nearly hard colloidal spheres. *Nature* **320**, 340–342 (1986).

8. Nedderman, R. M. *Statistics and Kinematics of Granular Materials* (Cambridge Univ. Press, 1992).

9. Pitman, E. B. & Schaeffer, D. G. Stability of time-dependent compressible granular flow in two dimensions. *Comm. Pur. Appl. Math.* **40**, 421–447 (1987).

10. Majmudar, T. S. & Behringer, R. P. Contact force measurements and stress-induced anisotropy in granular materials. *Nature* **435**, 1079–1082 (2005).

11. Brujić, J. et al. 3d bulk measurements of the force distribution in a compressed emulsion system. *Faraday Discuss.* **123**, 207–220 (2003).

12. Zhou, J., Long, S., Wang, Q. & Dinsmore, A. D. Measurement of forces inside a three-dimensional pile of frictionless droplets. *Science* **312**, 1631–1633 (2006).

13. Saadatfar, M., Sheppard, A. P., Senden, T. J. & Kabla, A. J. Mapping forces in a 3D elastic assembly of grains. *J. Mech. Phys. Solids* **60**, 55–60 (2012).

14. Vlahinić, I., Andò, E., Viggiani, G. & Andrade, J. E. Towards a more accurate characterization of granular media: extracting quantitative descriptors from tomographic images. *Granular Matter* **16**, 9–21 (2014).

15. Dijksman, J. A., Rietz, F., Lörincz, K., van Hecke, M. & Losert, W. Refractive index matched scanning of dense granular materials. *Rev. Sci. Instr.* **83**, 011301 (2012).

16. Mukhopadhyay, S. & Peixinho, G. Packings of deformable spheres. *Phys. Rev. E* **84**, 011302 (2011).

17. Andrade, J. E. & Avila, C. F. Granular element method (GEM): linking interparticle forces with macroscopic loading. *Granul. Matter* **14**, 51–61 (2012).

18. Hurley, R., Marteau, E., Ravichandran, G. & Andrade, J. E. Extracting inter-particle forces in opaque granular materials: beyond photoelasticity. *J. Mech. Phys. Solids* **63**, 154–166 (2014).

19. Babic, M. Average balance equations for granular materials. *Int. J. Eng. Sci.* **35**, 523–548 (1997).

20. Dijksman, J. A., Zheng, H. & Behringer, R. P. Imaging Soft Sphere Packings in a Novel Triaxial Shear Setup. in *Powders & Grains* 1542 AIP Conference Proceedings, 457–460 (2013).

21. Lachhab, T. & Weill, C. Compression of a soft sphere packing. *Eur. Phys. J. B* **9**, 59–69 (1999).

22. Landau, L. D. & Lifshitz, E. M. *Theory of Elasticity 7 – A Course of Theoretical Physics* (Pergamon Press, 1970).

23. Zhang, J., Behringer, R. P. & Goldhirsch, I. Coarse-graining of a physical granular system. *Progress Theor. Phys. Supp.* **184**, 16–30 (2010).

24. Ying, H. & Hong, Q. *in Proceedings of Geometric Modeling and Processing*, 2004 279–287 (IEEE, 2004).

25. Lörincz, K. A. & Schall, P. Visualization of displacement fields in a sheared granular system. *Soft Matter* **6**, 3044–3049 (2010).

**Acknowledgements**

This work has been supported by NASA grant NNX10AU01G, NSF grants DMR1206351, DMS1248071 and ARO grant W911NF-1-11-0110. J.A.D. acknowledges support from the Laboratory of Physical Chemistry and Colloid Science at Wageningen University.

**Author contributions**

N.B. and J.A.D. contributed equally to this work. N.B. wrote image analysis software, J.A.D. performed the experiments. N.B., J.A.D. and R.P.B. provided analysis, interpreted the data and wrote the paper.

**Additional information**

**Competing financial interests:** The authors declare no competing financial interests.

**Reprints and permission** information is available online at http://npg.nature.com/ reprintsandpermissions/

**How to cite this article:** Brodu, N. *et al*. Spanning the scales of granular materials through microscopic force imaging. *Nat. Commun.* **6**:6361 doi: 10.1038/ncomms7361 (2015).

<figure><img src="image_9.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0

International License. The images or other third party material in this

article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/