Capillary and Viscous Fracturing During Drainage in Porous Media

Detailed understanding of the couplings between fluid flow and solid deformation in porous media is crucial for the development of novel technologies relating to a wide range of geological and biological processes. A particularly challenging phenomenon that emerges from these couplings is the transition from fluid invasion to fracturing during multiphase flow. Previous studies have shown that this transition is highly sensitive to fluid flow rate, capillarity, and the structural properties of the porous medium. However, a comprehensive characterization of the relevant fluid flow and material failure regimes does not exist. Here, we used our newly developed Multiphase Darcy-Brinkman-Biot framework to examine the transition from drainage to material failure during viscously-stable multiphase flow in soft porous media in a broad range of flow, wettability, and solid rheology conditions. We demonstrate the existence of three distinct material failure regimes controlled by non-dimensional numbers that quantify the balance of viscous, capillary, and structural forces in the porous medium.


I. INTRODUCTION
Multiphase flow in deformable porous media is a ubiquitous phenomenon in natural and engineered systems that underlies key processes in water and energy resource engineering and materials science, including membrane filtration, soil wetting/drying, unconventional hydrocarbon recovery, and geologic carbon sequestration [1][2][3]. A key obstacle to more accurate representations of this phenomenon is our limited understanding of the transition from fluid invasion to flow-induced fracturing, i.e., material failure caused by multiphase flow. In large part, this limitation is caused by a lack of computational approaches capable of representing multiphase flow in fractured deformable porous media.
Previous work on multiphase flow within static porous media is extensive and includes detailed examinations of the influence of wettability, viscosity, and flow rate on flow in unsaturated porous media at multiple scales. In particular, existing studies have demonstrated how capillary forces give rise to differences between drainage and imbibition [4]; how the ratio of fluid viscosities controls the stability of the invading fluid front [5][6][7]; and how the magnitude of the capillary number delineates distinct flow regimes [8,9]. Each of the aforementioned controls is highly dependent on the system of interest. This complicates efforts to develop general relative permeability and capillary pressure models that apply to most systems of interest [10][11][12][13].
Flow of a single fluid phase through deformable porous media also has been studied in depth. Numerical modeling studies are largely based on the work of Biot and Terzaghi [14,15] and have been used to reproduce the behavior of arteries, boreholes, swelling clays, and gels [16][17][18][19]. In the last decade, fundamental studies have generated detailed information on the dynamics that arise from fluid-solid couplings beyond the ideal poroelastic regime, including fracturing, granular fingering, and frictional fingering [20][21][22]. In particular, these studies have shown that the main parameters controlling the deformation of a porous solid by single phase flow are the material softness and the magnitude of the fluid-solid momentum transfer [21].
The study of multiphase flow in a deformable porous medium is inherently more complex than the problems outlined above, as it requires simultaneous consideration of capillarity, wetting dynamics, fluid rheology, and solid deformation. Deformation modes associated with material failure (i.e., multiphase fracturing) are particularly challenging as they require simultaneous representation of multiphase flow in fractures and in the surrounding porous matrix. The existing detailed examinations of this phenomenon have focused exclusively on granular systems. Notably, Holtzman & Juanes [23,24] used experiments and discrete element models to demonstrate that the transitions between capillary fingering, viscous fingering, and fracturing during multiphase flow in granular media reflect two non-dimensional numbers: a fracturing number (ratio of fluid driving force to solid cohesive force) and a modified capillary number (the ratio between viscous and capillary pressure drops). Other discrete element approaches have shown that fracturing is highly dependent on the invading fluid's capillary entry pressure [25,26]. However, it is not clear how these conclusions translate to continuous non-granular systems.
To the best of our knowledge, no experimental or numerical investigation has simultaneously explored the effects of flow rate, wettability, and deformability during multiphase flow in deformable porous media at the continuum scale and identified the controlling parameters that relate all three properties within a single phase diagram. Here, we use simulations carried out with our new Multiphase Darcy-Brinkman-Biot (DBB) framework [27] to fill this knowledge gap and identify nondimensional parameters that govern viscously-stable fluid drainage and fracturing in deformable porous media. We also find that the fracturing dynamics predicted by our continuum-scale framework is consistent with those observed or predicted for granular systems. In other words, in systems with a large length scale separation between pores and fractures, volume-averaged properties are sufficient to capture the onset and propagation of fractures at the continuum scale.

