Crossover from viscous fingering to fracturing in cohesive wet granular media: a photoporomechanics study

Abstract

We study fluid-induced deformation and fracture of cohesive granular media, and apply photoporomechanics to uncover the underpinning grain-scale mechanics. We fabricate photoelastic spherical particles of diameter d = 2 mm, and make a monolayer granular pack with tunable intergranular cohesion in a circular Hele–Shaw cell that is initially filled with viscous silicone oil. We inject water into the oil-filled photoelastic granular pack, varying the injection flow rate, defending-fluid viscosity, and intergranular cohesion. We find two different modes of fluid invasion: viscous fingering, and fracturing with leak-off of the injection fluid. We directly visualize the evolving effective stress field through the particles' photoelastic response, and discover a hoop effective stress region behind the water invasion front, where we observe tensile force chains in the circumferential direction. Outside the invasion front, we observe compressive force chains aligning in the radial direction. We conceptualize the system's behavior by means of a two-phase poroelastic continuum model. The model captures granular pack dilation and compaction with the boundary delineated by the invasion front, which explains the observed distinct alignments of the force chains. Finally, we rationalize the crossover from viscous fingering to fracturing by comparing the competing forces behind the process: viscous force from fluid injection that drives fractures, and intergranular cohesion and friction that resist fractures.

---

1. Introduction

Multiphase flow through granular and porous materials exhibits complex behavior, the understanding of which is critical in many natural and industrial processes, and the design of climate-change mitigation strategies. Examples include infiltration of water into the vadose zone,¹ methane migration in lake sediments,² hillslope infiltration and erosion after forest fires,³ growth and deformation of cells and tissues,⁴ shale gas production,⁵ and geological carbon dioxide storage.⁶ This complexity is linked to the interplay between multiphase flow and granular mechanics, which controls the morphological patterns, evolution, and function of a wide range of systems.⁷ In many granular-fluid systems, the strong coupling among viscous, capillary, and frictional forces leads to a wide range of patterns, including desiccation cracks,^8,9 labyrinth structures,¹⁰ granular fingers,^11–13 corals, and stick-slip bubbles.¹⁴ In the context of interfacial flows, fracture patterns have been observed in loose systems—such as particle rafts as a result of surfactant spreading^15,16—as well as dense systems—such as colloidal suspensions as a result of drying.^17,18

While fracturing during gas invasion in fluid-saturated media has been studied extensively in experiments^12,15–22 and simulations,^13,23–30 the underlying grain-scale mechanisms behind the morphodynamics and rheologies exhibited by deformable granular media remain poorly understood. To tackle this challenge, Meng et al.³¹ adopt a recently developed experimental technique, photoporomechanics,³² to directly visualize the evolving effective stress field in a fluid-filled cohesionless granular medium during fluid-induced fracturing. The effective stress field exhibits a surprising and heretofore unrecognized phenomenon: behind the propagating fracture tips, an effective stress shadow, where the intergranular stress is low and the granular pack exhibits undrained behavior, emerges and evolves as fractures propagate.

Here we aim to extend our previous work³¹ to cohesive granular media. The mechanical and fracture properties of cohesive granular media are of interest for many applications, including powder aggregation,^33,34 stimulation of hydrocarbon-bearing rock strata for oil and gas production,³⁵ preconditioning and cave inducement in mining,^36,37 and remediation of contaminated soil.³⁸ Similar hydraulic fractures manifest naturally at the geological scale, such as magma transport through dikes^39–43 and crack propagation at glacier beds.^44,45 Following the early work on cemented aggregates,^46–48 Hemmerle et al.⁴⁹ created a well-defined cohesive granular medium with tunable elasticity by mixing glass beads with curable polymer. Due to the huge stiffness contrast between polymer bridges (kPa–MPa) and glass beads (GPa), the mechanical response of the material is dominated by the deformation of the bridges rather than the deformation of the beads.^50,51 There is limited experimental study on weakly sintered or cemented materials⁵² with bonds that are of comparable stiffness with that of the grains.

In this paper, we study fracturing in cohesive wet granular media and extract the evolving effective stress field via photoporomechanics. By mixing photoelastic particles with curable polymer of the same stiffness, we create a well-defined cohesive granular medium with tunable tensile strength. We uncover two modes of water invasion under different injection flow rate, defending fluid viscosity, and tensile strength of the granular pack: viscous fingering, and fracturing with leak-off of the injection fluid. Behind the invasion front, the granular pack is dilated with tensile effective stress in the circumferential direction, while ahead of the invasion front the granular pack is compacted with compressive effective stress in the radial direction. We develop a two-phase poroelastic continuum model to explain the observed distinct force-chain alignments. Finally, we conclude that the competition of intergranular cohesion and friction against viscous force dictates the crossover from viscous fingering to fracturing regime.

2. Materials and methods

Following the fabrication process in ref. 32, we produce photoelastic spherical particles with diameter d = 1.98 mm (with 3.5% standard deviation) and bulk modulus K_s = 1.6 MPa. The fabrication process is similar to “squeeze casting” for metals, but for polyurethane in our case. The process produces soft polyurethane spheres with an amber color. To test their sliding frictional properties, we build a shear box apparatus as follows. We prepare a thin acrylic plate in the size of 6 cm × 6 cm × 1.6 mm and punch 11 × 11 holes with diameter 2 mm into it. We squeeze particles into the holes, making sure they are integrated into the plate and can not roll against it. The bottom surface for the sliding test is either made of glass or cured from polyurethane. We then put a confining weight on the top of the acrylic plate, which varies from 2 N to 8 N. We immerse particles in the silicone oil, attach the side surface of the acrylic plate to a spring scale, and drag the plate at a constant velocity 1 mm s⁻¹. The spring scale measures the frictional force occurring between the particles and the bottom surface. After dividing it by the confining weight, we obtain the friction coefficient. When particles are immersed in the silicone oil, the friction coefficient between particles is μ_p = 0.2 ± 0.06, and the friction coefficient between the particle and the glass plate is μ_w = 0.05 ± 0.02. To prepare a cohesive granular pack, we mix a total mass M_g of cured photoelastic particles with a mass M_l of uncured, liquid-form polyurethane. We set M_g = 9.2 g to generate granular packs at a fixed initial two-dimensional packing density, 0.83, which is close to the random close packing density. We cast the solid–liquid mixture into a monolayer of granular pack inside a circular Hele–Shaw cell. The added polyurethane is imbibed directly into the granular pack and forms polymer bridges between particles that solidify once cured. Before the experiments, we peel the cured monolayer of particles out of the Hele–Shaw cell, eliminating bonds between particles and plates. We define the polymer content C as the mass of added polyurethane divided by the mass of particles, which determines the size of polymer bridges and thus the tensile strength of interparticle bonds. Fig. 1 shows a monolayer of cohesive photoelastic granular pack at a polymer content C = 4.9%, above which pendular bridges begin to merge and form clusters.⁴⁹

