Jamming transition and emergence of fracturing in wet granular media

We study fluid-induced deformation of granular media, and the fundamental role of capillarity and wettability on the emergence of fracture patterns. We develop a hydromechanical computational model, coupling a “moving capacitor” dynamic network model of two-phase flow at the pore scale with a discrete element model of grain mechanics. We simulate the slow injection of a less viscous fluid into a frictional granular pack initially saturated with a more viscous, immiscible fluid. We study the impact of wettability and initial packing density, and find four different regimes of the fluid invasion: cavity expansion and fracturing, frictional fingers, capillary invasion, and capillary compaction. We explain fracture initiation as emerging from a jamming transition, and synthesize the system’s behavior in the form of a phase diagram of jamming for wet granular media.

Immiscible fluid-fluid displacement in porous media is important in many natural and industrial processes, including the displacement of air by water during rainfall infiltration [1], storage of carbon dioxide in deep saline aquifers [2], contaminant soil remediation [3], enhanced oil recovery [4], and design of microfluidic devices [5]. While fluid-fluid displacement in rigid porous media has been studied in depth, fundamental gaps remain in our understanding of the interplay between multiphase flow in a granular medium and the displacement of the grain particles [6,7]. This interplay can lead to a wide range of patterns, including fractures [8][9][10][11][12][13][14], desiccation cracks [15,16], labyrinth structures [17], and granular and frictional fingers [18][19][20][21]. There are several controlling parameters behind the morphodynamics that govern the transition between the different regimes. A modified capillary number Ca * characterizes the crossover from capillary fingering to viscous fingering [22], and a transition from fingering to fracturing can be achieved either by decreasing frictional resistance [22], or setting the outer boundary as free [23]. The balance between frictional, viscous, and capillary forces has been studied in experiments [17,21,22] and simulations [10,24], and has helped understand the underlying mechanisms for a wide range of phenomena, including venting dynamics of an immersed granular layer [25][26][27], fractures in drying colloidal suspensions [8,12], and methane migration in lake sediments [28][29][30][31].
As one of the factors that influences multiphase flow in porous media, wettability (the relative affinity of the substrate to each of the fluids, and measured by the contact angle θ ) has * juanes@mit.edu been studied for decades. While much is now known about the role of wettability on multiphase displacements in porous media [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47], fundamental gaps remain in the context of grain-scale mechanisms and their macroscale consequences. Given the importance of capillarity on the fracture of granular packs [10,14,21,22,24], here we focus on the impact of wetting properties on the emergence of such fracture patterns. We also adopt packing density as a control parameter, since it can lead to a transition from Saffman-Taylor instability to dendritic (or ramified) fingering patterns [48], or from frictional fingering to stick-slip bubbles [21].
In this Rapid Communication, we uncover four fluidinvasion morphological regimes under different initial packing densities and substrate wettabilities: cavity expansion and fracturing, frictional fingers, capillary invasion, and capillary compaction. To rationalize these simulation outputs, we propose to analyze the evolution of the system as one approaching a jamming transition, which provides insights that allow us to map the wealth of behavior map the wealth of behavior onto a phase diagram of jamming for wet granular media.
We adopt a recently developed "moving capacitor" dynamic network model to simulate fluid-fluid displacement at the pore level [44] (see Supplemental Material [49]). The model employs an analog of the pore network geometry, where resistors, batteries, and capacitors are responsible for viscous, out-of-plane, and in-plane Laplace pressure drops, respectively. The fluid-fluid interface is represented as a moving capacitor-when the interface advances, the Laplace pressure increases until it encounters a burst (equivalent to a Haines jump), touch (touches the nearest particle), or overlap event (coalesces with a neighboring interface) [35,36,43]. These events determine how the interface advances, enlisting one or more new particles when a node on the interface reaches its filling capacity and becomes unstable. This model reproduces both the displacement pattern and the injection pressure signal under a wide range of capillary numbers and substrate wettabilities [43,44,50]. To capture particle motion, we couple the dynamic flow network model with a discrete element model (DEM), PFC2D® [49,51]. Hydromechanical two-way coupling is achieved from three perspectives: (1) The fluid pressures calculated from the moving-capacitor flow model exert forces on particles, and lead to particle rearrangement and deformation; (2) particle movements change the geometric configuration of the granular pack, which in turn changes the pore network topology and throat conductances and capillary entry pressures; and (3) expansion of the central cavity around the injection port "consumes" injected fluid, which decreases the flow of fluid permeating through the granular pack.
We simulate immiscible fluid-fluid displacement through a granular pack confined in a circular flow cell, by setting a constant injection rate at the center, and constant pressure at the perimeter. The invading and defending fluid viscosities are set to η inv = 8.9 × 10 −4 Pa s for water, and η def = 0.34 Pa s for oil, respectively, and the interfacial tension is set to γ = 13 × 10 −3 N/m. These parameters are chosen to mimic the experiments of Zhao et al. [41]. The granular pack has an outer and inner radius of R out = 13.25 mm, R in = 0.5 mm, and a height h = 330 μm. We adopt a simplified Hertz-Mindlin contact model [51] for particles in the granular pack, with the following properties: shear modulus G = 50 MPa, Poisson ratio ν = 0.5 (quasi-incompressible, as in Ref. [52]), coefficient of friction μ = 0.3 [22], density ρ = 1040 kg/m 3 , and mean diameter d = 300 μm with 10% standard deviation (the same polydispersity as in Ref. [52]). We choose an injection rate Q inj = 4.3 × 10 −11 m 3 /s, corresponding to a modified capillary number Ca * = η def Q inj R out /(γ hd 2 ) = 0.5 [22], for which viscous pressure gradients have time to relax between front movements, and capillary effects govern the displacement [53]. We conduct simulations in which we fix these parameters, and we vary the contact angle θ from 140 • (drainage) to 46 • (imbibition), and the initial packing density φ 0 from 0.68 (loose pack) to 0.84 (dense pack).
In Fig. 1, we show the fluid invasion morphologies that result from injection in the form of a visual phase diagram for different values of θ and φ 0 . The collection of patterns at breakthrough-when the invading fluid first reaches the outer boundary-exhibits four different regimes: (I) cavity expansion and fracturing, (II) frictional fingers, (III) capillary invasion, and (IV) capillary compaction.
To elucidate the conditions that lead to the emergence of each type of invasion pattern, we analyze the time evolution of the interface morphology and injection pressure for representative cases of each regime (see Fig. 3 of Supplemental Material and supplemental videos [49]).
Regime I: Cavity expansion and fracturing. When the injection pressure from fluid injection is sufficient to push particles outwards, the cavity keeps expanding until the energy input becomes insufficient to compact the granular pack further; the point at which fractures emerge [Supplemental Fig. 3(a) in Ref. [49]]. The wide range in P cap at breakthrough (t d → 1) confirms the vulnerability of fracture tips compared with other throats along the cavity perimeter.
Regime II: Frictional fingers. At only weakly wetting conditions, the injection pressure is positive but smaller than in drainage. In this case, the injected fluid pushes away particles in certain directions, preferably those with loosely packed particles, and develops frictional fingers [Supplemental Fig. 3(b) in Ref. [49]].
Regime III: Capillary invasion. When particles have been densely packed initially, a small injection pressure (either positive or negative) is insufficient to overcome the established chains of contact forces, and thus particles do not move. In this case, we observe patterns of capillary fluid invasion in rigid media [Supplemental Fig. 3(c) in Ref. [49]]. The crossover from capillary invasion to capillary fracturing can be triggered, as we demonstrate here, by increasing θ to increase capillary forces.
Regime IV: Capillary compaction. In strong imbibition the injection pressure is negative, and for sufficiently loose granular packs, particles are dragged into the invading fluid under the out-of-plane curvature effect, leading to capillary compaction [Supplemental Fig. 3(d) in Ref. [49]].
The temporal signal of the injection pressure encodes information needed to understand the interplay between particle movement and fluid-fluid displacement. Since we restrict our study to the case when capillary forces dominate and viscous dissipation is negligible, the injection pressure signal is determined by the capillary entry pressure P cap , which is a sum of in-plane and out-of-plane components. As a result, the injection pressure shows fluctuations in a stick-slip manner for all θ and φ 0 , as has been documented in slow drainage experiments [53][54][55] and simulations [44]. As θ decreases, indicating that the substrate becomes more wetting to the invading fluid, the fluid-fluid displacement is controlled by cooperative pore-filling events (touch and overlap) with smaller P cap compared with burst events [35,36,43,44]. This explains the general decreasing trend of injection pressure as θ decreases [ Fig. 2(a)].
In a drainage displacement, instead of fluctuating around a mean value [44], the injection pressure exhibits a surprising convex shape as a function of time, first decreasing and then increasing with time. This is a signature of the fluid-solid coupling: The particles around the cavity are separated (opening up the throats and decreasing P cap ) during the initial stages of expansion, and then brought closer together (narrowing the throats and increasing P cap ), as the granular pack is being compacted during the late stages [ Fig. 2(a)]. Figure 1 exhibits a surprising and heretofore unrecognized behavior of fluid injection into a granular pack: A decrease in θ -that is, transitioning from drainage to weak imbibitionleads to an earlier onset of fracturing, as evidenced by the smaller size of the fluid cavity at fluid breakthrough. This behavior cannot be explained by the evolving injection pressure level, or the evolving packing fraction outside the cavity, or This raises the question of how wettability impacts the onset of fracturing, and whether such dependence is amenable to prediction. To answer this question, we hypothesize that the emergence of fracturing is akin to a phase transition from liquidlike to solidlike behavior, and that therefore it can be understood as a jamming transition. Indeed, the jamming transition has proved instrumental in understanding mechanical integrity in a remarkably diverse range of systems [56]. Examples include colloidal suspensions [57], athermal systems such as foam and emulsions [58], and the glass transition in supercooled liquids [59,60]. The jamming transition also occurs in (dry) granular systems at a well-defined packing density φ c when the conditions of mechanical stability are satisfied [61][62][63][64][65]. Here, we explore whether the concept of jamming can be used to quantitatively explain the emergence of fractures in wet granular systems and, specifically, whether the onset of fracturing in our system arises from a jamming transition.
The jamming transition in a dry granular system occurs at a threshold packing fraction φ c when mechanical stability is achieved. For φ < φ c , the network of contact forces is constantly evolving and changing topology through particle rearrangement. For φ > φ c , in contrast, the force network locks in and its strength is enhanced through particle deformation [61,64]. Classic metrics that characterize the transition in frictionless systems are a discontinuous increase in the mean contact number Z, a rise in the mean isotropic stress P of the granular pack above its background value [61], or the emergence of a nonzero shear modulus [63].
We confirm that the behavior of our system responds in a manner consistent with a jamming transition. In particular, we compute at each stage of the granular pack deformation the Cauchy stress tensor for each particle in the system, where n c is the number of contacts for the particle. From the stress tensor we extract its isotropic component P = tr(σ i j ) and a measure of the shear stress, τ max = (σ max − σ min )/2, where σ max and σ min are the largest and smallest eigenvalues of σ i j , respectively. We observe that both quantities rise above a near-zero background as a function of the evolving mean packing fraction φ outside the central cavity [ Fig. 3(a)].
We determine the jamming transition φ c from the τ max profile as the intersection of two straight lines: one fitting the response of the background state, and one fitting a straight line to the asymptotic behavior in the highly compacted state [61,63] [Fig. 3(a), top inset]. For simulations with an initial packing density φ 0 = 0.77, the jamming transition occurs at a critical packing density φ c that takes increasing values (between 0.83 and 0.86) for increasing values of the contact angle (between θ = 75 • and 140 • ) [ Fig. 3(a)]. This result is consistent with our hypothesis of the emergence of fracturing being controlled by a jamming transition, in which the transition occurs earlier (at a smaller φ c ) in imbibition than in drainage. Previous studies of jamming transition in both frictionless [64,66,67] and frictional [61] systems show a power-law increase of the mean stress with packing fraction above jamming, P − P c ∼ (φ − φ c ) ψ , with an exponent slightly larger than 1, ψ ≈ 1.1. Our simulations for a wet granular system also show a power-law increase, with the exponent ψ in the range 1.06-1.39, larger values corresponding to drainage displacements and loose granular packs, and smaller values corresponding to imbibition displacements and dense granular packs [ Fig. 3(b), middle inset]. For our granular packings of finite μ = 0.3, Z c is expected to vary smoothly between Z c (μ = 0) = 4 and Z c (μ → ∞) → 3 [67,68]. Indeed, we find that Z c lies in the range of 3.49-3.96, and exhibits a power-law dependence with packing fraction above jamming, Fig. 3(c), bottom inset]. Earlier studies have found exponents at jam-ming in the vicinity of the jamming packing fraction and have shown that β ∼ 0.5 [61,64,66,67,69]. Here, we study the behavior of granular packs beyond the jamming transition, and therefore we conduct a correction-to-scaling analysis [70,71]: , with the leading correction-to-scaling exponent ω = 0.3 [70], and the prefactor a = 8.94 in the order of O(1), which validates the value of β obtained. The fact that fractures grow after the defined jamming transition φ c (as evidenced by a visual comparison of the interface morphology at jamming and at breakthrough [ Fig. 3(d)]) confirms our hypothesis that the onset of fractures emerges from a jamming transition.
A fundamental contribution to understanding jamming in (dry) granular systems was made in the form of a phase diagram that delineates the jammed state in the phase space of density, load, and temperature [72]. It shows that jamming can occur only at sufficiently high density, and that an increase in either load or temperature can unjam a system. We extend this description to wet granular systems by identifying quantities that determine the phase transition between jammed and unjammed states. We identify the packing fraction φ as the "density," and we posit that injection pressure P inj plays the role of the "load" during injection. Thus, we represent any generic evolution of our system as a trajectory in (P * inj , 1/φ) space (Fig. 4), where P inj is nondimensionalized by the characteristic capillary entry pressure in the system, γ /d.
Trajectories for regime I start with the prescribed φ 0 and move upwards in phase space as the granular pack is being compacted by the injected fluid. The injection pressure shows an initially decreasing and then increasing trend, as explained in Fig. 2(a). The transition from cavity expansion to fracturing corresponds to a transition from the unjammed state to the jammed state. We collect transition points φ c (shown as red markers in Fig. 4) for every simulation with a specific φ 0 and θ . These points collapse on a line in (P * inj , 1/φ) space, showing that under the same loading condition, the system jams at the same φ c , independently of θ or φ 0 . This transition line in the jamming phase diagram separates fundamentally different behaviors exhibited by our wet granular systems: fluidlike behavior (cavity expansion) in the unjammed state, and solidlike behavior (fracturing) in the jammed state (Fig. 4). This transition also helps explain the onset of fracturing: A larger energy input brought by the injection of a nonwetting fluid (larger value of the contact angle θ ) compacts the system to a denser state before jamming occurs, which, in turn, delays the onset of fracturing.
We also show in Fig. 4 the trajectories for regimes II, III, and IV. Frictional fingers (regime II) have a positive injection pressure. The trajectories corresponding to this regime move upwards in φ as the system is being compacted, with stick-slip fluctuations in P inj , but remain in the unjammed state for their entire evolution. Capillary invasion (regime III) occurs in an initially dense granular pack. The entire trajectory lies in the jammed state, with almost constant φ and stick-slip fluctuations in P inj . Capillary compaction (regime IV) occurs when the out-of-plane capillary pressure dominates and the granular pack is relatively loose initially. We calculate φ for the region inside the fluid-fluid interface. Since the negative dragging pressure is comparable for all our simulations in this regime (−50 to −10 Pa), the granular pack is compacted inwards up to approximately the same packing density (φ ≈ 0.83) above the jamming transition. At zero external load (P inj = 0), our system jams at the random close packing fraction φ c ≈ φ rcp ≈ 0.84 [73][74][75].
In summary, we have studied morphological transitions in granular packs as a result of capillary-dominated fluid-fluid displacement via a fully coupled model of twophase flow and grain mechanics. Simulations of fluid injection into a granular pack with different initial packing densities and substrate wettabilities have led us to uncovering four invasion regimes: cavity expansion and fracturing, frictional fingers, capillary invasion, and capillary compaction. In particular, we have identified the emergence of fracture, and its surprising and unexplored dependence on the system's wettability. We have shown that the onset of fracture can be explained as a jamming transition, as confirmed by the behavior of the classic metrics of jamming such as the mean isotropic stress. We have synthesized the system's response in the form of a phase diagram of jamming for wet granular media, on which the jamming transition for all different trajectories collapse on a single line in (P * inj , 1/φ) space, independently of the initial packing density φ 0 and contact angle θ . Due to the irreversible nature of friction during collective particle motion, pumping fluid back after injection-induced deformation will lead to a granular configuration very different from the initial packing, which lies outside the scope of this study.
Our study paves the way for understanding the impact of other key variables of a wet granular system, such as properties of the solid particles (rigidity, friction coefficient, cementation) or the fluid (viscosity contrast, capillary number). By tailoring the range of values of these variables, our analysis may provide fundamental insight on specific applications, from nanotechnology [76] to energy recovery [77], natural gas seeps [78,79], and geohazards [80,81]. This work was supported by the US Department of Energy (Grant No. DE-SC0018357).