This paper presents an immersed phase field model designed to predict the fracture-induced flow due to brittle fracture in vuggy porous media. Due to the multiscale nature of pores in the vuggy porous material, crack growth may connect previously isolated pores, which leads to flow conduits. This mechanism has important implications for many applications such as disposal of carbon dioxide and radioactive materials and hydraulic fracture and mining. To understand the detailed microporomechanics that causes the fracture-induced flow, we introduce a new phase field fracture framework where the phase field is not only used as an indicator function for damage of the solid skeleton but also used as an indicator of the pore space. By coupling the Stokes equation that governs the fluid transport in the voids, cavities, and cracks and Darcy’s flow in the deformable porous media, our proposed model enables us to capture the fluid–solid interaction of the pore fluid and solid constituents during crack growth. Numerical experiments are conducted to analyze how the presence of cavities affects the accuracy of predictions based on the homogenized effective medium during crack growth.

## I. INTRODUCTION

Geomaterials such as carbonate rocks, sandstone, or limestone often contain geometrical features such as cracks, joints, vugs, or cavities. When the defects are partially or fully saturated with the pore fluid, the geometry of the features may affect effective stiffness, permeability, water retention characteristics, and drained or undrained shear strength of the material (Juanes *et al.*, 2006; Sun *et al.*, 2011b; Kang *et al.*, 2016; Suh *et al.*, 2017; Selvadurai *et al.*, 2017; Wang and Sun, 2017; and Sun and Wong, 2018). Furthermore, brittle fracture in materials that possess geometrical features may lead to the pore fluid in cavities migrating into the flow channels and causing flow conduits that lead to often undesirable outcomes. Modeling geometrical features in porous media is, thus, highly important and, at the same time, a challenging subject for the hydromechanically coupled analysis in geomechanics problems such as hydrocarbon resource recovery or development of enhanced geothermal energy reservoirs (Paterson and Fermigier, 1997; Class *et al.*, 2002; Rutqvist *et al.*, 2007; Grant, 2013; Wagner *et al.*, 2015; Heider and Markert, 2017; and Suh and Sun, 2020).

One possible modeling choice is to consider a fictitious effective medium at a scale where representative elementary volume exists. In this case, the geometrical features of the material are not explicitly modeled, but the influences of these geometrical features are incorporated in the constitutive relations by treating defects as a different pore system that interacts with the matrix pores (Choo and Borja, 2015; Choo *et al.*, 2016; Liu and Abousleiman, 2017; and Wang and Sun, 2018). The upshot of the multi-porosity and multi-permeability models is mainly the simple numerical treatment, since there is no need for complex meshing techniques or embedded strong discontinuities, and the computational efficiency compared to pore-scale models that require extremely large domains in order to reproduce hydromechanical behavior at large scales (Ghaboussi and Barbosa, 1990; Spaid and Phelan, 1997; Blunt *et al.*, 2002; Tang *et al.*, 2005; Arson and Pereira, 2013; Pereira and Arson, 2013; and Suh and Yun, 2018). However, the drawback of this approach is that the homogenized effective medium may not sufficiently represent the microstructural details. This makes the identification of material parameters more complicated since the effective permeability of multiple interacting systems is not isotropic and the constitutive law for fluid mass exchange inherently depends on the microstructure.

Another common alternative to model the interaction between cavities and crack growth is to conduct simulations via a fracture network model (Ozkan *et al.*, 2010; Leung and Zimmerman, 2012; Fu *et al.*, 2013; and Hyman *et al.*, 2015). However, the obvious drawback is that the fracture in those models must be either a straight line (in the two-dimensional case) or a plane (in the three-dimensional case), and hence, the geometrical effect on the porous media cannot be captured precisely.

In this research, we introduce a phase field framework that allows us to enable a unified treatment to simulate the evolving geometry of cracks and cavities. By introducing the phase field as a unified representation of the void space that is not suitable to be treated as an effective medium, we introduce a framework that enables us to analyze how crack propagation in vuggy porous media may affect the flow mechanism differently than that in porous media with pores well distributed in the host matrix. Our result indicates that the interaction between the propagating cracks of the cavities is important for capturing the hydromechanical responses of the porous media, and that existing effective medium approach that characterizes the pore space with a single hydraulic model such as the cubic law and Kozeny–Carman model may not be sufficient to capture the cavity–crack–host matrix interactions.

As for notations and symbols, bold-faced and blackboard bold-faced letters denote tensors (including vectors that are rank-one tensors); the symbol “·” denotes a single contraction of adjacent indices of two tensors (e.g., ** a** ·

**=**

*b**a*

_{i}

*b*

_{i}or

**·**

*c***=**

*d**c*

_{ij}

*d*

_{jk}); the symbol “:” denotes a double contraction of adjacent indices of tensors of rank two or higher (e.g., $C:\epsilon $ =

*C*

_{ijkl}

*ε*

_{kl}); and the symbol “⊗” denotes a juxtaposition of two vectors (e.g.,

**⊗**

*a***=**

*b**a*

_{i}

*b*

_{j}) or two symmetric second order tensors [e.g., (

**⊗**

*α***)**

*β*_{ijkl}=

*α*

_{ij}

*β*

_{kl}]. We also define identity tensors as