Fig. 1 (a) Cohesive photoelastic granular pack at a polymer content C = 4.9%. (b) Close-up view of three particles connected by polymer bridges in the form of pendular rings.

We inject water into a monolayer of cohesive photoelastic particles saturated with silicone oil in a circular Hele–Shaw cell (Fig. 2). To observe the photoelastic response of the particles, we construct a darkfield circular polariscope by a white light panel together with left and right circular polarizers.⁵³ We clamp the Hele–Shaw cell vertically at a fixed height h = 1.92 mm with four internal spacers to achieve plane-strain conditions throughout the experiments. As the height of internal spacers is smaller than the particle diameter, the top plate applies vertical confinement on the top of the particles. To allow the fluids (but not the particles) to leave the cell, the disk is made slightly smaller than the interior of the cell (inner diameter L = 10.6 cm), resulting in a thin gap around the edge of the cell. A coaxial needle is inserted at the center of the granular pack for saturation, fluid injection and pore pressure measurement. We adopt a dual-camera system to record brightfield (camera A) and darkfield (camera B) videos. As water is injected into the cohesive granular pack, viscous forces from fluid injection promote the development of fractures, while intergranular cohesion within the granular pack resists it. To study these competing forces during the fluid invasion process, we vary the injection flow rate Q from 5 mL min⁻¹ to 220 mL min⁻¹, the silicone oil viscosity η from 30 kcSt, 100 kcSt, to 300 kcSt, and the polymer content C from 0% to 4.6% to stay in the pendular regime.

Fig. 2 Experimental setup to study fracturing in cohesive photoelastic granular media. (a) A monolayer of photoelastic particles (diameter d, polymer content C) is confined in a circular Hele–Shaw cell. The cell is uniformly clamped at a fixed height, h. Before the fracturing experiment, silicone oil (viscosity η) is slowly injected at the center of the cell with a coaxial needle to saturate the granular pack. After saturation, water is injected at a fixed flow rate Q with the injection pressure monitored by a pore-pressure sensor. A white-light panel, right and left circular polarizers form a darkfield circular polariscope. Brightfield and darkfield videos are captured by cameras placed underneath the cell. (b) Schematic of the circular Hele–Shaw cell (internal diameter L). A light panel, a polarizer, and a glass disk rest on top of the monolayer of photoelastic particles. The disk is slightly smaller than the cell to allow the fluids (but not particles) to leave the cell.

3. Representative experiments

In this section, we present two representative experiments with Q = 100 mL min⁻¹, η = 300 kcSt, C = 4.4% for viscous fingering, and C = 1.2% for fracturing with leak-off. For the dye color of the injected water, we need one that both visualizes the water invasion front in brightfield images and does not interfere with the photoelastic response in darkfield images. As a result, we dye the defending oil in black, and the invading water in light blue. By tracking the region in light blue color, we could easily identify the invading phase from brightfield images. To confirm this, we refer to the ESI† videos on the fluid morphology and the effective stress evolution for the two experiments.⁵⁴

We differentiate between viscous fingering and fracturing with leak-off regimes from the brightfield images. When water invades into the granular pack in viscous fingering patterns without noticeable relative motion between particles, the experiment is classified as viscous fingering (Fig. 3). When the injected water creates open channels with ensuing invasion into the pores, then the experiment is classified as fracturing with leak-off (Fig. 4). The darkfield images in Fig. 4 also confirm the formation of fractures where intergranular bonds exhibit strong photoelastic response and are torn apart under tension.

Fig. 3 For the viscous fingering experiment with Q = 100 mL min⁻¹, η = 300 kcSt and C = 4.4%, a sequence of snapshots shows the time evolution of (a) interface morphology from brightfield images, and (b) effective stress field from darkfield images. See ESI,† Video S1 for the evolution of the morphology and effective stress field in this regime.

Fig. 4 For the fracturing experiment with Q = 100 mL min⁻¹, η = 300 kcSt and C = 1.2%, a sequence of snapshots shows the time evolution of (a) interface morphology from brightfield images, and (b) effective stress field from darkfield images. See ESI,† Video S2 for the evolution of the morphology and effective stress field in this regime.

3.1 Viscous fingering

In Fig. 3, we show a sequence of snapshots for the viscous fingering experiment. We analyze the time evolution of the water–oil interface morphology from brightfield images, and the interparticle stresses of the granular pack from darkfield images. When particles have been strongly cemented initially, the injection pressure is insufficient to overcome the tensile strength of the intergranular bonds. In this regime, we observe patterns of viscous fingers without any significant relative motion between particles (Fig. 3(a)), and the intergranular bonds at finger tips endure tension without breakage (Fig. 3(b)).

3.2 Fracturing with leak-off

In Fig. 4, we show a sequence of snapshots for the fracturing experiment. The time evolution of the injection pressure P_inj is plotted in Fig. 5(a), which also indicates the times of the snapshots in Fig. 4. When particles have been densely packed initially, water firstly invades into the cohesive granular pack by expanding a small cavity around the injection port, with P_inj ramping up during this period. The pressure keeps building up until it becomes sufficient to overcome the tensile strength of the interparticle bonds; the point at which fractures emerge (in between Fig. 4(i)–(ii)). As injection continues, much of the injected water leaks off into the permeable granular media as the fractures propagate (Fig. 4(ii)–(v)). In this period of fracturing with leak-off, P_inj slightly drops from its peak value and reaches a plateau afterwards (Fig. 5(a)). In this regime, the effective stress field exhibits a surprising phenomenon: behind the water invasion front, a hoop effective stress region, where we observe tensile force chains in the circumferential direction, emerges and evolves as invasion front propagates. Ahead of the invasion front, we observe compressive effective stress in the radial direction (Fig. 4(b)). The phenomena regarding the effective stress (e.g. tensile hoop stress near the injection port) has been predicted by continuum theories, such as cavity expansion models for single-phase flow,^55,56 and tip asymptotics in fracture mechanics (Sections 2 and 3 in ref. 57). However, there is a lack of understanding of the effective stress evolution in a two-phase immiscible flow system, and our experiments visualize it for the first time.