II. MODELING FRAMEWORK
Our investigation is carried out through the use of the Multiphase DBB modeling framework, a new and flexible model used to simulate incompressible two-phase flow through and around deformable porous media [27]. It consists of five volume averaged fluid and solid conservation equations that are coupled by spatially-dependent momentum exchange and capillary force terms. The model is composed of a fluid mass conservation equation, a fluid saturation conservation equation, a fluid momentum conservation equation, a solid mass conservation equation, and a solid momentum conservation equation, In the previous equations, φ f is the fluid volume fraction, φ s is the solid volume fraction, α w is the wetting fluid saturation, α n is the non-wetting fluid saturation, U f is the single-field fluid velocity, U s is the solid velocity, U r is the relative velocity of the two immiscible fluids, p is the single-field fluid pressure, S is the volume averaged fluid viscous stress tensor, σ is the volume averaged solid stress tensor, µk −1 is the drag coefficient (a function of permeability k and single-field fluid viscosity µ), ρ s is the solid density, g is gravity, p c is the capillary pressure, and F c represents additional capillary terms. Here, "singlefield" refers to averaged variables that depend on the properties of both fluids. Lastly, ρ f = α w ρ w + α n ρ n and µ = α w µ w + µ n ρ n are the single-field fluid density and viscosity. The closed form representations for U r , µk −1 , p c , and F c can be found in the Supplemental Materials along with an in-depth description of the model. As indicated in Fig. 1, the system of equations presented above asymptotically approaches the Navier-Stokes multiphase volume-of-fluid [28] equations in solid free regions (where φ f = 1, k is very large, and viscous drag is negligible) and multiphase Biot Theory in porous regions (where φ f < 1, k is small, Re < 1, and drag dominates). This last point can be demonstrated by adding Eqs. 3 and 5 together within a porous domain, which results in the main governing equation used in Biot Theory [27,29,30]: A thorough discussion, derivation, and validation of this model can be found in Carrillo & Bourg 2020 [27] and related publications [18,31]. The two major limitation of the framework highlighted in these previous studies are as follows. First, there needs to be a clear length-scale separation between the averaging volume, the sub-REV heterogeneities, and the overall system [32]. This condition is sustained in most situations involving fractured porous materials, where fracture width is generally significantly larger than the pore width within the porous matrix, with the possible exception of microfractures. Second, closure of the system of equations necessitates the use of parametric models describing the average behaviour of the capillary pressure, permeability, and solid rheology within porous domains. As such, the accuracy of the overall model is inherently impacted by the limitations and assumptions of these parametric models. The complete numerical implementation of the solver, its validations, and the cases shown within this study can be found within the open-source simulation package "hybridBiotInterFoam" [33].

III. NUMERICAL SIMULATIONS
A. Crossover from Imbibition to Fracturing in a Hele-Shaw Cell In addition to the derivation and extensive quantitative validation of Eqs. 1-5, our recent work [27] included a qualitative validation of the ability of the Multiphase DBB model to predict the transition from invasion to fracturing during multiphase flow. Briefly, this validation replicated experiments by Huang et. al. [34] involving the injection of aqueous glycerin into dry sand at incremental flow rates within a 30 by 30 by 2.5 cm Hele-Shaw cell. As shown in Fig. 2, these experiments are inherently multiphysics as fluid flow is governed by Stokes flow in the fracture (aperture ∼cm) and by multiphase Biot Theory in the porous sand (pore width ∼ 100µm). As discussed in our previous work, the similarities between our model and the experimental results are evident: as the viscous forces imposed on the solid increase, so does the system's propensity to exhibit fracturing as the primary flow mechanism (as opposed to imbibition). Minor microstructural differences between our simulations and the experiments reflect the manner in which the implemented continuum-scale rheology model approximates the solid's granular nature. It is clear, however, that both systems are controlled by the balance between viscous forces and solid rheology at the scale of interest [27]. As such, these experiments present an ideal starting point for our investigation.