**=**

*I**δ*

_{ij}and $I=\delta ik\delta jl$, where

*δ*

_{ij}is the Kronecker delta. As for sign conventions, unless specified, the directions of the tensile stress and dilative pressure are considered positive.

## II. THE MODEL PROBLEM

We consider a fully saturated Biot–Stokes system (Fig. 1) that consists of two regions (intact porous matrix $BD$ and cracks or cavities $BS$) separated by the sharp interface Γ*, where we assume that the solid phase in $BD$ forms a deformable porous matrix while solid particles in $BS$ are in suspension. In this case, both the solid and fluid phases coexist in both regions. By considering our material of interest as a multi-phase continuum, we utilize the effective stress principle for the intact porous matrix where the fluid flow is modeled with Darcy’s law, while the motion of the solid–fluid mixture is modeled by the Stokes equation (Li *et al.*, 2018). Two distinct regions are then coupled by properly imposing three transmissibility conditions at the interface. The model problem with the sharp interface will be later on extended to a diffuse Biot–Stokes model by introducing the phase field in Sec. III.

### A. Continuum representation

Although the Biot–Stokes system only contains two immiscible solid and fluid phases, for mathematical convenience, we idealize the material of interest as a three-phase continuum where each constituent [i.e., solid (*s*), pore fluid (*f*_{D}), and free fluid (*f*_{S})] occupies a fraction of volume at the same material point. By letting *dV* = *dV*_{s} + *dV*_{f} denote the representative elementary volume of the material, we define the volume fractions of the constituents as

where the index *s* refers to the solid phase and *f* indicates the fluid phase. Since the sharp interface separates our system of interest into two regions, the volume fraction of the pore and free fluids can be expressed as

where $H\Gamma *$ is the Heaviside function that satisfies

In addition, by letting *ρ*_{s} and *ρ*_{f} denote the intrinsic mass densities of the solid and fluid, respectively, the partial mass densities of each constituent (*ρ*^{α}, where *α* = *s*, *f*_{D}, *f*_{S}) are given by

where *ρ* is the mass density of the entire system. In this study, we assume that both the solid and fluid phases are incompressible, so that intrinsic mass densities *ρ*_{s} and *ρ*_{f} are considered constants.

### B. Governing equations

#### 1. Conservation laws for an intact porous matrix

For the region where the solid forms an intact porous matrix, we adopt the effective stress principle (Lade and De Boer, 1997; Borja, 2006) so that the external loading imposed on the matrix is assumed to be carried by both the solid skeleton and the pore fluid. In this case, the region $BD$ is governed by the following system of equations (Borja and Alarcón, 1995; White and Borja, 2008; and Sun *et al.*, 2013):

where ** σ**′ is the effective stress,

*B*= 1 −

*K*/

*K*

_{s}is Biot’s coefficient,

*M*is Biot’s modulus, $pfD$ is the pore pressure,

**is the gravitational acceleration, $v$**

*g*_{α}is the intrinsic velocity of the constituent

*α*, and $wfD=\varphi fD(vfD\u2212vs)$ is the Eulerian relative flow vector of the pore fluid (i.e., Darcy’s velocity). Here, we assume that

*B*≈ 1 and 1/

*M*≈ 0 to simplify the formulation. Note that Biot’s coefficients of many sandstone and shale specimens are often less than one, whereas it is more reasonable to assume Biot’s coefficient equal to 1 for granite (e.g., Westerly granite) (Berryman and Wang, 1995; Zimmerman, 2000; and Jaeger

*et al.*, 2009). In either case, the damage of the solid skeleton may reduce the elastic bulk modulus of the solid skeleton. Therefore, both Biot’s coefficient and modulus may evolve according to the solid deformation. This nonlinear effect is not considered in this study but will be considered in the future. We also assume that the behavior of the intact matrix in $BD$ is linear and isotropic elastic, and hence, only two independent elastic moduli are needed to capture the elastic response. The constitutive relation for the solid skeleton can, therefore, be written as follows:

where $\sigma 0\u2032$ indicates the effective stress of the undamaged matrix. The actual and undamaged effective stresses are related by a degradation function, which will later be discussed in Sec. III B. Furthermore, $\epsilon =(\u2207us+\u2207usT)/2$ is the infinitesimal solid strain tensor that depends on the solid displacement *u*_{s}, and parameters *λ* and *μ* are the Lamé constants. For the constitutive equation that describes the laminar pore fluid flow in $BD$, we use the generalized Darcy’s law that linearly relates the relative velocity $wfD$ and pore pressure gradient $\u2207pfD$, i.e.,