Fig. 5 Modeling results for the fracturing experiment with Q = 100 mL min⁻¹, η = 300 kcSt and C = 1.2%. (a) Time evolution of the injection pressure P_inj. The solid curve represents the experimental measurement, and the dashed curve represents the model prediction. The cross markers indicate water breakthrough through the cell edge. The circular markers indicate times for the snapshots shown in Fig. 4–8: t = 0.2 s, 1.2 s, 2.2 s, 3.2 s, 4.2 s in sequence. Modeling results of the time evolution of (b) pore pressure p(r, t), and (c) water saturation S_w(r, t). We show the solution at 8 times, linearly spaced from t = 0.75 s (light gray) to t = 6 s (black).

4. Two-phase poroelastic continuum model

We model the immiscible injection of water into a cohesive granular pack saturated with silicone oil. To rationalize the experimental findings in Section 3, we develop a two-phase poroelastic continuum model focusing on the fracturing with leak-off regime. The wetting phase is the defending oil, and the nonwetting phase is the invading water. Under the experimental conditions explored, the modified capillary number²⁰ Ca* = ηQL/(γhd²) ∼ 10⁶, which indicates that viscous forces outweigh capillary forces so we can safely neglect capillary effects. In the following, we present the extension of Biot's theory⁵⁸ to two-phase flow.^59,60 In our model, we assume radial symmetry and small deformations.

4.1 Fluid flow equations

For the two-phase immiscible flow system, the conservation of fluid mass can be written as follows:where ϕ is the porosity and ρ_α and S_α are the density and saturation of the fluid phase α (water w or oil o), respectively. The phase velocity v_a is related to the Darcy flux q_α in a deforming medium by the following relation:where v_s is the velocity of the solid skeleton, k₀ is the intrinsic permeability of the granular pack, g is the gravity vector, and η_α, k_rα = k_rα(S_α) and p_α are the dynamic viscosity, relative permeability, and fluid pressure for phase α, respectively. Since capillary pressure is negligible here, p_c = p_w − p_o = 0, the two phases have the same fluid pressure p. The relative permeability functions are given as Corey-type power law functions:⁶¹where the fitting parameters are the critical water saturation for water to flow, S_wc = 0.2, the residual oil saturation, S_ro = 0.2, and the exponents a_w = 2 and a_o = 5.

Considering the mass conservation of the solid phase:

where ρ_s is the density of the solid constituents of the porous medium. Assuming isothermal conditions and using the equation of state for the solid, the following expression for the change in porosity is obtained:⁶²where b is the Biot coefficient of the saturated porous medium, and c_s is the compressibility of the solid phase. We use eqn (2), (4), and (5) to expand eqn (1) as follows:where ε_kk is the volumetric strain of the solid phase. The Biot modulus of the saturated granular pack, M, is related to fluid and rock properties as .⁶³ Adding eqn (6) for water and oil phases, and imposing that S_o + S_w ≡ 1 for the saturated granular pack, we obtain the pressure diffusion equation:where q_t is the total Darcy flux for water and oil phases, q_t = q_w + q_o.

4.2 Geomechanical equations

Under quasi-static conditions, the balance of linear momentum of the solid-fluid system states that:

∇·σ + ρ_bg = 0, (8)

where σ is the Cauchy total stress tensor, and , is the bulk density for the solid–fluid system. For axisymmetric deformation in plane-strain condition, the force balance equation becomes:Following,⁶³ the poroelasticity equation states thatwhere represents the identity matrix. Terzaghi's effective stress tensor σ′ is the portion of the stress supported through deformation of the solid skeleton, and where we adopt the convention of tension being positive. We adopt isotropic linear elastic theory for the granular pack; the constitutive equation for stress–strain is:where K, ν are the drained bulk modulus, and the drained Poisson ratio of the granular pack, respectively. The strain tensor is defined as , where u is the displacement vector. For the axisymmetric deformation in plane-strain condition, the strains are written as:Using eqn (10), (11) and (12), the force balance eqn (9) can be expressed as a function of radial displacement u_r(r,t) and pore pressure p(r,t).

4.3 Summary of governing equations

The model has three governing equations, two derived from conservation of fluid mass [eqn (7) for the water–oil fluid mixture and eqn (6) for the water phase] and one derived from conservation of linear momentum (eqn (9)). The model solves the time evolution of three unknowns: (1) pore pressure field p(r,t); (2) water saturation field S_w(r,t); and (3) radial displacement field u_r(r,t) of the cohesive granular pack, all of which are assumed to be radially symmetric. The governing equations are summarized and written in radial coordinates as follows:For the axisymmetric deformation in plane-strain condition, the stresses and strains are written in radial coordinates as:

ε_kk = ε_rr + ε_θθ + ε_zz. (19)

We initialize the model by specifying u_r(r, 0) = p(r, 0) = S_w(r, 0) = 0. As for the boundary conditions, the inner boundary is free to move, subject to injection pressure as the total stress:

σ_rr(r₀, t) = −p(r₀, t) = −P_inj(t), (20)

where r = r₀ is the radius of the injection port, and P_inj(t) is the injection pressure at time t. At the injection port, the total Darcy flux is the same as the Darcy flux of water. Since the injection system is composed of plastic syringe and tubing, we take the system compressibility into account for the inner flow boundary condition:where Q is the injection flow rate, and c_sys and V_sys are the compressibility and volume of the injection system, respectively. At the outer boundary r = R, the pressure is atmospheric, and particles have zero radial displacement:

p(R, t) = u_r(R, t) = 0. (22)

