An immersed phase field fracture model in fluid-infiltrating porous media with evolving Beavers-Joseph-Saffman condition

. This study presents a phase field model for brittle fracture in fluid-infiltrating vuggy porous media. While the state-of-the-art in hydraulic phase field fracture considers Darcian fracture flow with enhanced permeability along the crack, in this study, the phase field not only acts as a damage variable that provides diffuse representation of cracks or cavities, but also acts as an indicator function that separates the domain into two regions where fluid flows are governed by Stokes and Darcy equations, respectively. Since the phase field and its gradient can be respectively regarded as smooth approximations of the Heaviside function and Dirac delta function, our new approach is capable of imposing interfacial transmissibility conditions without explicit interface parametrizations. In addition, the interaction between solid and fluid constituents is modeled by adopting the concept of mixture theory, where the fluid velocities in Stokes and Darcy regions are considered as relative measures compared to the solid motion. This model is particularly attractive for coupled flow analysis in geological materials with complex microstructures undergoing brittle fracture often encountered in energy geotechnics problems, since it completely eliminates the needs to generate specific enrichment function, integration scheme, or meshing algorithm tailored for complex geological features.


Introduction
Defects, such as cracks, joints, vuggy pores and cavities and impurities are important for the hydro-mechanical coupling encountered in porous media in energy geotechnics problems like enhanced oil recovery or the development of enhanced geothermal energy reservoirs.When the defects are partially or fully filled with pore fluid, the size and the geometric features of the defects may both impose significant changes on both the effective stiffness and permeability of the materials as well as the Biot's coefficient [30] and the undrained and drained shear strength [25,32].For instance, carbonate rocks and limestone often contain pores of profoundly different sizes [5,7,17,27,28].
One possible modeling choice is to not explicitly model each defect and imperfects but instead incorporate the influences of these defects in the constitutive laws of an imaginary effective medium at a scale where a representative elementary volume exists.In this case, defects may be treated as a different pore system that may interact with the prime pore space through fluid mass exchanges, as shown in the multi-porosity and multipermeability models in the literature [6,11,12,33].The upshot of this model is mainly the relatively simple numerical treatment.Due to the absence of defects in the homogenized effective medium, there is no need for complex meshing techniques.The absence of embedded strong discontinuities or sharp gradient of phase field or level set also make it more feasible to employ standard finite element or finite element/finite volume solver, provided that the numerical stability is ensured [26,29,34].However, the drawback of these approaches is that the homogenized effective medium may not be sufficiently representing the microstructure.Furthermore, the identification of material parameters can become more complicated as 1) the effective permeability of the multiple interacting pore systems are generally not anisotropic and not co-axial with each other and 2) the constitutive law for the fluid mass exchanges are inherently dependent on the geometric attributes such as surface areas and the mechanisms may change over time as cracks grow, coalescence and heal and interact with pores exciting in the porous media.
Explicitly capturing the geometrical evolution of these defects may resolve these issues.However, this approach has been an ongoing challenge for the poromechanics community in the past decades [5,8,31].From the modeling prospective, it is essential to employ a proper representation for the defects either explicitly (e.g.insertion of cohesive zone element, enriching basis function or assumed strain methods) or implicitly (e.g. the phase field, level set or eigenfracture approaches).While both approaches have achieved a level of success in the past decades for solid mechanics problems, the modeling effort of the hydraulic responses of defects are often limited to cubic law models that relates hydraulic aperture to hydraulic conductivity [18,21].While the cubic law can be applied to onedimensional case where local aperture does not fluctuate, the actual flow pattern inside the crack is not necessarily one dimensional.In fact, the macroscopic flow pattern inside the crack is also dominated by the crack geometry and hence the cubic law or simply enhancing the permeability is not universally feasible and certainly not appropriate in highly heterogeneous domain.
In this work, our contribution is to introduce the coupling of the Stokes-Darcy flow that enable a unified treatment for flow in the defects.By leveraging the phase field not only as an indicator function for the location of cracks but also other defects such as cavities and large or geometrically complicated voids that does not fit for computational homogenization, we introduce a phase field framework that may efficiently couple the Stokes flow in the defect regions that interact with the pore fluid infiltrating in the intact porous matrix while the Stokes and Darcy regions are evolving due to the crack growth.By explicitly modeling the fluid-filled defects, we bypass the need to introducing additional phenomenological models for the mass exchanges and enable modelers to capture how different geometries of these defects, meso-scale voids and cracks affect the macroscopic outcome.Numerical examples are provided to showcase the potential applications of this proposed models.

