Fluid-driven deformation of a soft granular material

Compressing a porous, fluid-filled material will drive the interstitial fluid out of the pore space, as when squeezing water out of a kitchen sponge. Inversely, injecting fluid into a porous material can deform the solid structure, as when fracturing a shale for natural gas recovery. These poromechanical interactions play an important role in geological and biological systems across a wide range of scales, from the propagation of magma through the Earth's mantle to the transport of fluid through living cells and tissues. The theory of poroelasticity has been largely successful in modeling poromechanical behavior in relatively simple systems, but this continuum theory is fundamentally limited by our understanding of the pore-scale interactions between the fluid and the solid, and these problems are notoriously difficult to study in a laboratory setting. Here, we present a high-resolution measurement of injection-driven poromechanical deformation in a system with granular microsctructure: We inject fluid into a dense, confined monolayer of soft particles and use particle tracking to reveal the dynamics of the multi-scale deformation field. We find that a continuum model based on poroelasticity theory captures certain macroscopic features of the deformation, but the particle-scale deformation field exhibits dramatic departures from smooth, continuum behavior. We observe particle-scale rearrangement and hysteresis, as well as petal-like mesoscale structures that are connected to material failure through spiral shear banding.

Poromechanical deformations are poroelastic when they are controlled by the reversible storage and release of elastic energy. The classical theory of poroelasticity couples linear elasticity with Darcy's law for fluid flow through a porous medium, and the hallmark of these systems is diffusive propagation and dissipation of fluid pressure with characteristic time scale T pe = µL 2 /(Kk), where µ is the viscosity of the fluid, L is a characteristic length scale, and K and k are the elastic modulus and permeability of the solid skeleton. This approach is valid for small deformations, but many real systems feature large deformations, small-scale microstructure, or physical mechanisms such as damage, growth, or swelling, that lead to a strongly nonlinear coupling between the pore structure and the fluid flow [5,19].
Many poroelastic deformations of practical interest are driven by fluid injection. Injection-driven deformations * christopher.macminn@eng.ox.ac.uk involve radial dilation (outward expansion), which is particularly interesting and challenging because it leads to a nontrivial state of stress and strain (e.g., [20][21][22][23]). Indeed, fluid injection into granular materials can lead to spectacular damage patterns when the injection pressure exceeds the inter-particle friction or the external confining stress (e.g., [24][25][26][27]). However, the deformation in these examples is almost completely irreversible because the solid skeleton is stiff and the fluid pressure is low, so the majority of the input energy is dissipated through frictional sliding and rearrangement, and almost none is stored elastically [28]. Fluid injection can only drive significant storage of elastic energy when the fluid pressure becomes comparable to the stiffness of the solid skeleton.
True poroelastic deformation requires either much larger pressures or much softer materials. As a result, it has proven difficult to study in a laboratory setting. Experiments with rocks and sands have been limited to postmortem inspection after high-pressure injection (e.g., [29,30]), providing useful insight into the failure of realistic geomaterials, but at a very coarse level in time and space. This limitation has been avoided in a onedimensional geometry with soft, open-cell polymer foams (kitchen sponges; e.g., [31][32][33]). However, these materials have proven to be experimentally challenging for a variety of reasons, with unsatisfying comparison between experiment and theory.
Here, we study the poromechanical deformation of a system with granular microstructure by injecting fluid into a confined monolayer of spherical particles. By using particles that are extremely soft, we construct a system that exhibits striking poroelastic phenomena at relatively low working pressures. High-resolution imaging and particle-tracking provide experimental access to the full, multi-scale deformation field. We show that the smooth, quasi-reversible macroscopic deformation can be captured in part by a minimal continuum model, despite the presence of complex shear banding and structural rearrangement.