We now summarize the model in dimensionless form, denoting dimensionless quantities with a tilde. We adopt characteristic scales for length, time, stress/pressure, viscosity and permeability, non-dimensionalizing the governing equations via the scalingwhere is the characteristic poroelastic timescale. We can then rewrite eqn (13)–(15) in dimensionless form,where the dimensionless stresses are written in radial coordinates as:We initialize the model by specifying ũ_r(r̃, 0) = p̃(r̃, 0) = S_w(r̃, 0) = 0. The boundary conditions are as follows:where , . Both of these quantities compare the viscous pressure due to injection with the Biot modulus of the granular pack.

4.4 Model parameters

The four poroelastic constants in the model are the drained bulk modulus K, the drained Poisson ratio ν, the Biot coefficient b, and the Biot modulus M of the granular pack. We obtained the drained and undrained bulk modulus K, K_u from a separate consolidation experiment.³² We calculate the Biot coefficient from the relationship ,⁶³ and then obtain the Biot modulus via.⁶⁴ To obtain the drained Poisson ratio of the granular pack, we build a discrete element model and conduct a biaxial test.⁶⁵ The permeability of the granular pack, k, is measured experimentally during the initial oil saturation process. A summary of the modeling parameters is given in Table 1.

Table 1 Modeling parameters for the two-phase poroelastic continuum model

Symbol	Value	Unit	Variable
r 0	4	mm	Injection port radius
R	5.3	cm	Hele–Shaw cell radius
d	2	mm	Diameter of the photoelastic particles
h	1.92	mm	Height of the Hele–Shaw cell
Q	100	mL min^−1	Water injection flow rate
c sys	6 × 10^−8	Pa^−1	Injection system compressibility
V sys	30	mL	Injection system volume
K	200	kPa	Drained bulk modulus of the pack
K u	1.35	MPa	Undrained bulk modulus of the pack
ν	0.4		Drained Poisson ratio of the pack
b	0.88		Biot coefficient of the pack
M	1.49	MPa	Biot modulus of the pack
η w	0.001	Pa s	Injecting water viscosity
η o	291.3	Pa s	Defending silicone oil viscosity
ϕ	0.4		Porosity of the pack
k 0	10^−10	m^2	Intrinsic permeability of the pack

---

4.5 Numerical implementation

We use a finite volume numerical scheme to solve the three coupled governing equations (eqn (13)–(15)). We place the pressure and saturation unknowns (p(r, t), S_w(r, t)) at volume centers, and displacement unknowns (u_r(r, t)) at nodes. We partition the coupled problem and solve two sub-problems sequentially: the coupled flow and mechanics, and the transport of water saturation. We first fix the water saturation, and solve the coupled flow and mechanics equations (eqn (13) and (15)) simultaneously using implicit time discretization. Then we solve the water transport equation (eqn (14)) with prescribed pressure and displacement fields.

5. Results and discussion

5.1 Pore pressure and water saturation

We compare the experimental and modeling results of the time evolution of the injection pressure P_inj for the case with Q = 100 mL min⁻¹, η = 300 kcSt and C = 1.2% (Fig. 5(a)). By taking the injection system compressibility into account, the model captures the initial P_inj ramp-up measured in the experiment (t = 0–0.3 s). Before t = 3.5 s, P_inj keeps increasing, with the diffusion of pore pressure (Fig. 5(b)) and propagation of water invasion front (Fig. 5(c)). The transient pressure response comes from the compressibility of the granular pack, the timescale of which is , where η′ and k′ are the effective viscosity and permeability of the pore fluid: a water-oil mixture. As the pore pressure diffuses to the cell boundary, P_inj approaches its steady state value, .

The cross markers in Fig. 5(a) represent the moment when water reaches the cell boundary for the experiment and the model. The breakthrough predicted by the model is faster than that of the experiment by around 1.2 s. The observation that the water invasion front propagates faster in the model hints at an overestimation of the Biot modulus M; in other words, the model underestimates the granular pack compressibility/storativity. It reveals two underlying model limitations: (1) the storativity in the model, , is assumed to be a constant without spatiotemporal variations, which in the experiment increases with porosity in the region where the granular pack dilates; and (2) by assuming linear elastic granular packs with small deformations, the model cannot capture the significant increase in storativity arising from the opening of fracture, where the porosity of the granular pack locally increases to 1.

Solving the water transport equation (eqn (14)), we obtain the time evolution of the water saturation field (Fig. 5(c)). The water invasion front propagates with the injection until its breakthrough at t = 2.6 s. After breakthrough, the radial profile of water saturation is nonmonotonic, exhibiting an increase of S_w with r and then a decrease. The position of the local maximum of the saturation profile moves towards the center of the cell as time evolves, and the value of the maximum saturation increases with time. This unusual behavior of water saturation contrasts that of fluid–fluid displacement in a rigid porous medium^66,67 and highlights the strong coupling between fluid flow and medium deformation in our system.

5.2 Displacement and volumetric strain

To probe into the granular mechanics behind the fracturing experiment in Fig. 4, we first measure the internal deformation of the pack via particle tracking, which provides a direct measure of the displacement field. 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 u_i = (x_i − X_i, y_i − Y_i), with magnitude and radial component . The deformation is primarily radial because of the axisymmetric boundary conditions, so we focus on u_r.

Fig. 6 shows a sequence of snapshots of the experimental radial displacement field, corresponding to t = 0.2 s, 1.2 s, 2.2 s, 3.2 s, 4.2 s sequentially. We find that the radial displacement is large near the injection port and fades to zero at the rigid outer edge, with a petal-like mesoscale structure as reported by MacMinn et al.⁶⁸ and Zhang and Huang⁶⁹ (Fig. 6(a)). The radial displacements of particles are plotted as black dots in Fig. 6(b), and the red dashed line is the prediction from the continuum model. As P_inj increases from snapshots (i) to (iii), particles move radially outwards. From snapshots (iii) to (v), P_inj reaches a plateau, and particles relax and recover part of the deformation. The model captures the general trends in particle displacement behavior, with the notable exception of the compaction front near r ∼ 0.5R between times (iii) and (iv). Between this time period, the experimental data shows that particles with r < 0.5R are compacted to a similar extent, as evidenced by their similar u_r values, which we refer to as a compaction front. The model underestimates the displacements there due to our assumption of linear elastic behavior: it cannot capture the plasticity-induced compaction front brought by bond breakage and particle rearrangements. As a result, the model fails to capture the compaction front exhibited in the experiment, which we define as the plasticity-induced compaction front.