Diffuse representation of Stokes-Darcy system
In this study, we model coupled free and porous medium flow by considering two subdomains (i.e., cracks or vugs and porous matrix , where fluid flow is governed by Stokes and Darcy's equations, respectively) separated by a sharp interface (Fig. 1). is the porous matrix while represents a subdomain for vuggy regions.
Similar to [22], we attempt to employ a diffuse approximation for the sharp interface between two regions via implicit function.By adopting a phase field approach which is widely used in modeling fracture [3,14,23,24], we approximate the interfacial area as , which can be expressed in terms of volume integration of surface density over : (1) where is the phase field which varies from 0 in Darcy region ( ) to 1 in Stokes region ( ).This study considers the density functional originally used to introduce elliptic regularization of the Mumford-Shah functional [16] which possesses quadratic local dissipation function, i.e., Here, is the regularization length that controls the size of diffuse interface.By solving a differential equation which can be recovered by seeking the stationary point where the functional derivative of the density functional in Eq. ( 2) vanishes, we can obtain diffuse representation of the Stokes-Darcy system (e.g., Fig. 2).Since -convergence requires that the sharp interface is recovered by reducing the regularization length to zero [1,4], the phase field and its gradient can be regarded smooth approximations of the Heaviside function and Dirac delta function , respectively [22].Therefore, the volume integral of an arbitrary function over Stokes and Darcy regions can be approximated as follows: (3) Likewise, the surface integral of the function over the interface is approximated as, (4) where the magnitude of the gradient of the phase field approximates a Dirac delta distribution at the sharp interface .We can also obtain the normal ( ) and tangent ( ) vectors (Fig. 1) from the phase field representation of the interface.The normal vector can be computed as: E3S Web of Conferences 205, 03009 (2020) (5) In the general three-dimensional case, the tangent vectors can then be obtained as, (6)

Governing equations for hydraulic phase field fracture
In this study, the interaction between solid and fluid constituents is modeled by adopting the mixture theory, where the fluid velocities and in both regions become relative measures compared to the solid motion.For simplicity, we assume that the vuggy porous material is composed of incompressible solid and fluid constituents under quasi-static conditions.We also assume that effective stress principle is only valid in while solid constituent inside is fully liquefied.In this case, we may write local forms of the governing equations as: (7) (8) (9) where is the solid displacement; and are fluid pressures in Stokes and Darcy regions, respectively; indicates fluid viscosity; is the effective stress; and is the critical fracture energy.Here, notice that Eqs.(7.1) and (7.2) describes the solid and fluid motions in and , respectively.By adopting the effective stress principle, we incorporate the damage in the solid matrix from the phase field evolution equation [Eq.( 9)] as follows: (10) where is the degradation function, and are the tensile and compressive parts of the fictitious undamaged effective stress [13,14].In order to enforce crack irreversibility, we introduce a history function in Eq. ( 9), which is a pseudo-temporal maximum of the tensile part of the effective strain energy density ( ) stored in solid skeleton, i.e., (11) In addition, the undamaged solid matrix in is considered to be linear elastic, i.e., (12) and fluid flux in is described by Darcy's law: (13) where is the effective permeability of the porous matrix that depends on Kozeny-Carman equation.
The Stokes flow and the Darcy flow can properly be coupled along the interface by enforcing three transmissibility constraint therein.The transmissibility conditions govern the interaction between free and porous-medium flows.The first interface condition is the continuity equation that ensures the fluid mass conservation across : (14) Here, notice that ; where is the absolute fluid velocity in region and is the porosity.The second condition is the force equilibrium equation that balances the normal stresses across the interface: (15) where indicates the fluid traction vector acting on .The third condition is the well-known Beavers-Joseph-Saffman law [2,10,15,20] which correlates the slip velocity with the shear stress along the interface : (16) where i = {1, 2}, and is the dimensionless Beavers-Joseph slippage coefficient which can be determined experimentally [2].
Starting from the strong form, we first follow the standard procedure to recover the variational form of Eqs. ( 7)-( 9) that contains transmissibility conditions [Eqs.( 14)-( 16)] while employing the low-order finite elements for , , , and .By using Eqs.( 3) and (4), we then reformulate the sharp interface variational form into the variational form of the hydromechanically coupled diffuse Stokes-Darcy problem.The solution procedure is based on the operator-split scheme to successively update the field variables.In this operator-split setting, the phase field is updated first while all other field variables are held fixed, and then the solver holds the phase field and advances the variables , , and .For brevity, we omit the detailed variational formulation.Interested readers are referred to [22], where they adopt similar approach to replace integrals over sharply defined volumes and surfaces by diffuse integrals formulated in terms of phase field.

Numerical example
This section presents a numerical example to showcase the applicability of the proposed phase field model for fluid-infiltrating vuggy porous media.For simplicity, we limit our attention to two-dimensional case.As illustrated in Fig. 3, our problem domain consists of a square-shaped porous matrix that possesses an initial notch and two cavities at the middle.The material parameters for this example are chosen as follows: Young's modulus GPa, Poisson's ratio , initial permeability m 2 , fluid viscosity , initial porosity , Beavers-Joseph slippage coefficient , critical energy release rate J/m 2 , and phase field regularization length mm.Within the same computational domain, two different types of simulations are conducted: the pure tensile test with prescribed vertical displacement mm/s; and the pure shear test with prescribed horizontal displacement mm/s.In both cases, the displacements are prescribed along the entire top boundary, while the bottom part of the domain is fixed.Furthermore, zero pore pressure conditions ( ) are applied at the top and the bottom, while no-flux boundary conditions are applied at the left and the right boundaries.Fig. 4 illustrates the predicted crack trajectories for both tension and shear tests, while Fig. 5 shows their global force-displacement responses.Since notch tip experiences higher stress concentration compared to matrix-cavity interface, cracks tend to initiate from the tips of the initial notch, and then coalescence with cavities.Notice that, in pure shear test, inclined crack propagation occurs at the beginning, but eventually coalesces with the cavities due to an increase of the singularity [19].After coalescence, pure tensile loading exhibits the horizontal crack pattern due to the symmetry of the boundary value problem, whereas pure shear test yields inclined crack pattern.Fig. 6 illustrates the resultant pressure and velocity fields for tension and shear tests at t = 0.91 sec and t = 2.03 sec, respectively.Compared to the permeability enhancement models [9,13], in our new approach, fluid pressure inside the fracture and cavity remains low while high negative pressure builds up in surrounding Darcy region due to the transmissibility conditions along the interface.In addition, compared to the undamaged porous matrix, fluid velocity inside the fracture and cavity exhibits higher values, since it can be interpreted as a region that possesses extremely high permeability.The results highlight that our new approach is capable of modeling fracture-cavity interaction in fluid-infiltrating porous media by explicitly modeling the interface conditions.

Conclusion
This study presented a phase field method for brittle fracture in fluid-infiltrating vuggy porous media.To the best of the authors' knowledge, this is the first ever mathematics model that employs the phase field fracture framework combined with coupled Stokes-Darcy flow.By adopting the mixture theory, we started with sharp interface formulation with particular emphasis on the interfacial transmissibility conditions, and then transferred it into a diffusive problem by assuming the existence of a phase field representation of the interface.We believe that our approach is a significant step towards hydraulic fracturing simulation capability, since this model allows one to explicitly model friction at the crack surface and to explore fracture-cavity interaction without any numerical treatment.

Fig. 1 .
Fig. 1.Schematic representation of computational domain of interest.isthe porous matrix while represents a subdomain for vuggy regions.

Fig. 2 .
Fig. 2. Diffuse representation of the interface where exemplary 1D domain consists of Stokes region , sandwiched between Darcy regions.

Fig. 5 .
Fig. 5.The force-displacement curves obtained from tension and shear tests.