where *μ*_{f} is the dynamic viscosity of the pure fluid phase and *k* is the effective permeability of the porous matrix. Additionally, in order to incorporate the effect of deformation of the matrix on the porous medium flow (Mauran *et al.*, 2001; Schutjens *et al.*, 2004), this study adopts the Kozeny–Carman equation to empirically capture the porosity–permeability relation (Chapuis and Aubertin, 2003; Costa, 2006; and Wang and Sun, 2017). Note that the Kozeny–Carman equation is often considered a rough approximation of the porosity–permeability relation. A more precise prediction of permeability may require new geometrical attributes such as tortuosity (Sun *et al.*,2011b; 2011a), formation factor (Worthington, 1993; Sun and Wong, 2018), and percolation threshold (Mavko and Nur, 1997). This extension is out of the scope of this study but will be considered in the future work.

Recall Sec. II A that $\varphi s+\varphi fD=1$ in $BD$. Then, by letting $\varphi \u2254\varphi fD$ denote the porosity of the matrix, the Kozeny–Carman equation reads

where *k*_{0} and *ϕ*_{0} denote the reference permeability and porosity, respectively.

#### 2. Conservation laws for solid–fluid mixture

This study attempts to model suspension flow in $BS$, where mass and linear momentum balances for both solid and fluid phases should be satisfied. We, therefore, write the governing balance equations for $BS$ as

where *σ*^{α} is the Cauchy stress tensor of the *α* constituent, and the relative flow vector of the free fluid can be defined as $wfS=\varphi fS(vfS\u2212vs)$. By assuming that the free fluid resides in $BS$ with low Reynolds number (i.e., *Re* ≪ 1), we adopt a simplified version of the Navier–Stokes model, i.e., the Stokes equation. The Stokes model for the steady-state motion of an incompressible fluid yields the following relationship for the free fluid stress tensor $\sigma fS$:

where $pfS$ is the free fluid pressure and *μ*_{eff} is the effective viscosity of the solid–fluid mixture (Mooney, 1951; Cheng and Law, 2003), i.e.,

where $c\u22541\u2212\varphi fS$ indicates the solid particle concentration and *c*_{max} denotes its upper bound. Again, notice that we introduce only one solid constituent for the entire system since it is convenient for us to later on impose interface conditions and further adopt the phase field fracture model that simulates the evolving interface. This approach may not be suitable for modeling complete suspension flow where $vs=vfS$. However, we assume that solid particles in $BS$ follow the same constitutive relations as the free fluid in order to replicate the suspension flow as close as possible, i.e.,

#### 3. Conservation laws for the sharp interface between intact matrix and solid–fluid mixture

In order to properly model the interaction between the porous matrix ($BD$) and the vugs or cavities ($BS$), complete mass conservation and force equilibrium for the entire system should be satisfied. Since we have two different constituents for the same type of fluid (*f*_{D} and *f*_{S}) while considering only one solid constituent (*s*), coupling two subsystems, thus, requires the enforcement of fluid transmissibility conditions at the sharp interface Γ^{*} that models the coupled Stokes–Darcy flow (Arbogast and Lehr, 2006; Arbogast and Brunson, 2007; Badia *et al.*, 2009; Wu and Mirbod, 2018; and Bergkamp *et al.*, 2020).

The first interface condition is the fluid continuity that ensures mass conservation. Since we assume that the fluid phase is incompressible, the interfacial fluid fluxes for each subsystem ($MfD*$ and $MfS*$) can be expressed as follows:

where $nD*$ and $nS*$ denote the outward-oriented normal vectors from $BD$ and $BS$, respectively. From Eq. (16), mass continuity ($MfD*+MfS*=mfD*+mfS*=0$) yields the following transmissibility condition:

where we take $n*=nS*=\u2212nD*$ for notational convenience (Fig. 1). Here, Eq. (16) implies that the normal component of the fluid velocities ($wfS$ and $wfD$) should be identical in order to guarantee that the exchange of fluid mass between $BS$ and $BD$ is conservative.

The second condition is the force equilibrium at the interface Γ^{*}. From each subsystem, total forces acting on the interface ($\mathcal{F}fD*$ and $\mathcal{F}fS$) may be written as

where $tfD*$ and $tfS*$ indicate the tractions at the interface. The force equilibrium requires $\mathcal{F}fD*+\mathcal{F}fS*=tfD*+tfS*=0$, implying that the normal and shear components should be balanced at the same time. By decomposing the traction vectors as

where $m1*$ and $m2*$ are the interfacial tangent vectors, we get two more transmissibility conditions that describe normal and shear force equilibrium, respectively,

Equation (21) is the Beavers–Joseph–Saffman condition (Beavers and Joseph, 1967; Saffman, 1971; Layton *et al.*, 2002; and Arbogast and Brunson, 2007). This idealized condition relates the slip velocity and the shear stress through the dimensionless slippage coefficient *α*_{SD}, which depends on the microstructural attributes of the interfaces, such as surface roughness and irregular patterns, as well as the flow velocity (Beavers and Joseph, 1967; Terzis *et al.*, 2019; and Guo *et al.*, 2020). The validity and limitations of the Beavers–Joseph–Saffman condition are documented in a number of literature studies such as Auriault (2010), Mikelic and Jäger (2000), and Monchiet *et al.* (2019) and will not be repeated here. Possible extensions of the interface conditions to turbulent and multiphase flows are an active research area that is clearly out of the scope of this study but will be considered in the future.