Fig. 6 For the fracturing experiment with Q = 100 mL min⁻¹, η = 300 kcSt and C = 1.2%, we analyze the sequence of snapshots shown in Fig. 4(i)–(v) corresponding to t = 0.2 s, 1.2 s, 2.2 s, 3.2 s, 4.2 s, respectively. The sequence of snapshots shows the time evolution of (a) experimental radial displacement field, and (b) radial displacement of the particles (black dots) compared with the continuum model prediction (dashed line).

We use the particle positions to calculate a best-fit local strain field. At time t during the water injection, we compute the closest possible approximation to a local strain tensor ε in the neighborhood of any particle with a sampling radius r_s = 3d.^68,70 The local strain for the particle is determined by minimizing the mean-square difference D² between the actual displacements of the neighboring particles relative to the central one and the relative displacements that they would have if they were in a region of uniform strain ε_ij. That is, we define

where the index n runs over the particles within the sampling radius and n = 0 for the reference particle. We then compute ε for the reference particle at time t that minimizes D²(t). With this method, we obtain the local strain tensors for all particles in the granular pack.

We present a sequence of snapshots of the volumetric strain field in Fig. 7. The granular pack dilates (positive ε_kk) in the water-invaded region, and compacts (negative ε_kk) in the oil-saturated region (Fig. 7(a)). Such injection-induced dilation has also been reported for cohesionless granular packs and cohesive poroelastic cylinders.^55,68,71Fig. 7(b) shows that the model captures the granular dilation and compaction, but cannot account for the plastic dilation near fractures brought by bond breakage and particle rearrangements.

Fig. 7 For the fracturing experiment with Q = 100 mL min⁻¹, η = 300 kcSt and C = 1.2%, we analyze the sequence of snapshots shown in Fig. 4(i)–(v) corresponding to t = 0.2 s, 1.2 s, 2.2 s, 3.2 s, 4.2 s, respectively. The sequence of snapshots shows the time evolution of (a) experimental volumetric strain field, and (b) volumetric strain of the particles (black dots) compared with the continuum model prediction (dashed line).

The injection-induced deformation also feeds back to the fluid flow, as evidenced by the observed nonmonotonic water saturation curves (Fig. 5(c)). The granular pack dilation near the injection port increases the local porosity, and results in a smaller value of S_w. As encoded in eqn (14), the coupling between fluid flow and medium deformation becomes strong when the deformation term is comparable to the flow term, .

5.3 Effective stress

The photoelastic response offers a unique opportunity to gain additional understanding of the coupled pore-scale flow and particle mechanics during fluid-induced fracturing of the cohesive granular pack. To interpret the photoelastic response, we rely on the results of calibration experiments,³² which have shown that, for the range of interparticle forces expected in our fracturing experiments, the relation between light intensity and force is monotonically increasing and approximately linear. From two-dimensional photoelasticity theory,⁷² the stress-optic law states that in this “first-order” region, the photoelastic response is approximated to be linearly proportional to the principal effective stress difference with a constant coefficient: , where and are the maximum and minimum principal effective stresses, respectively.

To quantify the photoelastic response into the principal effective stress difference, we conduct a separate calibration experiment to obtain the coefficient, F = 0.29 from the blue channel light intensity (see Appendix). To differentiate the direction of force chains, we set an ad hoc sign convention for the principal effective stress difference δσ′, which should otherwise always be positive, as follows: δσ′ is positive for tensile force chains in circumferential/hoop direction, and negative for compressive force chains in radial direction. After converting I into δσ′, and assigning its sign from the force chain direction, we present a sequence of snapshots of the effective stress field (Fig. 8(a)). Behind the water invasion front, a hoop effective stress region, where we observe tensile force chains in the circumferential direction, emerges and evolves as the invasion front propagates. Ahead of the invasion front, we observe radial compaction of the granular pack.

Fig. 8 For the fracturing experiment with Q = 100 mL min⁻¹, η = 300 kcSt and C = 1.2%, we analyze the sequence of snapshots shown in Fig. 4(i)–(v) corresponding to t = 0.2 s, 1.2 s, 2.2 s, 3.2 s, 4.2 s, respectively. The sequence of snapshots shows the time evolution of (a) experimental effective stress field, and (b) the radial distribution of the averaged effective stress (solid line) compared with the continuum model prediction (dashed line). To differentiate the direction of force chains, we set a sign convention manually to the principal effective stress difference δσ′, which should otherwise always be positive, as follows: δσ′ is positive for tensile force chains in circumferential/hoop direction, and negative for compressive force chains in radial direction.

In the model, we found that always holds, where and are the hoop and radial components of the effective stress, respectively. As the force chain direction aligns with the effective stress direction with larger absolute magnitude, we calculate δσ′ numerically with the aforementioned sign convention as follows:

We compare the experimental and numerical radial distribution of δσ′ in Fig. 8(b). The model captures the hoop effective stress region and radial compaction delineated by the invasion front. As mentioned in our previous discussion on the displacement field, the model cannot capture the plasticity-induced compaction front, resulting in an underestimation of compressive effective stress between times (iii) and (iv).

5.4 Phase diagram of fluid invasion patterns in cohesive granular media

We observe two invasion patterns when varying the experimental parameters η, Q, and C: (I) pore invasion in the form of immiscible viscous fingering, and (II) fracturing with leak-off of the injection fluid. In a granular medium, fractures open when forces exerted by the fluids exceed the mechanical forces that resist particle rearrangements. Here the competing forces are the viscous force that drives fractures, and intergranular cohesion and friction that resist fractures. For a fixed domain geometry and granular medium (particle size and packing fraction), the viscous force is expected to increase with the product of the fluid viscosity η and the injection rate Q. We use a dimensionless flow rate to characterize the viscous force. The resisting force is expected to have a friction-dependent component that will be constant across our experiments, and a cohesion-dependent component that will increase with the polymer content C. We use a dimensionless tensile strength to characterize the resisting force. Thus, we plot an empirical phase diagram of all our experiments, indicating whether they are either “fracturing” or “viscous fingering” (not fracturing) on the axes q̃ vs. C (Fig. 9(a)). This empirical plot shows a transition from viscous fingering at low ηQ and high C to fracturing at high ηQ and low C.