B. Creation of Fracturing Phase Diagrams
Here, we use the same simulation methodology developed in [27] and illustrated in Figure 2 to identify the general non-dimensional parameters that control the observed transitional behavior between invasion and fracturing in a plastic porous medium. To do so, we systematically vary the solid's porosity (φ f = 0.4 to 0.8), density-normalized plastic yield stress (τ yield = 1.5 to 24 m 2 /s 2 ), capillary entry pressure (p c,0 = 100 to 50, 000 Pa), and permeability (k = 1 × 10 −13 and were taken from Huang et. al. [34] and numerically replicated using equivalent conditions (D, E, F). Black lines represent the advancing saturation front. Additional cases can be found in [27].
5 × 10 −9 m 2 ) as well as the invading fluid's viscosity (µ n = 0.5 to 50 cP) and injection rate (U f = 1 × 10 −4 to 8 × 10 −2 m/s). As in our previous work, the solid's porosity was initialized as a normally-distributed field, the deformable solid was modeled as a Hershel-Bulkley-Quemada plastic [35,36], the porosity-dependence of permeability was modeled through the Kozeny-Carman relation, and relative permeabilities where calculated through the van-Genuchten model [13]. Further details regarding the base numerical implementation of this model can be found in [27], the accompanying code [33], and the Supplementary Materials. The only major differences relative to our previous simulations are that we now include capillary effects and represent viscously-stable drainage as opposed to imbibition (i.e., the injected glycerin is now non-wetting to the porous medium). A representative sample of the more than 400 resulting simulations is presented in the phase diagrams shown in Fig.  3.
Overall, the results make intuitive sense. Figure 3A shows that, ceteris-paribus, less permeable solids are more prone to fracturing. This is due to the fact that, given a constant flow rate, lower permeability solids experience greater drag forces. Our results also show that solids with lower plastic yield stresses fracture more readily, as their solid structure cannot withstand the effects of relatively large viscous or capillary forces. The y-axis behavior of Fig. 3B further shows that systems with higher entry pressures are more likely to fracture, i.e., the capillary stresses are more likely to overwhelm the solid's yield stress, in agreement with grain scale simulations [25]. Finally, Fig. 3B also shows that higher injection rates lead to more fracturing, as these increase viscous drag on the solid structure.