## III. THE PHASE FIELD BIOT–STOKES MODEL WITH EVOLVING FRACTURES

This section introduces the mathematical model that uses smooth implicit function, i.e., the phase field, to approximate evolving sharp interfaces due to damage. We first review the general procedure that employs an implicit function to approximate sharp interfaces (Sec. III A) shown in Fig. 2. Since the phase field is a smooth representation of the Heaviside function, we derive the corresponding mathematical model that approximates interfacial transmissibility conditions suitable for the diffuse representation of the interface. To capture crack growth according to Griffith’s theory, we adopt the classical variational fracture model to allow for crack growth represented by the evolution of the phase field defined over the spatial domain (Sec. III B). These techniques are then applied to the derivation shown in Sec. III C in which a mathematical model to capture the hydromechanical coupling of pore fluid flows in both the host matrix and evolving interfaces in brittle porous media is discussed. The resultant model does not require locally defined enrichment function or remeshing and can be implemented in a standard finite element or finite element/volume solver.

### A. Diffuse interface approximation

This study employs a diffuse approximation for the sharp interface Γ^{*} by introducing a phase field variable *d* ∈ [0, 1] that varies smoothly from 0 in $BD$ to 1 in $BS$. Specifically, we approximate the interfacial area $A\Gamma *$ as $A\Gamma d*$, which can be expressed in terms of volume integration of the surface density functional $\Gamma d*(d,\u2207d)$ over $B=BD\u222aBS\xaf$ (Miehe *et al.*, 2010; Borden *et al.*, 2012; Suh and Sun, 2019; and Suh *et al.*, 2020),

Here, the size of the diffusive zone [i.e., transition zone where *d* ∈ (0, 1)] is controlled by the regularization length scale parameter *l** such that $A\Gamma d*\u2009\Gamma $-converges to $A\Gamma *$ (Mumford and Shah, 1989), i.e.,

Based on this approach, the phase field *d* and its gradient ∇*d* can be regarded as smooth approximations of the Heaviside function $H\Gamma *$ and the Dirac delta function $\delta \Gamma *$, respectively (Stoter *et al.*, 2017; Suh and Sun, 2020). Therefore, the volume integrals of an arbitrary function $G\u0303$ over $BD$ and $BS$ can, respectively, be approximated as

Similarly, the surface integral of the function $G\u0303$ along the sharp interface Γ^{*} can be approximated as

and we also approximate the normal vector *n*^{*} as

### B. Crack growth approximated by evolving phase field

For completeness, this section reviews the phase field model for brittle fracture. We consider the following surface density functional, which is widely used in modeling brittle or quasi-brittle fracture (Bourdin *et al.*, 2008; Miehe *et al.*, 2010; Borden *et al.*, 2012; Bryant and Sun, 2018; and Suh *et al.*, 2020) that possesses quadratic local dissipation function:

At this point, we highlight that the evolution of the phase field (i.e., propagation of cavities or cracks) is a mechanical process driven by the effective stress ** σ**′. In other words, we assume that the solid skeleton is completely damaged in the liquefied zone $BS$, whereas, in $BD$, the solid skeleton remains undamaged. We, thus, omit the terms that are unrelated to the deformation and fracture in this section. Having critical energy $Gc$ that is required to create new free surfaces, the potential energy density

*ψ*reads

where *ψ*_{bulk}(** ε**,

*d*) is the degrading elastic bulk energy and

*g*(

*d*) = (1 −

*d*)

^{2}is the degradation function that induces energy dissipation. Following Amor

*et al.*(2009), we adopt the additive decomposition scheme that splits the elastic energy

*ψ*

_{e}into compressive ($\psi e\u2212$) and tensile and deviatoric ($\psi e+$) modes, where we only degrade $\psi e+$ in order to avoid crack propagation under compression (Na and Sun, 2018; Bilgen and Weinberg, 2019; and Heider and Sun, 2020), i.e.,

where *K* = *λ* + 2 *μ*/3 is the bulk modulus of the porous matrix and $\u2022\xb1=(\u2022\xb1|\u2022|)/2$ is the Macaulay bracket operator. In this case, the effective stress tensor ** σ**′ can also be decomposed as follows:

where $\sigma 0\u2032\xb1=\u2202\psi e\xb1/\u2202\epsilon $ is the fictitious undamaged effective stress, in which we previously assumed *σ*_{0}′ to be linear elastic [Eq. (7)].

Based on the fundamental lemma of calculus of variations, the damage evolution equation can be obtained by seeking the stationary point where the functional derivative of Eq. (29) with respect to *d* vanishes, i.e.,

where

Here, the superposed prime denotes the derivative with respect to *d* and ∇^{2}(•) = ∇ · ∇(•) is the Laplacian operator. Furthermore, by following the treatment used by Miehe *et al.* (2010), we introduce a history function $H$ that is the pseudo-temporal maximum of the positive energy density ($\psi e+$) in order to ensure the crack irreversibility constraint,

By replacing $\psi e+$ in Eq. (34) with $H$, Eq. (33) finally yields the following phase field equation that governs the evolution of the interface:

Note that we can obtain the diffuse representation of the interface by solving Eq. (36), as shown in Fig. 2.

In this study, we leverage the phase field as an indicator function not only for the location of cracks but also for the location of other defects such as cavities or geometrically complicated voids that do not fit for computational homogenization. This approach may efficiently couple the Stokes flow inside the vugs ($BS$) that interact with the pore fluid in the intact porous matrix ($BD$) while both regions are evolving due to crack growth. A major advantage of this work is that free flow inside the fracture is explicitly replicated, and hence, there is no need to introduce permeability enhancement models (e.g., cubic law) (Witherspoon *et al.*, 1980; Konzuk and Kueper, 2004; and Jin and Arson, 2020). This explicit treatment enables the simulations to remain physical even in the situations (e.g., high Reynolds number, rough surface, aperture variation) where the validity of the cubic law is questioned (Witherspoon *et al.*, 1980; Miehe and Mauthe, 2016; Heider and Markert, 2017; Wang and Sun, 2017; Sun *et al.*, 2017; Choo and Sun, 2018; Chukwudozie *et al.*, 2019; and Wang and Sun, 2019).

Verification and experimental validation of the phase field fracture model for brittle solids have been well documented in the literature. For brevity, similar studies are not provided in this paper. Interested readers may refer to, for instance, studies by Nguyen *et al.* (2016) and Pham *et al.* (2017).

### C. Variational formulation of the phase field Biot–Stokes model

We present an immersed phase field Biot–Stokes model designed to simulate the coupled hydro-mechanical behaviors of flow of vuggy porous media with evolving fractures in the brittle regime. This section omits the gravitational effects for brevity (i.e., ** g** =

**0**).

The model problem with the sharp interface (Sec. II) is the system that possesses two distinct boundaries $\u2202BD$ and $\u2202BS$ that can both be decomposed into Dirichlet ($\u2202BDu$, $\u2202BDp$, $\u2202BSw$, and $\u2202BSp$) and Neumann ($\u2202BDt$, $\u2202BDq$, $\u2202BSt$, and $\u2202BSq$) boundaries satisfying

where the union of $BS$ and $BD$ is $B$ and the boundary domain follows the same treatment. Here, we capture the transition of the constitutive responses of the solid constituent in the intact and liquefied states through a partition of unity argument in the local constitutive responses. As such, we adopt only one solid constituent, and the balance equations of linear momentum in the sub-domains $BD$ and $BS$ [Eqs. (5) and (11)] are combined into one set of equations over the domain $B$. The governing equations for the model problem are summarized as follows:

where the natural and essential boundary conditions are not included for brevity. Following the standard weighted-residual procedure, we multiply Eqs. (39)–(43) with proper weight functions (*η*_{s}, $\xi fD$, $\eta fS$, $\xi fS$, and *ζ*) and integrate over their corresponding domain. The resultant weighted-residual statement reads (Badia *et al.*, 2009; Stoter *et al.*, 2017)

where $t^D$ and $q^D$ are the prescribed traction and flux in the porous matrix, respectively, and $t^S$ is the fluid traction. Then, we directly impose the interfacial transmissibility conditions [Eqs. (44)–(46)] into the field equations [Eqs. (48) and (50)]. Due to the fluid mass continuity (i.e., $mfD*+mfS*=0$), the fourth term on the left-hand side of Eq. (48) becomes

while normal and shear force equilibrium (i.e., $tfD*+tfS*=0$) can be imposed at the third term on the left-hand side of Eq. (49), i.e.,

Finally, we apply Eqs. (24)–(27) in order to convert subdomain integrals ($BD$ and $BS$) into an integral over the entire domain ($B$) and to also transform the interface equations [Eqs. (52) and (53)] into a set of immersed boundary conditions. As a result, we get the weak statements for a phase field Biot–Stokes model, which is to find {$us,pfD,wfS,pfS,d$} such that, for all {$\eta s,\xi fD,\eta fS,\xi fS,\zeta $},

where

Here, as pointed out by Stoter *et al.* (2017), the Γ-convergence ensures that the immersed boundary conditions imposed in Eqs. (56) and (57) are consistent with the interface conditions [Eqs. (44)–(46)] if *l*^{*} → 0, which, in turn, confirms the mass conservation and force equilibrium for the entire system $B$.

## IV. NUMERICAL EXAMPLES

This section highlights the capacities of the immersed phase field model to capture the hydromechanical interactions among the pore fluids in the cavities, the cracks, and the homogenized pore space and the host matrix in two numerical experiments. Our focus is on modeling the problems that involve the mechanically driven pore fluid migration due to deformation and crack growth inside the solid skeleton. The first example simulates the consolidation process of the porous material that contains a semi-circular cavity at the bottom that serves as a pore fluid outlet, while the second problem showcases the fracture-induced Stokes–Darcy flow in a vuggy porous medium.