Fig. 9 Phase diagrams of fluid–fluid displacement patterns in the experiments. Diagram (a) shows the invasion patterns for all experiments, ranging in oil viscosity η from 30 kcSt, 100 kcSt, to 300 kcSt, water injection rate Q from 5 mL min⁻¹ to 220 mL min⁻¹, and polymer content C from 0 to 4.6%. Diagram (b) shows the modeling prediction of the dimensionless maximum effective hoop stress at the injection port, , as a function of the dimensionless flow rate, . The dashed line shows the fitted power law, . Diagram (c) shows the dimensionless tensile strength of the granular pack against fracturing, , increases with polymer content C. The dashed line shows the fitted linear relationship, .

In the fracturing experiment (Section 3.2), the photoelastic response reveals that fractures initiate as tensile cracks near the injection port, where intergranular bonds break under tensile stress in the circumferential direction. Shear failure also occurs during fracture propagation, as evidenced by the classic slip line fracture pattern. To rationalize the crossover from viscous fingering to fracturing regimes quantitatively, we focus on the fracture initiation and assume the tensile failure mode. We adopt a fracturing criterion for cohesive granular media: the maximum hoop effective stress should exceed the tensile strength between particles to break interparticle bonds and generate fractures. To theoretically predict , we run the model with different flow conditions, and obtain at the injection port. We then obtain the dimensionless maximum hoop effective stress by . Fig. 9(b) shows that increases with q̃ approximately as a power law, .

To construct the relationship between and C, we record the injection pressure at the onset of fracturing when interparticle bonds break as P^frac_inj. We obtain the dimensionless tensile strength , and find a linear relationship, (Fig. 9(c)). It does not pass through the origin because of the frictional resistance against fracturing for a cohesionless granular pack. After entering these relationships into the fracturing criterion, we obtain a condition involving q̃ and C for the transition into the fracturing regime, q̃ ≥ (38.67C + 0.51)^6.0. The theoretical prediction on the crossover from viscous fingering to fracturing regime agrees well with the experimental results (Fig. 10).

Fig. 10 Phase diagram of fluid–fluid displacement patterns in the experiments. The dashed line is the theoretical prediction on the crossover from viscous fingering to fracturing regime, q̃ = (38.67C + 0.51)^6.0.

6. Conclusions

In summary, we have studied the morphology and rheology of injection-induced fracturing in cohesive wet granular packs via a recently developed experimental technique, photoporomechanics, which extends photoelasticity to granular systems with a fluid-filled connected pore space.³² Experiments of water injection into cohesive photoelastic granular packs with different tensile strength, injection flow rate, and defending fluid viscosity have led us to uncover two invasion regimes: viscous fingering, and fracturing with leak-off of the injection fluid. Contrary to the observed effective stress shadow for cohesionless granular packs,³¹ here we discover a hoop effective stress region behind the water invasion front. We developed a two-phase poroelastic continuum model that captures the transient pressure response arising from the granular pack compressibility. Behind the water invasion front, the granular pack is dilated with tensile force chains in the circumferential direction. Ahead of the water invasion front, the granular pack is compacted with compressive force chains in the radial direction. Finally, we rationalize the crossover from viscous fingering to fracturing across our suite of experiments by comparing the competing forces behind the process: viscous force from fluid injection that drives fractures, and intergranular cohesion and friction that resist fractures.

The developed two-phase continuum model assumes linear elasticity, which is insufficient to capture bond breakage and particle rearrangements. In spite of its limitations, our minimal-ingredients model still sheds insight and explains some of the key features observed in the experiments. The model reveals that the transient pressure response comes from the compressibility of the granular pack. It also captures granular pack dilation and compaction with the boundary delineated by the invasion front, which explains the observed distinct alignments of the force chains. Furthermore, the model predicts the injection-induced hoop stress at the injection port where tensile cracks emerge, which is the key to rationalizing the crossover from viscous fingering to fracturing regimes quantitatively.

An interesting next step would be to account for these irreversible processes by means of a poroelastoplastic or poroviscoplastic model, possibly in large deformations to reflect the substantial variations in porosity during the fluid injection. Then the poroelastic constants could be taken to be porosity-dependent. One could start with extending previous work from Auton and MacMinn⁵⁶ to two-phase flow. To gain more insights on the fluid-induced fracturing, the radially symmetric model could be extended to a two-dimensional model that takes fracture morphology into account. Motivated by our experiments, Guevel et al.⁷³ develop a Darcy–Cahn–Hilliard model coupled with damage to describe multiphase-flow and fluid-driven fracturing in porous media. The model adopts a double phase-field approach, regularizing both cracks and fluid–fluid interfaces. The damage model allows for control over both nucleation and crack growth, and successfully recovers the flow regime transition from fingering to fracturing with leak-off observed in our experiments. Lastly, by adding capillarity in the fluid flow equations, the model would be able to explore a wider range in η and Q, and possibly explains more invasion regimes, such as capillary fingering and fracturing.

Our study paves the way for understanding the mechanical and fracture properties of cohesive porous materials that are of interest for applications in various fields of research and industry, such as rock mechanics,^46,74,75 the fracture of concrete and biomaterials,^76,77 and geoengineering.⁷⁸

Conflicts of interest

There are no conflicts to declare.

Appendix: Calibration experiment for photoelastic response

The stress-optic law states that in the first order, the photoelastic response is approximated to be linearly proportional to the principal effective stress difference with a constant coefficient: I = Fδσ′, where δσ′ is the principal effective stresses difference.⁷² To obtain the coefficient F, we conduct a calibration experiment in the same Hele–Shaw cell where we conduct the fracturing experiments.