II. FLUID INJECTION INTO A MONOLAYER OF SOFT PARTICLES
We pack a single layer of about 25,000 spherical, polyacrylamide hydrogel particles between two glass discs and we saturate the packing with the working fluid, a mixture of water and glycerol (Figure 1a; see Appendix A). The discs are separated by a permeable spacer that confines the particles, but allows fluid to leave freely around the edge. The particles are soft (having a Young modulus of ∼20 kPa), nearly incompressible (having a Poisson ratio of ∼1/2), Hertzian (exhibiting Hertz-like contact mechanics), elastic (allowing order-one elastic strains), non-cohesive, and very slippery (having low friction at particle-particle and particle-wall contacts) [34,35]. The particles have mean diameter d ≈1.2 mm with about 10% polydispersity ( Figure 5). The packing has an apparent area void fraction of ∼0.14 when viewed from above, and an actual volume void fraction of ∼0.51. The former is denser than random close packing in 2D (∼0.18) due in part to the softness of the particles, as often occurs in fluid emulsions, and also in part to their polydispersity.
To perform an experiment, we inject more of the same working fluid into the cell at a constant volume rate Q. This fluid enters the cell via an injection port in the center of the lower disc, flows radially outward through the packing, and exits through the spacer at the outer edge. The resulting fluid pressure gradient within the porous layer (large pressure at the center dropping to atmospheric at the edge) deforms the packing, driving the particles radially outward and opening a cavity in the center (Figure 1b,c; Video S1). This coupling of fluid pressure and solid deformation is the core idea behind poromechanics. The deformation eventually reaches a steady state (here, after ∼100 s) in which the gradient of elastic stress in the solid skeleton balances the gradient of pressure in the fluid. We then stop injecting, at which point the elastic stress relaxes as the pressure gradient dissipates and the cavity closes. The relaxation of the packing highlights the macroscopically elastic nature of the system, demonstrating that the packing stores elastic energy during the injection phase and releases it during the relaxation phase. We repeat this injectionrelaxation cycle several times in the same packing. We image the deformation and subsequent relaxation of the packing with a digital still camera, detecting the particle positions in each image to within about 0.01d and then tracking the particles from image to image (Figure 1b; see Appendix A).