In order to solve Eqs. (55)–(59) numerically, we adopt the standard finite element method where the solution procedure is based on operator-split (Miehe *et al.*, 2010; Heister *et al.*, 2015; and Suh *et al.*, 2020) that successively updates the field variables. In other words, the phase field *d* is updated first by solving *G*^{d} = 0, while all other field variables are held fixed, and the solver then advances the remaining variables by solving {$Gu,GDp,GSw,GSp$}^{T} = **0**. The implementation of our proposed model including finite element discretization and the solution scheme relies on the finite element package FEniCS (Logg *et al.*, 2012a; 2012b; and Alnæs *et al.*, 2015). It is noted that there exist multiple different strategies to solve the same system of equations. Since the exploration of different solution schemes is out of the scope of this study, we omit the details of the implementation for brevity.

### A. Consolidation of porous matrix with a semi-circular cavity

We first simulate a consolidation problem, which has always been one of the key problems in geotechnical engineering. While the classical consolidation problem considers time-dependent water expulsion from the homogeneous porous material, as illustrated in Fig. 3, this numerical example explores the case where the system includes a cavity at the bottom that serves as a pore fluid outlet. This specific setting is designed to simulate the mechanically driven Stokes–Darcy flow without significant changes in microstructural attributes.

The problem domain is a water-saturated 1 × 2 m^{2} sized rectangular porous matrix ($BD$) that contains a semi-circular cavity ($BS$) whose diameter is 0.2 m. We prescribe a 1 kPa compressive mechanical traction at the top, while zero pressure boundary is imposed at the bottom of the cavity so that the time-dependent dissipation of pore pressure can be observed. The material parameters for this example are chosen as follows: intrinsic mass densities of the solid and fluid: *ρ*_{s} = 2700 kg/m^{3} and *ρ*_{f} = 1000 kg/m^{3}; Young’s modulus and Poisson’s ratio of the solid skeleton: *E* = 100 MPa and *ν* = 0.25; initial permeability and initial porosity of the matrix: *k*_{0} = 1.0 × 10^{−8} m^{2} and *ϕ*_{0} = 0.4; dynamic viscosity of the fluid phase: *μ*_{f} = 1.0 × 10^{−3} Pa s; slippage coefficient *α*_{SD} = 0; and regularization length of the interface *l*^{*} = 0.002 m. Furthermore, we assume that the solid constituent remains intact in $BD$ throughout the simulation, while the free fluid inside the cavity has zero particle concentration (i.e., *c* = 0).

Figure 4 shows the spatial distributions of the prime variables at *t* = 1.0 × 10^{−3} s. Here, we compute fluid pressure and relative fluid velocity of the entire system as $pf=(1\u2212d)pfD+dpfS$ and $wf=(1\u2212d)wfD+dwfS$, respectively, since we have separate degrees of freedom for pore and free fluids residing in regions $BD$ and $BS$. The results imply that the applied mechanical load $t^D$ at *t* = 0 builds up the pore pressure, which, in turn, affects the pore fluid to migrate toward the cavity. Furthermore, the free fluid inside the cavity tends to exhibit higher velocity and lower pressure compared to those of the pore fluid because of different constitutive relations (i.e., Stokes equation and Darcy’s law) in each region. As illustrated in Fig. 5, we also investigate the time-dependent response of the system that clearly describes the consolidation process and, at the same time, highlights the continuous pressure and velocity fields along the *y*-axis (i.e., from the center point of the cavity to the top-central point of the external boundary). At *t* = 0, the entire load is taken by the incompressible pore water, which triggers the fluid flow inside the medium. This fluid flow is accompanied by a dissipation of pore pressure over time and an increase in the compression of the entire system, which is consistent with previous studies on homogeneous materials (Li *et al.*, 2004; White and Borja, 2008; Kim *et al.*, 2009; and Wang and Sun, 2016). In addition, the continuous pressure and velocity profiles imply that our model is capable of imposing mass continuity and force equilibrium at the interface as a set of immersed boundary conditions, which confirms the validity of the model.

### B. Comparison studies on fracture-induced flow in vuggy porous media

In the second set of experiment, we conduct numerical simulations within two different types of domains that possess horizontal edge crack (Fig. 6): one explicitly captures the geometry of the large cavities in the porous media; another one captures the influence of the cavities by increasing the porosity of the homogenized effective medium. While the former approach adopts a more explicit representation of the pore geometry and, hence, may provide more detailed information on the interactions between the vugs and the propagating cracks, the latter approach could be numerically more efficient. Our objective is to demonstrate, quantitatively, the difference in the two approaches such that a fuller picture on the trade-off between computational efficiency, accuracy, and precision of the predictions can be established.

#### 1. Modeling vuggy porous media