We prepare a monolayer of photoelastic particles at a polymer content C = 3%. The particle diameter and initial packing density are the same as in the fracturing experiments. We saturate the granular pack with silicone oil of viscosity 5 cSt, which lubricates the particle–particle and particle-wall contacts. After saturation, we slowly inject water at Q = 2 mL min⁻¹ into a sealed, elastic membrane that is connected to the injection port, and we monitor the injection pressure during injection. As injection proceeds, P_inj increases and drives the outward compaction of the granular pack quasi-statically. The membrane expands in size without any water leakage. We present a sequence of snapshots of the blue-channel light intensity field from darkfield images (Fig. 11(a)). The injection takes place under drained conditions, where the pressure in the defending fluid has time to fully dissipate, and the solid skeleton takes all the load from the water pressure at the inner boundary. The process is the same as a classical linear elastostatic example: a cylindrical vessel subject to an internal pressure and fixed outer boundary.⁷⁹ For any specific P_inj and size of the inner boundary, we obtain the theoretical prediction on δσ′, which helps us to calibrate the conversion factor F between experimental light intensity (Fig. 11b)) and effective stress difference (Fig. 11(c)). The calibration shows that F = 0.29 under our experimental conditions.

Fig. 11 For the calibration experiment with increasing water injection pressure, a sequence of snapshots shows the time evolution of (a) experimental light intensity field from the blue channel, (b) the radial distribution of the averaged light intensity I(r, t), and (c) the radial distribution of the averaged effective stress difference (solid line) compared with the continuum model prediction (dashed line). The conversion factor between light intensity and effective stress difference is calibrated to be .

Acknowledgements

We thank Chris MacMinn (University of Oxford), John Dolbow (Duke University), Alex Guevel (Duke University) and Ken Kamrin (MIT) for helpful discussions. We acknowledge funding from the U.S. National Science Foundation (grant no. CMMI-1933416).

Notes and references