III. MULTI-SCALE DEFORMATION FIELD
One striking aspect of the deformation is the cavity that opens and then closes in the center of the packing. Despite the irregular shape of the cavity (Figure 2b), we find that the macroscopic dynamics of its expansion and collapse are smooth and relatively reproducible across repeated injected-relaxation cycles (Figures 2a and 7). In contrast, the shape of the cavity varies from cycle to cycle (Figure 2b; Video S2). The size of the cavity increases with injection rate roughly in accordance with the prediction of a minimal continuum model (described in more detail below), but repeating an experiment at a given injection rate after "resetting" the packing by completely rearranging the particles leads to large variability (Figure 2c). This implies that the macroscopic properties of the packing are a strong function of particle arrangement.
We measure the internal deformation of the packing via particle tracking, which provides a direct measure of the displacement field. For this purpose, we define a rectangular coordinate system centered at the injection port, where (x i , y i ) is the position of particle i at time t and (X i , Y i ) is its initial position. The displacement of particle i is then The deformation is primarily radial because of the axisymmetric boundary conditions, so we focus on u r . The difference u − |u r | is a measure of the non-radial component of the displacement, which we find to be a few percent or less of the radial component (see Appendix A).
We find that the radial displacement is large near the cavity and fades to zero at the rigid edge, with a petallike mesoscale structure (Figure 3a; Video S3). Similar petal-like features have been observed in simulations of fluid injection into an initially dry packing of frictional particles [36], indicating that these structures are not an artifact of our low-friction system. Additionally, similar but much more regular features have been observed in quasi-static "grain injection" experiments [37], where they were identified with preferential directions in the far-field crystal structure. However, our packings are isotropic due to the polydispersity of the particles.
Each petal represents a group of particles that move radially outward further than their neighbors, implying that the edges of each petal are bands of localized shear failure. We confirm this by calculating the local strain field ( [38]; see Appendix A), revealing a network of spiral shear bands that span the entire system (Figure 3b). Shear bands following logarithmic spirals are a well-known feature of failure in radial dilation and hydraulic fracture [21,29]. We see spirals with a pitch of roughly of π/4, which implies that the packing has a very low shear strength (i.e., a very low friction angle). Correlations between shear strain and positive volumetric strain are evidence of shear dilation, a well-known feature of deformation in granular materials (Figure 3b FIG. 2. Fluid injection into the center of the packing drives outward compaction, opening a cavity. When injection stops, the packing relaxes and the cavity closes. Here, (a-b) we inject at a steady rate of Q = 16 mL/min for about 135 s and then stop, allowing the packing to relax. We repeat this injection-relaxation cycle two more times (black, then blue, then red).
(a) Comparing the area-vs-time curves from the three cycles (offset by 12 s for comparison; cf. Figure 7) shows that the process is macroscopically smooth, with similar dynamics and maximum area in each cycle. The first cycle (black) is somewhat different than subsequent cycles due to heterogeneities in the initial particle distribution. We also show the prediction of a continuum model using best-fit material properties (dashed gray curve); we discuss the model in §V. (b) The steady-state cavity shape is neither smooth nor repeatable, indicating the presence of irreversible micromechanics. The cavity does not open symmetrically about the injection port (gray circle with diameter ∼2.5 mm). The scalloped edges of the cavity profiles are particle-scale roughness (∼1 mm). (c) We repeat this experiment at different injection rates, showing here the maximum (i.e., steady-state) cavity area against injection rate Q. Each group (color) of circles indicates a series of least two cycles after "resetting" the packing by removing the particles, cleaning the apparatus, and replacing the same particles (the cross is a single cycle). The variability between cycles is relatively small except in one case (red), indicating that cycle-to-cycle irreversibilities have a weak impact on the macroscopic mechanics. The variability between series at the same injection rate is much larger, indicating that the particle arrangement has a strong impact on the macroscopic mechanics. All packings have initial porosity between 0.506 and 0.515, and we do not observe any clear correlation between initial porosity and cavity size in this range. The black curve is the prediction of the continuum model using the same mechanical properties for all points. irrev.
FIG. 3. The displacement field is characterized by a petal-like mesoscale structure that corresponds to spiral bands in the shear strain field. Here we show (a) the radial displacement, which reveals detailed mesoscale structure reminiscent of flower petals. We show (b) the magnitude of the shear strain with several logarithmic spirals for comparison (dashed black). The spirals have a pitch of π/4, which is consistent with the shear failure of a material with negligible shear strength. We also show (c) the volumetric strain. The volumetric strain shows expansion in the inner region and compression in the outer region, and is relatively axisymmetric. Lastly, we include (d) the typical non-affine displacementDmin (see Appendix A). The nonaffine displacement is much smaller than the radial displacement except near the cavity, where the kinematics become strongly non-affine. The strain field is not necessarily a good measure of the deformation in areas with large amounts of non-affinity (near the cavity). We plot these quantities from an Eulerian perspective (against current radial position r). We decompose the deformation field into (e-h) reversible and (i-l) irreversible components. For comparison, we project all fields onto the deformed configuration at steady-state. The reversible component is larger in magnitude than the irreversible component, and contains most of the mesoscale structure. The irreversible component contains much of the non-affine displacement.

IV. ELASTICITY, PLASTICITY, (IR)REVERSIBILITY, AND DISSIPATION
Macroscopically, elastic deformations involve the reversible storage and release of strain energy, whereas plastic deformations are dissipative and irreversible. In crystalline solids, there is a clear distinction between elasticity and plasticity at the particle scale: Elastic deformations involve stretching or compressing bonds between atoms or molecules, whereas plastic deformations involve breaking and/or rearranging bonds. This distinction is less clear in amorphous or granular materials, where deformations often involve a combination of reversible and irreversible rearrangement events [39][40][41].
Here, the fact that the cavity closes completely upon relaxation implies that the deformation is macroscopically reversible. However, the hysteresis in cavity shape is evidence of particle-scale irreversibility, and the shear bands are evidence of plastic failure. To investigate the apparent contradiction of strongly irreversible micromechanics coexisting with smooth, quasi-reversible macroscopic mechanics, we decompose the deformation field into reversible and irreversible components. We calculate these by considering the transformation between three configurations: the initial (before the deformation), the deformed (at steady-state), and the final (relaxed, after the deformation). The total strain E is that which transforms the initial configuration into the deformed one. The irreversible strain E irr is that which remains after the deformation relaxes (i.e., the residual strain); this transforms the initial configuration to the final relaxed one. The reversible strain E rev is that which dissipates as the deformation relaxes; this would transform the final relaxed configuration back to the deformed one. For infinitesimal deformations, these three strains are related by superposition, E = E rev + E irr ; we calculate them independently since the deformation is large. Since we calculate strain as the locally affine best fit to the actual deformation field, we also calculate the root-mean-square difference between the affine field and the actual deformation field,D min ( [38]; see Appendix A). This is a measure of the typical non-affine displacement, which is indicative of the amount of particle-level rearrangement.
Comparing the reversible and irreversible components of the strain field (Figure 3f,g vs. j,k), we find that the inner region is dominated by a combination of reversible and irreversible volumetric expansion (positive volumetric strain) and shear. Volumetric expansion indicates that particles have traveled away from their neighbors, which is expected near the cavity since the particles move radially outward by several diameters. Since the packing cannot support tension, this leads to local collapse or "unjamming" of the packing structure, which leads to large amounts of both reversible and irreversible rearrangement (Figure 3d,h,l). In contrast, the outer region is dominated by relatively smooth, axisymmetric, reversible volumetric compression (negative volumetric strain). The displacement and shear strain are much smaller, and there is much less rearrangement. The spiral shear bands span the entire system and, surprisingly, are primarily reversible (Figure 3b,f,g).
All rearrangements play a strong role in the dynamics since reconfiguration of the packing takes time and dissipates energy. Macroscopically, dissipative deformations that are reversible are known as "viscoelastic", whereas those that are irreversible are known as "viscoplastic". Unlike elasticity, which is quasi-static, these viscous processes are rate dependent.

V. POROVISCOELASTIC CONTINUUM MODEL
The steady-state deformation is set by the balance between the gradient in fluid pressure within the packing and the roughly axisymmetric elastic compression of the outer region of the packing. Motivated by this, we next derive a minimal axisymmetric model for this system based on the theory of poroelasticity [5]. Our model is intended to capture four main features of this system: (1) conservation of volume, (2) poromechanical coupling between pressure gradients in the fluid and stress gradients in the solid skeleton, (3) elastic energy storage in the solid skeleton, and (4) viscous dissipation due to reversible and irreversible rearrangements. We do not attempt to capture the evolution of effective material properties due to irreversible rearrangements (c.f., Figure 2c). We emphasize that we are not attempting to develop a general model for deforming granular materials, but rather a minimal one that captures the leading order behavior of our poroelastic system. We outline the core assumptions of the model here and present a detailed derivation in Appendix B.
We assume that the packing is homogeneous, and that the flow and deformation fields are axisymmetric. We also assume that the fluid and the solid are individually incompressible which, for the solid, implies that the beads can rearrange and deform without changing volume. This is justified here because the working pressure is low relative to the bulk moduli of the fluid and the particles (∼5 kPa vs. ∼2 GPa). This allows for a simple but exact kinematic relationship between the volumetric strain and the local porosity (fluid or void fraction) φ f (r, t) (Equation B10). We assume that the porosity is initially uniform and equal to φ f (r, 0) = φ f,0 .
We assume that the elastic stress in the solid skeleton is isotropic, meaning that the skeleton stores elastic energy due to volumetric compression but not due to shear. This is justified by the fact that similar packings are known to have an anomalously low ratio of shear to bulk modulus due to the extreme softness and slipperiness of the contacts [42]. This allows us to link the stress directly to the volumetric strain, and therefore to φ f . This is useful because although the displacements are large, the changes in porosity are small. We therefore take the elastic component of the stress to be a power law in the change in porosity with exponent 3/2, as appropriate for a granular material consisting of Hertzian particles [43], with an effective (drained) bulk modulus K.
The shear strain field indicates that the skeleton experiences shear failure near the cavity and along the spiral shear bands. The shear failure of granular materials is typically modelled with a frictional (Mohr-Coulomb) failure criterion, which states that the material will yield (fail plastically) when the shear stress anywhere exceeds some fixed fraction of the local normal stress. After yielding, the structure of the material rearranges according to a suitable (visco)plastic flow law. This transition to plastic flow (unjamming) has been studied extensively from a variety of perspectives (e.g., [38,[43][44][45]).
Although the dynamics of unjamming can be extremely important in many systems, we do not treat this behavior here because our system is compactioncontrolled. That is, the rate of shear strain near the cavity is fundamentally limited by the rate of volumetric strain in the outer region since the cavity can only expand as fast the outer region compresses. This is consistent with the fact that we expect the effective viscosity of rearrangement in the compacting outer region to be much larger than that in the unjamming inner region [45]. In an expansion-controlled system where the outer region is not yet jammed, we would expect the cavity shape to exhibit regular, sharp, triangular cracks [46]. In contrast, we see irregular and relatively smooth cavity shapes.
We model viscous dissipation due to rearrangement in the outer region in a very simplistic way by assuming that this contributes a transient component to the volumetric stress that is linear in the rate of change of porosity with an effective viscosity η. This linear, Kelvin-Voigt-like representation is often used in viscoelasticity, but here it is intended to capture both reversible and irreversible rearrangements. The linear viscous term embodies the rate-dependence of these processes, introducing a characteristic time scale for viscous rearrangement, T vr = η/K. Although the viscosity should itself be a function of the local volumetric stress [45], we ignore this for simplicity.
Finally, we assume that fluid flows relative to the solid skeleton according to Darcy's law with a constant permeability k. The assumption of constant permeability is justified by the fact that the changes in porosity are relatively small, at most a few percent.
These assumptions lead to a nonlinear conservation law for the local porosity as a function of time: is the normalized change in porosity, which is a measure of the Eulerian volumetric strain;t andr are time and radial position scaled by the poroelastic time scale T pe = µb 2 /(Kk) and the outer radius b, respectively; and is the stress in the solid skeleton (the effective stress) scaled by the bulk modulus K. The model has two dimensionless parameters, α = µQ/(2πhKk) and β = ηk/(µb 2 ). The former compares the pressure gradient in the fluid with the stiffness of the solid skeleton, while the latter compares the viscous time scale to the poroelastic one. Because the mechanical properties of the packing are difficult to measure and strongly dependent on the particle arrangement, we use α and β as fitting parameters. The former influences the rate of deformation and sets the steady-state; the latter influences only the rate. We solve Eq. (1) subject to two boundary conditions, which are that the cavity wall is a free surface where the stress in the solid skeleton vanishes, and that the outer edge is a rigid boundary where the solid is stationary. We also solve simultaneously an evolution equation for the position of the cavity wall. This model is a radial version of those that have been developed for the rectilinear deformation of kitchen sponges [33]. A similar radial model was developed in the context of blood-vessel pressurization [47], but our model is kinematically exact for finite axisymmetric deformations (see Appendix B) and we incorporate a nonlinear elastic stiffness as well as viscous dissipation.
In steady state, the model can be simplified to an integral coupled with a nonlinear algebraic equation for the cavity radius; this is straightforward to evaluate numerically. Despite the simplicity of the model and the large grain-scale variability in the experiments, we find that, after choosing α to match the cavity area, the steadystate displacement field predicted by the model agrees surprisingly well with the azimuthal average from the experiments (Figure 4a). We also calculate an azimuthally averaged porosity field from the displacement field via Equation (B10). The model does not agree very well with this (Figure 4b), although it does capture certain aspects: The initial porosity in the experiment is φ f,0 ≈ 0.51, and both the model and the experiment show a concave-up trend from a value near 0.51 near the cavity down to a value near 0.50 at the wall. The measured porosity field is non-monotonic with a minimum value at r/b ≈ 0.5, whereas the model is monotonic with its minimum at r/b = 1, but the minimum values are similar.
The model is able to capture the displacement field relatively well because the displacement field is an integral measure of the deformation, dominated by the mean features and relatively insensitive to the local details. In contrast, the porosity field is a direct measure of the local strain, strongly reflecting the local details of the deformation. In addition, use of Equation (B10) assumes smoothness, axisymmetry, and an initially uniform porosity field, but the experiment does not necessarily observe any these except in an averaged sense. A more elaborate model might take advantage of plastic-failure theory from soil mechanics [20], the theory of shear-transformation zones in amorphous solids [38,48], or constitutive laws from the rheology of suspensions [45], although any of these approaches would lead Despite the large and structured variability in the particle displacements, the continuum model agrees well with the azimuthally averaged displacement profile after choosing α to match the cavity radius. Here we show (a) the steady-state radial displacement of the particles (gray dots), the radius of the cavity (red circle), and two different measures of the azimuthally averaged radial displacement (red curve: direct azimuthal average; blue curve: integrated volumetric strain). We also show the azimuthally averaged non-affine displacement (solid gray curve) and the displacement field predicted by the continuum model (solid black curve). The displacement becomes strongly non-affine near the cavity, so we choose a threshold to the left of which affine quantities are poor measures of the deformation: The dashed portions of the red and blue curves are where the nonaffine displacement accounts for more than 10% of the total. We calculate (b) two corresponding measures of the azimuthally averaged porosity by differentiating the averaged displacement according to Eq. (B10) and smoothing the result (red and blue curves). This is a poor measure of porosity near the cavity because it assumes a smooth, axisymmetric displacement field. The model (solid black curve) does not capture the porosity field very well. We plot these quantities from an Eulerian perspective.
to a substantial increase in complexity and several additional constitutive parameters. Solving the time-dependent model numerically, we find that it captures the dynamics of cavity expansion after choosing β accordingly (Figure 2a). The fit (β ≈ 5) yields a poroelastic time scale of ∼1.5 s and a viscous one of ∼5 s, implying that viscous rearrangement controls the overall rate of cavity expansion. Given the complex nature of this process, it is surprising that our linear Kelvin-Voigt model captures the dynamics of cavity expansion as well as it does. The model does not capture relaxation very well, relaxing much more slowly than the experiment. This is not surprising since the amount of viscous dissipation due to rearrangement is likely different during relaxation than during cavity expansion. Hysteresis in effective material properties is common in plasticity, soil mechanics, and granular materials, but our minimal model does not include this.

VI. CONCLUSIONS
Fluid injection into a soft granular material drives deformation that is macroscopically poroelastic, despite rich micromechanical complexity. We found that the deformation in the inner region (near the injection port) is dominated by irreversible structural plasticity that leads to strong variations in the cavity shape, whereas the outer region deforms smoothly and reversibly. The latter ultimately supports the radially outward loading and controls the macroscopic mechanics of the steady state, such that the leading order features of the deformation can be captured relatively well with an axisymmetric continuum model. We expect this coexistence of microscopic irreversibility with macroscopic reversibility to be a strong feature of elastic dilational deformation in any system with a low ratio of shear strength to bulk stiffness. Many intriguing problems remain, such as connecting changes in grain-scale structure with the evolution of macroscopic properties, examining more sophisticated material models based on plasticity theory, and exploring the dynamics of fluid-driven shear failure.
More broadly, this system is a promising platform for high-resolution measurement of the dynamics of poromechanical deformation, and it has several additional features that we did not take advantage of here. For example, polyacrylamide hydrogel is closely index-matched with water, making it extremely well-suited to visualization in three-dimensional systems [34,35,49,50]. Polyacrylamide hydrogel is also sensitive to both temperature and dissolved salt concentration. This allows for precise system-wide tuning of the size and stiffness of the particles; it also enables studies of the dynamics of swelling or shrinking in response to local or global stimuli, which has particular relevance to biophysical systems [7,9,11].
In much the same way that granular monolayers and rafts of bubbles have served as indispensable model systems for developing fundamental concepts in the mechanics of both crystalline and amorphous materials [51][52][53][54], so too can packings of soft particles provide unique in-sight into the deformation and failure of materials under nontrivial poromechanical loading, from the propagation of poroelastic waves to the coupling of deformation with flow and transport. In addition to serving as a tool for benchmarking numerical simulations [36,55,56], this system offers an avenue into the experimental exploration of other fundamental problems of poroelasticity that have previously existed only as theoretical predictions or inferences from field observations. The discs are 19 mm thick and 212.7 mm in diameter, and they are separated by a plastic spacer that defines a working area of diameter b=210.5 mm and thickness h=1.44 mm. This is two SDs greater than d, confining the particles to a two-dimensional monolayer without restricting their in-plane motion. The spacer is permeable, allowing fluid to exit while confining the particles. We inject fluid at a fixed volume flow rate using a syringe pump (New Era NE-4000). The fluid is a mixture of water and glycerol (61% glycerol by mass) with a viscosity of µ=0.012 Pa·s at 20 • C. We acquire images with a digital camera (Canon EOS Rebel T2i).
b. Particle detection and tracking. We process the experimental images in MATLAB, detecting particle positions in each image via centroid-finding after applying an intensity threshold. We track the particles from image to image using a standard particle-tracking algorithm [57]. The images have a resolution of about 7.9 pixels per 1 mm (about 9.4 pixels per particle diameter). Frame-to-frame detection noise for particle centers is about 0.1 pixels (13 µm or about 1% of one diameter; Figure 6).
c. Deformation field. We use the particle positions to calculate a best-fit local strain field following [38]. However, the quantity described in [38] as a local strain tensor is more correctly identified as a local displacement gradient tensor ∂u/∂X, where X is the undeformed configuration and u is the displacement field. The distinction is important here, where the displace- ments are large. We calculate the displacement gradient tensor and then use it to calculate the Green-Lagrange strain tensor E = 1 2 (F F − I), where F = I + ∂u/∂X is the deformation gradient tensor and I is the identity tensor. The strains we report above are Green-Lagrange strains. We also calculate the root-mean-square (RMS) non-affine displacementD min from their quantity D 2 min by dividing each local value by the number of neighbors to form a mean (their quantity is a sum) and then taking the square root. The result has dimensions of length. We estimate this detection noise by calculating (a) total particle displacement between two frames taken before injection has started, when the particles are at rest. The noise has a mean of about 0.006d particle diameters with a standard deviation of about 0.005d, indicating a detection threshold of ∼0.01d (vertical blue line). At steady-state deformation, (b) the average radial displacement is about ∼1.4d with a standard deviation of about 1.6d, whereas (c) the average nonradial displacement is about 0.03d with a standard deviation of 0.06d, indicating that particle motion is almost entirely radial.
where v f (r, t) and v s (r, t) are the velocity of the fluid and of the solid, Q is the volume rate of fluid injection, and h is the thickness of the gap. We assume that fluid flows relative to the solid skeleton according to Darcy's law, where k and µ are the permeability of the solid skeleton and the viscosity of the fluid, respectively, and p(r, t) is the fluid pressure. Poroelastic theory dictates that the internal gradient in fluid pressure acts as a body force on the solid skeleton, and mechanical equilibrium requires that this must be supported by the divergence of the stress in the solid skeleton. We expect that the packing has an extremely low ratio of shear to bulk modulus, so for simplicity we assume that the solid cannot support shear or tensile stresses. This implies that the stress tensor is isotropic, and we write mechanical equilibrium as where σ is the effective stress (i.e., the stress in the solid skeleton). Combining all of the above, we obtain a conservation law for the evolution of the porosity, is the Eulerian volumetric strain, or the normalized change in porosity, and φ f,0 is the initial (relaxed) porosity.
We must now specify a constitutive relationship between effective stress and strain, but note that Equation (B4) is valid for any stress-strain relationship that yields an isotropic stress tensor. Here, we take the stress to be Hertzian elastic and linearly viscous in the volumetric strain. This can be written as where K and η are the effective bulk modulus and viscosity of the solid skeleton, respectively. Note that the stress vanishes at φ f = φ f,0 →φ f = 0. The linear viscous term introduces a simple "rearrangement" timescale, without which the mechanics would be quasi-static. We assume that fluid is injected into a cavity in the solid of radius a(t) and initial radius a(0) = a 0 . The behavior is independent of a 0 for a 0 b. The fluid and the solid are initially at rest, v f (r, 0) = v s (r, 0) = 0, and the initial porosity field is φ f (r, 0) = φ f,0 →φ f (r, 0) = 0. Injection begins at t = 0.
At the inner boundary, r = a, the normal component of the effective stress must vanish since the cavity is a free surface of the solid skeleton. For an isotropic stress field, this implies that σ (a, t) = 0 and therefore that φ f (a, t) = 0 (B6) forφ f (r, 0) = 0. At the rigid outer boundary, r = b, we have that u r (b, t) = v s (b, t) = 0. This can be written We also require an evolution equation for the position of the moving inner boundary, Equations (B4)-(B8) constitute a one-dimensional, timedependent, moving-boundary problem. The solid displacement field does not appear explicitly, but we can calculate it at any time through a kinematic relationship, where J is the Jacobian determinant, or the Jacobian of the deformation, and F is the deformation gradient tensor (see Appendix A). For an axisymmetric, planestrain deformation, Eq. (B9) becomes where u r (r, t) is the r-component of the (Eulerian) solid displacement field. In contrast with linear poroelasticity, this problem is nonlinear for four reasons: We account rigorously for (1) the moving boundary, (2) the solid velocity, and (3) the exact relationship between porosity and displacement, and we use (4) a nonlinear elastic law. A similar axisymmetric model was developed by [47] for small elastic deformations. Our model incorporates a nonlinear elastic law and linear viscous effects, and is kinematically exact-it is valid for arbitrarily large deformations as long as the constitutive laws remain valid.

Dimensionless form
We present Equations (B4) and (B5) in dimensionless form in the main text (Equations 1 and 2, respectively), wherer = r/b is the radial coordinate scaled by the outer radius,t = t/T pe is time scaled by the poroelastic time scale,σ = σ /K is the effective stress scaled by the bulk modulus, and the two dimensionless parameters are α = µQ/(2πhKk) and β = ηk/(µb 2 ). The dimensionless boundary conditions arẽ r [α ln(r/ã ss )] 2/3 dr (B16) and the subscript "ss" refers to the value of a quantity at steady state. The steady-state cavity radius is determined implicitly by the definitionũ r,ss (ã ss ) =ã ss −ã 0 , which leads toã ss = ã 2 0 + 2I(ã ss ).
We solve Equation (B17) numerically using a standard root-finding technique, evaluating the integral I(ã ss ) from Equation (B16) using the trapezoidal rule. The mechanical properties of the packing are difficult to measure and strongly dependent on the particle arrangement (Figure 2c), so we use the mechanical properties as fitting parameters. To compare the model with the experiment at steady-state, as in Figure 4, we choose a value of the product Kk to match the steady-state cavity radius. We use a single value of Kk for all of the experiments in Figure 2c. To compare the dynamics, as in Figures 2a and 7, we choose Kk to match the steadystate cavity radius and then ηk to match the rate of deformation and relaxation (ηk plays no role in the steady state).  Figure 2a of the main text, but without the 12 s offset. The second and third cycles (blue and red) and the model (dashed gray) agree very well during cavity expansion, but the model relaxes more slowly than the experiment. The first cycle of the experiment (black) is somewhat different due to heterogeneities in the initial condition.