As illustrated in Fig. 6(a), we first consider a domain that consists of the porous matrix with explicitly modeled cavities. Our first representation consists of totally nine cavities with different major and minor radii (Table I) that share the same aspect ratio of 2:1 and are tilted by 45°, such that the volume fraction of the cavities *θ*_{cav} is 0.056. Here, we assume that the solid skeleton inside the cavities is completely damaged (i.e., *d* = 1), while the porous matrix initially remains completely undamaged (i.e., *d* = 0). The material properties for this case are chosen as follows: *ρ*_{s} = 2700 kg/m^{3}, *ρ*_{f} = 1000 kg/m^{3}, *E* = 20 GPa, *ν* = 0.2, *k*_{0} = 1.0 × 10^{−12} m^{2}, *μ*_{f} = 1.0 × 10^{−3} Pa s, *α*_{SD} = 0.1, $Gc=20$ J/m^{2}, and *l*^{*} = 0.125 × 10^{−3} m. In addition, the initial particle concentration is chosen as *c*_{0} = 0.6 and its upper bound as *c*_{max} = 0.7, in order to mimic the mudflow inside the cracks or cavities (O’Brien and Julien, 1988; Iverson, 1997).

Index . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . | 9 . |
---|---|---|---|---|---|---|---|---|---|

Major radius r_{a} (mm) | 0.400 | 0.600 | 0.500 | 0.500 | 0.820 | 0.800 | 0.580 | 0.700 | 0.650 |

Minor radius r_{b} (mm) | 0.200 | 0.300 | 0.250 | 0.250 | 0.410 | 0.400 | 0.290 | 0.350 | 0.325 |

Index . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . | 9 . |
---|---|---|---|---|---|---|---|---|---|

Major radius r_{a} (mm) | 0.400 | 0.600 | 0.500 | 0.500 | 0.820 | 0.800 | 0.580 | 0.700 | 0.650 |

Minor radius r_{b} (mm) | 0.200 | 0.300 | 0.250 | 0.250 | 0.410 | 0.400 | 0.290 | 0.350 | 0.325 |

In contrast, our second domain in Fig. 6(b) is a homogenized representation of Fig. 6(a), where all the cavities are considered a part of matrix pores. In this case, the porosity of the homogenized medium is determined as *ϕ*_{hom} = (1 − *θ*_{cav})*ϕ*_{0} + *θ*_{cav} = 0.433. It is noted the correct homogenized effective properties often depend on the geometry of vugs or inclusions, which can be determined from computed tomographic images or directly obtained from the experiment (Sun *et al.* 2011a; 2011b; Kim *et al.*, 2016; and Lee *et al.*, 2017). Since the micro-structural attributes are not always available, this study adopts an alternative approach where the effective material properties are determined by using the equivalent inclusion method (Hashin, 1960; Zimmerman, 1991; and Ramakrishnan and Arunachalam, 1993). Following Ramakrishnan and Arunachalam (1993) and by assuming that the matrix shares the same material properties as those chosen in Fig. 6(a), the effective bulk modulus (*K*_{hom}) and shear modulus (*μ*_{hom}) for the homogenized representation [Fig. 6(b)] are determined as follows:

so that the effective Young’s modulus *E*_{hom} = 16.50 GPa and Poisson’s ratio *ν*_{hom} = 0.206. Assuming that the inclusion permeability is much higher than the matrix permeability, we approximate the effective permeability (*k*_{hom}) by following Markov *et al.* (2010), which is obtained based on Maxwell’s formula, i.e.,

In addition, since all the cavities in Fig. 6(a) are completely isolated, we adopt the following effective critical energy $Gc,hom$ proposed by Jelitto and Schneider (2018) for the homogenized representation, which depends on the volume fraction of the cavities:

#### 2. Mechanically driven fracture-induced flow

As illustrated in Fig. 6, we conduct two different types of simulations within each domain: the tension tests with the prescribed vertical displacement rate of 0.01 × 10^{−3} m/s and the shear tests with the horizontal displacement rate of 0.01 × 10^{−3} m/s. In both tension and shear tests, the displacements are prescribed at the upper boundary, whereas the bottom part of the domain is held fixed. We also impose hydraulically insulated boundary conditions for the left and right boundaries, while we permit water intake from the upper and lower boundaries by imposing $p^fD=0$.

Figure 7 illustrates the evolution of the phase field for both tension and shear tests in a computational domain where the cavities are explicitly modeled, compared with the crack trajectories obtained from the homogenized domain. The domain without cavities exhibits the crack patterns that are similar to the results of previous studies on homogeneous solids (Miehe *et al.*, 2010; Borden *et al.*, 2012; Bryant and Sun, 2018; and Suh *et al.*, 2020), while the domain with explicitly modeled cavities exhibits distinct crack patterns. More importantly, Figs. 8 and 9 reveal that neglecting the interaction between the cavity and the crack in the homogenized model may lead to over-simplified global responses that lack the distinctive characteristics of the cavity–crack coalescence.

During the numerical experiments, the porous matrix initially undergoes linear elastic deformation until the crack nucleation takes place. At this point, since tensile loading directly influences the volume change of the material, both specimens under tensile load exhibit higher fluid influx at the top, compared to those measured from the shear tests. After the first peaks shown in Figs. 8(a) and 8(b), cracks start to initiate from the tips of the pre-existing flaw since they experience higher stress concentration compared to the matrix–cavity interface.