IV. CHARACTERIZATION OF FRACTURING MECHANISMS
The deformation regimes observed in the previous experiments can be delineated by defining two simple nondimensional parameters that quantify the balance between viscous pressure drop, solid softness, and capillary entry pressure.
N cF = p c,0 τ yield ρ s = 2γ r pore τ yield ρ s Here, the viscous fracturing number (N vF ) represents the ratio between the viscous pressure drop and the solid's structural forces. It embodies the question: Does fluid flow generate sufficient friction to induce fracturing? As shown in Fig. 4, the answer is no if N vF < 1 and yes if N vF > 1. This number is the continuum scale analog to the fracturing number presented by Holtzman et. al. [24] for granular solids. It also explains the experimental finding by Zhou et. al. [37] that fracture initiation is only a function of the resulting fluid pressure drop, irrespective of the injection rate or fluid viscosity used to create it. Furthermore, it illustrates why increasing the injection rate and decreasing the permeability have similar effects in Fig. 3.
Complementarily, the capillary fracturing number (N cF ) represents the ratio between the capillary entry pressure and the solid's structural forces; it embodies the question: Does multiphase flow generate sufficient capillary stresses to fracture the solid? Figure 4 shows that when N cF < 1 drainage is the preferential flow mechanism and when N cF > 1 fracturing becomes the dominant phenomenon.
This analysis yields the rudimentary conclusion that fracturing should occur if either of the fracturing numbers is greater than unity, as confirmed by our simulations. However, our simulations further demonstrate the existence of three distinct fracturing regimes (Figs. 3-4). The first regime, referred here as non-invasive fracturing (N vF > 1 and N cF > 1) is characterized by fracturing of the porous solid with minimal fluid invasion, where fractures precede any invasion front. In the second regime, referred to here as the viscous fracturing transition (N vF > 1 and N cF < 1), only the viscous stresses are sufficiently large to fracture the solid. This leads to the formation of relatively wide fractures enveloped and preceded by a non-uniform invasion front. Finally, in the third regime, referred to here as the capillary fracturing transition (N vF < 1 and N cF > 1), only the capillary stresses are sufficiently large to fracture the solid. Given a constant injection rate, this leads to the formation of fractures preceded by an invasion front, as in the viscous fracturing transition regime, but with a more uniform saturation front (due to lower viscous stresses) and less solid compaction (hence narrower fractures). We note that the crossover between each of the four regimes is continuous, meaning that systems with N vF or N cF ∼ 1 can share elements of neighboring regimes.
Although N vF and N cF are fairly intuitive numbers, their impacts on fracture propagation mechanisms are not. For this reason, we also studied the dynamics of fracture nucleation and growth and the evolution of the solid's strain for all three fracturing regimes. As seen in Fig. 5, fracturing in the two transition zones is characterized by the initial formation of non-flow-bearing failure zones (hereafter referred to as cracks), which function as nucleation sites for propagating flow-bearing fractures. These cracks are formed by the simultaneous movement of large contiguous sections of the porous medium in different directions, a process induced by uniform fluid invasion into the porous medium. However, the similarities between both transition zones end here. In the viscous fracturing transition regime, fractures quickly become the dominant deformation mechanism, localizing the majority of the stresses and solid compaction around the advancing fracture tip. Conversely, in the capillary fracturing transition regime, fluid-invasion continues to serve as the main flow mechanism and source of deformation, where fractures and cracks are slowly expanded due to the more evenly-distributed capillarity-induced stresses localized at the advancing invasion front. Finally, noninvasive fracturing follows a different process, where there is little-to-no crack formation and fracture propagation is the main source of deformation and flow. Here, the co- advancing fracture and saturation fronts uniformly compress the solid around and in front of them until this deformation reaches the outer boundary of the simulated system (see the "jet" like-structures at fracture tips in Fig. 5C.) Pressure profiles that further showcase these behaviours can be found in the Supplementary Materials.

V. INFLUENCE OF LOCALIZED AND UNIFORM DEFORMATION
So far we have explored how independently changing k, p c , and τ yield (among others) can affect the fracturing of plastic materials. However, our results also have implications for situations in which these variables are all varied simultaneously, such as during the compaction of soils, sediments, or viscoplastic sedimentary rocks (i.e. mudstones or clay-shales). In such situations, with increasing compaction, k −1 , p c , and τ yield should all increase, although at different rates. As such, we now study the effects of local and uniform deformation on the outlined fracturing regimes.

A. Localized Deformation
The simulations presented above were carried out using the simplifying assumption that p c is invariant with φ f (whereas k and τ yield are not). To evaluate the impact of this simplification on the results shown in Figs. 3-4, we carried out additional simulations for all four regimes with a deformation-dependent capillary entry pressure based on a simplified form of the Leverett J-function where p c,0 = p * c,0 (φ s /φ avg s ) n , p * c,0 is the capillary pressure at φ s = φ avg s , and n > 0 is a sensitivity parameter [38,39]. The results show that non-zero values of n promote the creation of finger-like instabilities and the nucleation of cracks at the fluid invasion front, particularly in the capillary fracturing transition regime. Simulation predictions with different n values are shown in Fig. 6 in the capillary fracturing transition regime and in Supplemental Materials in other regimes.
Despite the additional complexity of the resulting fluid invasion and fracturing patterns, results with n > 0 conform to the overall phase diagram presented in Fig. 4. The results at n = 0 are therefore highlighted in the previous sections due to the greater simplicity of their fluid and solid distribution patterns. 6. Influence of the φ f -dependence of pc on fluid invasion (red and blue) and fracturing patterns (white) in the capillary fracturing transition regime. Here, n represents the sensitivity parameter in the Leverett J-function analogue presented above.