1. D. E. Hill and J. Y. Parlange, Soil Sci. Soc. Am. J., 1972, 36, 697–702.
2. B. P. Scandella, L. Pillsbury, T. Weber, C. Ruppel, H. F. Hemond and R. Juanes, Geophys. Res. Lett., 2016, 43, 4374–4381.
3. J. Mataix-Solera, V. Arcenegui, N. Tessler, R. Zornoza, L. Wittenberg, C. Martnez, P. Caselles, A. Pérez-Bejarano, D. Malkinson and M. M. Jordán, Catena, 2013, 108, 6–13.
4. G. T. Charras, J. C. Yarrow, M. A. Horton, L. Mahadevan and T. J. Mitchison, Nature, 2005, 435, 365–369.
5. T. W. Patzek, F. Male and M. Marder, Proc. Natl. Acad. Sci. U.S.A., 2013, 110, 19731–19736.
6. M. L. Szulczewski, C. W. MacMinn, H. J. Herzog and R. Juanes, Proc. Natl. Acad. Sci. U.S.A., 2012, 109, 5185–5189.
7. R. Juanes, Y. Meng and B. K. Primkulov, Phys. Rev. Fluids, 2020, 5, 110516.
8. A. Groisman and E. Kaplan, Europhys. Lett., 1994, 25, 415–420.
9. H. Shin and J. C. Santamarina, Geotechnique, 2011, 61, 961–972.
10. B. Sandnes, H. Knudsen, K. J. Måløy and E. G. Flekkøy, Phys. Rev. Lett., 2007, 99, 038001.
11. X. Cheng, L. Xu, A. Patterson, H. M. Jaeger and S. R. Nagel, Nat. Phys., 2008, 4, 234.
12. H. Huang, F. Zhang, P. Callahan and J. Ayoub, Phys. Rev. Lett., 2012, 108, 258001.
13. F. Zhang, B. Damjanac and H. Huang, J. Geophys. Res. Solid Earth, 2013, 118, 2703–2722.
14. B. Sandnes, E. Flekkøy, H. Knudsen, K. Måløy and H. See, Nat. Commun., 2011, 2, 288.
15. C. Peco, W. Chen, Y. Liu, M. M. Bandi, J. E. Dolbow and E. Fried, Soft Matter, 2017, 13, 5832–5841.
16. D. Vella, P. Aussillous and L. Mahadevan, Europhys. Lett., 2004, 68, 212.
17. E. R. Dufresne, E. I. Corwin, N. A. Greenblatt, J. Ashmore, D. Y. Wang, A. D. Dinsmore, J. X. Cheng, X. S. Xie, J. W. Hutchinson and D. A. Weitz, Phys. Rev. Lett., 2003, 91, 224501.
18. L. Goehring, W. J. Clegg and A. F. Routh, Phys. Rev. Lett., 2013, 110, 024301.
19. H. Shin and J. C. Santamarina, Earth Planet. Sci. Lett., 2010, 299, 180–189.
20. R. Holtzman, M. L. Szulczewski and R. Juanes, Phys. Rev. Lett., 2012, 108, 264504.
21. J. M. Campbell, D. Ozturk and B. Sandnes, Phys. Rev. Appl., 2017, 8, 064029.
22. Z. Sun and J. C. Santamarina, J. Geophys. Res. Solid Earth, 2019, 124, 2274–2285.
23. A. K. Jain and R. Juanes, J. Geophys. Res. Solid Earth, 2009, 114, B08101.
24. R. Holtzman and R. Juanes, Phys. Rev. E, 2010, 82, 046305.
25. B. Carrier and S. Granet, Eng. Fract. Mech., 2012, 79, 312–328.
26. B. Lecampion and J. Desroches, J. Mech. Phys. Solids, 2015, 82, 235–258.
27. A. Mikelic, M. F. Wheeler and T. Wick, Multiscale Model. Simul., 2015, 13, 367–398.
28. D. Santillán, R. Juanes and L. Cueto-Felgueroso, J. Geophys. Res. Solid Earth, 2018, 123, 2127–2155.
29. Y. Meng, B. K. Primkulov, Z. Yang, C. Y. Kwok and R. Juanes, Phys. Rev. Res., 2020, 2, 022012.
30. F. J. Carrillo and I. C. Bourg, Phys. Rev. E, 2021, 103, 063106.
31. Y. Meng, W. Li and R. Juanes, Phys. Rev. Appl., 2022, 18, 064081.
32. W. Li, Y. Meng, B. K. Primkulov and R. Juanes, Phys. Rev. Appl., 2021, 16, 024043.
33. W. Pietsch, E. Hoffman and H. Rumpf, Ind. Eng. Chem. Prod. Res. Dev., 1969, 8, 58–62.
34. K. Kendall and C. Stainton, Powder Technol., 2001, 121, 223–229.
35. M. J. Economides, K. G. Nolteet al., Reservoir Stimulation, Prentice Hall Englewood Cliffs, NJ, 1989, vol. 2.
36. A. van As and R. G. Jeffrey, Proceedings of the 5th North American Rock Mech. Symposium and the 17th Tunneling Association of Canada Conference (NARMSTAC conference), 2002, pp. 7–10.
37. R. G. Jeffrey, Z. Chen, K. W. Mills and S. Pegg, ISRM International Conference for Effective and Sustainable Hydraulic Fracturing, 2013.
38. L. C. Murdoch, J. Geotech. Geoenviron. Eng., 2002, 128, 488–495.
39. D. A. Spence, P. W. Sharp and D. L. Turcotte, J. Fluid Mech., 1987, 174, 135–153.
40. J. R. Lister and R. C. Kerr, J. Geophys. Res. Solid Earth, 1991, 96, 10049–10077.
41. A. M. Rubin, Annu. Rev. Earth Planet. Sci., 1995, 23, 287–336.
42. S. M. Roper and J. R. Lister, J. Fluid Mech., 2005, 536, 79–98.
43. S. M. Roper and J. R. Lister, J. Fluid Mech., 2007, 580, 359–380.
44. V. C. Tsai and J. R. Rice, J. Geophys. Res., 2010, 115, F03007.
45. C. Y. Lai, J. Kingslake, M. G. Wearing, P. C. Chen, P. Gentine, H. Li, J. J. Spergel and J. M. van Wessem, Nature, 2020, 584, 574–578.
46. J. Dvorkin, G. Mavko and A. Nur, Mech. Mater., 1991, 12, 207–217.
47. J. Dvorkin, A. Nur and H. Yin, Mech. Mater., 1994, 18, 351–366.
48. J. Dvorkin, J. Berryman and A. Nur, Mech. Mater., 1999, 31, 461–469.
49. A. Hemmerle, M. Schröter and L. Goehring, Sci. Rep., 2016, 6, 35650.
50. A. Schmeink, L. Goehring and A. Hemmerle, Soft Matter, 2017, 13, 1040–1047.
51. A. Hemmerle, Y. Yamaguchi, M. Makowski, O. Bäumchen and L. Goehring, Soft Matter, 2021, 17, 5806–5814.
52. R. Affes, J. Y. Delenne, Y. Monerie, F. Radjai and V. Topin, Eur. Phys. J. E, 2012, 35, 117.
53. K. E. Daniels, J. E. Kollmer and J. G. Puckett, Rev. Sci. Instrum., 2017, 88, 051808.
54. See supplementary videos (ESI†).
55. L. C. Auton and C. W. MacMinn, Proc. R. Soc. A: Math. Phys. Eng. Sci., 2018, 474, 20180284.
56. L. C. Auton and C. W. MacMinn, J. Mech. Phys. Solids, 2019, 132, 103690.
57. E. Detournay, Annu. Rev. Fluid Mech., 2016, 48, 311–339.
58. M. A. Biot, J. Appl. Phys., 1941, 12, 155–164.
59. B. Jha and R. Juanes, Water Resour. Res., 2014, 50, 3776–3808.
60. T. I. Bjørnarå, J. M. Nordbotten and J. Park, Water Resour. Res., 2016, 52, 1398–1417.
61. R. H. Brooks, Hydraulic properties of porous media, Colorado State University, 1965.
62. R. W. Lewis and B. A. Schrefler, The finite element method in the static and dynamic deformation and consolidation of porous media, John Wiley & Sons, 1998.
63. O. Coussy, Mechanics of porous continua, Wiley, 1995.
64. H. F. Wang, Theory of linear poroelasticity with applications to geomechanics and hydrogeology, Princeton University Press, 2000, vol. 2.
65. ITASCA, PFC2D, v3.1-Theory and Background, Itasca Consulting Group, Inc., Minneapolis, MN, 2004.
66. S. E. Buckley and M. C. Leverett, Trans. AIME, 1942, 146, 107–116.
67. M. J. Blunt, Multiphase flow in permeable media: A pore-scale perspective, Cambridge University Press, 2017.
68. C. W. MacMinn, E. R. Dufresne and J. S. Wettlaufer, Phys. Rev. X, 2015, 5, 011020.
69. F. Zhang and H. Huang, 45th US Rock Mechanics/Geomechanics Symposium, 2011.
70. M. L. Falk and J. S. Langer, Phys. Rev. E, 1998, 57, 7192.
71. L. C. Auton and C. W. MacMinn, Proc. R. Soc. A: Math. Phys. Eng. Sci., 2017, 473, 20160753.
72. M. Frocht, Photoelasticity, John Wiley & Sons, 1941.
73. A. Guével, Y. Meng, C. Peco, R. Juanes and J. E. Dolbow, 2023, arXiv, preprint, arXiv:2306.16930 DOI:10.48550/arXiv.2306.16930.
74. J. C. Jaeger, N. G. W. Cook and R. Zimmerman, Fundamentals of Rock Mechanics, John Wiley & Sons, 2009.
75. R. Holtzman, Int. J. Numer. Anal. Methods Geomech., 2012, 36, 944–958.
76. O. Buyukozturk and B. Hearing, Int. J. Solids Struct., 1998, 35, 4055–4066.
77. V. Topin, F. Radjai, J. Y. Delenne, A. Sadoudi and F. Mabille, J. Cereal Sci., 2008, 47, 347–356.
78. D. L. Turcotte, E. M. Moores and J. B. Rundle, Phys. Today, 2014, 67, 34.
79. L. Anand and S. Govindjee, Continuum mechanics of solids, Oxford University Press, 2020.

---

Footnotes

† Electronic supplementary information (ESI) available: 2 Movie files and a text document with captions for ESI files. See DOI: https://doi.org/10.1039/d3sm00897e

‡ Present address: Stanford University, Stanford, CA, USA.

§ Present address: Stony Brook University, Stony Brook, NY, USA.

---