In both tensile and shear experiments performed on the vuggy specimen, the crack nucleation increases the surface influx rate at the permeable boundaries as the pore fluid starts to leak from the intact matrix to the damaged regions regardless of the spatial homogenization. The two numerical specimens, nevertheless, begin to behave differently when the cracks propagate toward the adjacent vugs and coalesce with each other in both tension and shear tests in the vuggy specimen (Fig. 7) (Qinami *et al.*, 2019; Suh and Sun, 2020). These changes in surface influx cannot be replicated in the homogenized porous specimen as homogenization takes away the possibility of simulating the coalescence between the cavity and the crack.

After the coalescence of the cavity and the crack in the vuggy specimen, the reaction force in both cases increases again with lower influx rates until it reaches the second peak (i.e., where crack nucleation takes place at the matrix–cavity interface), and the crack eventually reaches both ends of the specimens. This result implies that the existence of vugs or cavities has a profound impact on the material behavior that cannot be easily replicated in the homogenized effective medium. Consequently, either a more effective macroscopic theory or a suitable multiscale technique is needed to incorporate the cavity–crack interaction into the predictions.

Figure 10 illustrates the pressure (*p*_{f}) and *x*-directional velocity ($wfS$) fields from a domain with explicitly modeled cavities under tensile and shear loadings, at *t* = 0.476 s and *t* = 1.012 s, respectively, where cracks start to propagate from the cavities. Here, the superimposed arrows in Fig. 10 indicate the direction of the velocity vector $wf=(1\u2212d)wfD+dwfS$. In both cases, the leakage of pore fluid takes place toward the interconnected cracks and cavities at the middle, while the free fluid in $BS$ tends to migrate toward the center, the region where a large crack opening displacement occurs. However, it is worthy to note that the fluid flow occurs from the region that has negative pressure to the damaged zone where $pfS\u22480$. Unlike previous studies that use the cubic law to predict the hydraulic responses of the flow conduit (Mauthe and Miehe, 2017; Heider and Markert, 2017; Wang and Sun, 2017; and Chukwudozie *et al.*, 2019), the pore pressure distribution inside the void space is governed by the Stokes equation directly. This set of numerical experiments again highlight that our proposed model is capable of simulating the fracture–cavity interaction with the evolving interface, which may not be easily captured either by using hydraulic phase field fracture models or by adopting the classical Biot–Stokes model with a sharp interface.

To assess the computational efficiency of the proposed model, we record the CPU time for both simulations. A laptop with an Intel Core i9-9880H Processor CPU with 16 GB memory at 2667 MHz (DDR4) is used to run both simulations on a single core. Both simulations are solved by the same Scalable Nonlinear Equation Solver (SNES) available in FEniCS. In the case where vuggy pores are explicitly modeled, the time taken to assemble the system of equations is 1.13 s and the average time taken to advance one time step with (on average) five Newton–Raphson iterations is 35.69 s. Meanwhile, in the homogenized case, it takes 1.17 s to assemble the system of equations and 33.34 s to advance one time step also with (on average) five Newton–Raphson iterations. In this example, simulations with the explicitly captured vuggy pores require about 7% more CPU time to run the same simulation.

The future work may consider flow with a higher Reynolds number suitable for the Navier–Stokes equation in the fluid domain. Such an extension is, nevertheless, out of the scope of the current study.

## V. CONCLUSION

This article presents a new immersed phase field model that captures the hydro-mechanical coupling mechanisms in vuggy porous media where brittle cracks filled with water may coalescence with pores that trigger both redistribution of flow and macroscopic softening that cannot be captured without the Stokes–Darcy flow. By generalizing the phase field as an indicator of defects, we introduce a simple and unified treatment to handle the evolving geometries due to crack growth and the resultant changes in constitutive responses without the need of re-meshing or introduction of enrichment functions. By directly simulating the flow inside the cracks, we bypass the need of introducing a phenomenological permeability enhancement model to replicate the flow conduit. This explicit approach can be advantageous over the embedded discontinuity approach when there are a substantial crack opening and a flow near the locations with void–crack interactions where a homogenized pore pressure jump would not be sufficient to capture the pattern of the pore pressure field in the defects. The future work may include the extension of the proposed model to three-dimensional cases, as well as extending the Stokes–Darcy flow model to the generalized Navier–Stokes–Darcy flow for injection and other problems with higher Reynolds numbers.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

The suggestions of the two anonymous reviewers are gratefully acknowledged. The authors are thankful for the fruitful discussion with Professor Robert Zimmerman while the second author visited London in 2017. The first author was supported by the Earth Materials and Processes program from the US Army Research Office under Grant No. W911NF-18-2-0306. The second author was supported by the NSF CAREER grant from the Mechanics of Materials and Structures program at National Science Foundation under Grant No. CMMI-1846875 and the Dynamic Materials and Interactions program from the Air Force Office of Scientific Research under Grant No. FA9550-17-1-0169. These supports are gratefully acknowledged. The views and conclusions contained in this article are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the sponsors, including the Army Research Laboratory or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for governmental purposes notwithstanding any copyright notation herein.

## REFERENCES

_{2}storage

_{2}using coupled fluid flow and geomechanical fault-slip analysis

_{2}-acidized water flow through stressed carbonate rocks