B. Uniform Deformation
Having verified that the applicability of the fracturing numbers holds for systems were k, τ yield , and p c all vary with φ f , we now examine the effects of uniform compaction on said numbers. A direct analysis using the widely-used porosity-parameter relationships implemented above (the Kozeny-Carman relation for k, Leverett J-Function for p c , and Quemada model for τ yield [35,36,38]) yields the following fracturing numberporosity dependence: where D is a rheological parameter based on the solid's fractal dimension (common values range for 1.7-2.9 for different clayey sediments [35]) and φ f,min is the maximum possible degree of compaction. Through these relations, we can see that uniform compaction (or expansion) has a highly non-linear effect on fracturing. Equations 9-10 indicate that whereas N cF tends to consistently decrease with increasing compaction, N vF is considerably more susceptible to changes in φ f and exhibits multiple changes in the sign of its first derivative when D > 2, non-intuitively suggesting that fracturing can be either induced or suppressed through uniform compres-sion. Plots of N vF and N cF as a function of solid fraction are reported in the Supplementary Materials.

VI. CONCLUSIONS
In this article, we used the Multiphase DBB modeling framework to create a phase diagram that identifies two non-dimensional parameters that categorize the crossover between viscously-stable fluid drainage and fracturing as a function of wettability, solid deformability, and hydrodynamics. To the best of our knowledge, our results are the first to relate all three of these properties to characterize multiphase flow in viscoplastic porous media. As expected intuitively, we observe that fracturing occurs if the viscous and/or capillary stresses are sufficient to overcome the solid's structural forces. Thus, when it comes to systems with multiple fluids, it is necessary to consider the effects of surface tension, wettability, and pore size on the fluids' propensity to fracture or invade the permeable solid. Lastly, we find that the two non-dimensional fracturing numbers identified above delineate the existence of three fracturing regimes with distinct fracture propagation mechanisms.

Hershel-Bulkley Plasticity
A Bingham plastic is a material that deforms only once it is under a sufficiently high stress. After this yield stress is reached, it will deform viscously and irreversibly. The Herschel-Bulkley rheological model combines the properties of a Bingham plastic with a power-law viscosity model, such that said plastic can be shear thinning or shear thickening during deformation. In OpenFOAM this model is implemented as follows: where µ e f f s is the effective solid plastic viscosity, which is then modeled through a power law expression: where µ 0 s is the limiting viscosity (set to a large value), τ is the yield stress, µ s is the viscosity of the solid once the yield stress is overcome, n is the flow index (n = 1 for constant viscosity), and η is the shear rate.

Quemada Rheology Model
The Quemada rheology model ; Spearman (2017) is a simple model that accounts for the fact that the average yield stress and effective viscosity of a plastic are functions of the solid fraction. These two quantities are large at high solid fractions and small at low solid fractions, as described by the following relations here, φ max s is the maximum solid fraction possible (perfect incompressible packing), τ 0 is the yield stress at φ s = φ max s /2 , µ 0 is the viscosity of the fluid where the solid would be suspended at low solid fractions (high fluid fractions), and D is a scaling parameter based on the solid's fractal dimension. Figure 2 As shown in our previous work (Carrillo and Bourg, 2020b), the simulations shown in Figure 2 where implemented as follows: Numerical parameters were set to the known properties of aqueous glycerin, air, and sand (ρ gly = 1250 kg/m 3 , µ gly = 5 to 176 cP, ρ air = 1 kg/m 3 , µ air = 0.017 cP, ρ s = 2650 kg/m 3 ), the air-glycerin surface tension γ = 0.063 kg/s 2 , and the sand grain radius r s = 100 µm. To mimic the existence of sub-REV scale heterogeneity, the solid fraction was initialized as a normally-distributed field φ s = 0.64 ± 0.05. To account for