Phase-field modeling of geologic fracture incorporating pressure-dependence and frictional contact

Geologic fractures such as joints and faults are central to many problems in energy geotechnics. Notable examples include hydraulic fracturing, injection-induced earthquakes, and geologic carbon storage. Nevertheless, our current capabilities for simulating the development and evolution of geologic fractures in these problems are still insufficient in terms of efficiency and accuracy. Recently, phase-field modeling has emerged as an efficient numerical method for fracture simulation which does not require any algorithm for tracking the geometry of fracture. However, existing phase-field models of fracture neglected two distinct characteristics of geologic fractures, namely, the pressure-dependence and frictional contact. To overcome these limitations, new phase-field models have been developed and described in this paper. The new phase-field models are demonstrably capable of simulating pressure-dependent, frictional fractures propagating in arbitrary directions, which is a notoriously challenging task.


Introduction
Geologic fractures such as joints and faults are central to a variety of problems in energy geotechnics. Examples range from stimulation and operation of unconventional reservoirs and geothermal systems, to safe sequestration of CO2 and hazardous wastes into deep underground. These problems commonly require engineers the ability to understand, predict, and manage the development of new fractures in geomaterials as well as the activation of existing fractures. Physics-based numerical modeling is an indispensable tool to achieve this ability.
Nevertheless, it remains very challenging to simulate the development and evolution of geologic fractures with sufficient efficiency and accuracy. Much of the difficulty is attributed to the complex geometry of fracture, which demands sophisticated algorithms for its tracking. This is particularly true for simulation of a propagating fracture like that in hydraulic fracturing and fault rupture.
In the computational mechanics community, phasefield modeling has recently emerged as a highly efficient method for simulating propagating fractures without any algorithm for tracking crack geometry (e.g. [1][2][3]). In essence, the phase-field method diffusely approximates a discontinuous interface (e.g. a crack) using a continuous field variable called the phase field, and it models the propagation of a fracture as evolution of the phase-field variable. By doing so, the phase-field method allows one to simulate fracture propagation without explicit tracking their geometry, making it much easier to model complex fracturing processes such as kinking and branching than other types of numerical methods.
Due to this appealing feature, phase-field modeling is increasingly used to simulate fractures in geomaterials especially fluid-driven fractures (e.g. [4][5][6]). In these works, more or less standard phase-field models (e.g. [1]) -developed for other types of solids -were used. Fractures in geologic materials, however, exhibit several characteristics that distinguish themselves from fractures in other solids. The two most critical characteristics are perhaps (i) the pressure-dependence of the fracturing mode, and (ii) frictional contact along a fracture. The first characteristic, which is pressure-induced brittleductile transition, is demonstrated in Fig. 1. The second characteristic is a consequence of that geologic fractures are usually subjected to significant confining pressure in underground. These important characteristics, however, are not incorporated in the existing phase-field models of fracture. Therefore, new phase-field formulations for pressure-dependent, frictional fractures are necessary to realistically simulate important fracturing processes in geologic materials, such as the activation and rupture of shear fractures in unconventional energy [7] as well as induced seismicity [8]. This paper briefly describes new phase-field models recently developed to efficiently simulate pressuredependent and frictional fracture in geologic materials. Herein, the essence of the new models' formulations and verification and validation results are presented. For more detailed descriptions of the models, please refer to journal publications [10,11].
2 Phase-field formulations for pressuredependent, frictional fracture This section outlines the formulations of the phase-field models incorporating pressure-dependence and frictional contact.

Basics of phase-field modeling
In essence, phase-field modeling approximates a discrete (discontinuous) fracture as a diffusive (continuous) zone where the phase-field variable (d) attains a certain value. The value of d ranges from 0 to 1. Usually, d = 1 denotes a fully cracked region and d = 0 denotes an intact region. See Fig. 2 for illustration. In phase-field modeling, the propagation of a fracture is equivalent to evolution of the phase-field variable. The phase-field evolution is governed by a partial differential equation derived from a fracture mechanics theory. In the literature, several different governing equations have been proposed. Among them, the most common form is derived from Griffith's brittle fracture theory and can be written as , W + is the crack driving force, Gc is the critical energy release rate, and L is the length parameter for phase-field approximation. It is noted that although Eq.
(1) will be used in the sequel, other forms of phase-field evolution equations can also be employed.

Incorporating pressure-dependence
The standard phase-field model described above assumes that the fracturing behavior is independent of confining pressure. Yet -as demonstrated in Fig. 1 -confining pressure has a strong impact on the fracturing mode of geologic materials.
To incorporate pressure-dependence into phase-field modeling, it has been proposed to couple the phase-field method with pressure-sensitive plasticity widely used in geomechanics (e.g. Mohr-Coulomb and Drucker-Prager plasticity) [10]. The major challenge for this coupling is how to intertwine phase-field fracture and plasticity in a physically meaningful and numerically efficient way. To address this challenge, it has been proposed to (i) use a phase-field effective stress to consider the damage and fracture effects on plasticity, and (ii) augment the stored energy from plastic deformation (cold work) to the crack driving force. These two approaches can be well justified theoretically and reproduce pressure-induced transition of the fracturing mode as observed from experiments (see Section 3).

Incorporating frictional contact
Frictional contact is ubiquitous in cracks and interfaces in natural and engineered systems including geologic deposits. Therefore, the numerical modeling of interfaces with frictional contact has long been an important topic in computational mechanics. Nevertheless, the majority of phase-field models in the literatures have focused only on opening fractures, and it remained uncertain as to the treatment of frictional contact inside a fracture diffusely approximated by the phase-field variable.
Recently, the authors have developed the first way to incorporate frictional contact into a phase-field fracture [11]. The key idea is to compute the stress tensor inside an interface region, which is diffusely approximated by the phase field, according to the contact condition. The identification of contact condition and the calculation of the interface stress tensor are performed in an interfaceoriented region as depicted in Fig. 3. When the contact condition is found to be a slip condition, the slip-related components of the interface stress tensor are calculated according to a frictional contact law, while all the other stress components are made compatible with the stress tensor in the bulk (intact) region. In this way, one can effectively impose contact behavior and no-penetration constraints without any algorithm. This algorithm-free feature is a crucial advantage of the phase-field method because standard discrete methods for frictional contact problems require special algorithms for imposing contact constraints and oftentimes, they suffer from oscillations in contact normal pressures. So the phase-field method allows for unprecedently effective modeling of frictional fracturing such as the activation and rupture of faults.  Fig. 3. An interface-oriented coordinate system inside a diffusely approximated region. Source: Fei and Choo [11].

Transition from extension to hybrid to shear fractures by confining pressure
The first numerical example is fracturing of laboratoryscale rock specimens under tension. Ramsey and Chester [12] showed that under the same tension test, fracturing mode of Carrara marble shows transition from extension to hybrid to shear fractures by an increase in confining stress. This transition can be evidenced by changes in the fracture angle, as shown in Fig. 4. The fracture angle of extension (Mode I) fracture is close to zero, whereas the angle of shear fracture (Mode II) is 45° minus the half of the friction angle. The new contribution of Ramsey and Chester [12] was, however, the finding of a previously unknown mode -which they called hybrid fracturethat give rise to fracture angles in between the angles of extension and shear fractures. This is the reason why hybrid fracture is highlighted in Fig. 4. Nevertheless, no computational model had been reported to be capable of simulating hybrid fracture, let alone its manifestation during the transition from extension to shear fracture. To investigate whether the coupled phase-field and plasticity model captures the pressure-induced transition of fracturing modes including hybrid fracture, tension tests on dogbone-shape laboratory-scale specimens -as done by Ramsey and Chester [11] were simulated at different confining pressures. The simulation results are shown in Fig. 5. When the confining pressure was 20 MPa, the specimen showed extension fracture without any plastic strain, which indicates purely brittle fracture. Conversely, when the confining pressure was 50 MPa, the same specimen failed by shear fracture developed inside a localized plastic zone. But when the confining pressure was 25 MPa, the developed fracture was not fully inscribed by a localized plastic zone. This indicates that the fracture in this case was a mix of extension and shear fracture, which was called a hybrid fracture in [11]. To confirm this, the fracture angles of the three cases are also plotted in Fig. 5. Comparing these angles with the experimental results in Fig. 4, it can be affirmed that the three fractures are extension, hybrid, and shear fractures, respectively, in the order of increasing pressure.
Through this example, it can be concluded that the proposed coupling of phase-field fracture and pressuresensitive plasticity enables the simulation of all the fracturing mode of geomaterials -extension, hybrid, and shear fractures, and their transition by confining pressure.

Sliding along a frictional interface
The second example is concerned with sliding along a frictional interface. The specific setup of this example is depicted in Fig. 6. This example is a benchmark problem in contact mechanics, which was simulated by a number of previous studies including Simo and Laursen [13] and Oden and Pires [14]. The same problem was revisited by the newly proposed phase-field method incorporating frictional contact [11].   7 compares the deformed shape simulated by the phase-field method of Fei and Choo [12] (color) with that by a standard finite element method of Simo and Laursen (mesh). One can see that both results are nearly undistinguishable, which means that the phase-field method provides highly accurate solutions to frictional contact problems. Furthermore, the contact pressures and shear stresses along the frictional interface are plotted and compared in Fig. 8. It can again be observed that the results of the phase-field method and the standard finite element method are nearly identical. Note, however, that the phase-field method does not require any algorithm for imposing contact constraints as well as for geometry tracking, which dramatically simplifies the problem as compared with discrete methods for the same purpose.
In summary, this example has verified the phase-field method for modeling sliding along a frictional interface. Fig. 7. Comparison of the deformed shapes simulated by the phase-field method (color) and a standard finite element method (mesh). Reproduced from Fei and Choo [11].  [13] and Oden and Pires [14]. Source: Fei and Choo [11].

Propagation of a frictional shear fracture
The third and last example simulates the propagation of a shear fracture from a pre-existing crack. Fig. 9 depicts the setup of this example wherein a rectangular specimen possessing a pre-existing crack is subjected to vertical compression. This setup was borrowed from a numerical example simulated in Borja [15] using three different numerical methods. Here, the problem is revisited by the phase-field method incorporating frictional contact. Fig. 9. Setup of the shear fracture propagation example. The red line denotes the pre-existing discontinuity. Source: Fei and Choo [11]. Fig. 10 presents how the phase-field variable evolves during the compression. Under a small displacement, the phase-field was nearly the same as the initial state due to the frictional resistance of the pre-existing crack. Under further compression, however, the crack propagates and eventually leads to the failure of the overall specimen.
The vertical displacements in the specimen during the fracturing process are presented in Fig. 11, after being normalized by the maximum vertical displacement inside the specimen. The displacement patterns at the vertical displacements of 0.010 m and 0.015 m are very similar to those presented in Borja [15]. This agreement affirms that the new phase-field method incorporate frictional contact accurately even for a propagating fracture.
It is further noted that the geometry of the propagating fracture was neither tracked nor aligned with the element boundaries. This means that the phase-field method can capture a fracture of any geometry very easily. This is a crucial advantage of the phase-field method over other numerical methods in the literature, as geologic fractures often exhibit very complex geometry which is difficult to be captured algorithmically.  Fig. 11. Normalized vertical displacements at the three instances in Fig. 10. Source: Fei and Choo [11].

Closure
This paper has briefly described two recently developed phase-field models of geologic fracture that incorporate pressure-dependence and frictional contact. These two characteristics are ubiquitous in geologic fractures such as joints and faults, playing critical roles in a number of problems in energy geotechnics. Nevertheless, they were notoriously challenging to be simulated especially when the fractures propagate. The demonstrated capabilities of the phase-field models will significantly help us address fractures in geo-technologies for energy production and its environmental impact. Future work includes coupling the phase-field models with fluid flow in porous geomaterials, using advanced numerical methods for coupled poromechanics (e.g. [16][17][18][19]). If successful, the future work will greatly enhance our ability to tackle some of the most pressing problems at present, such as fluid-driven fractures and injectioninduced earthquakes.
Portions of the materials presented in this paper were supported by the Research Grants Council of Hong Kong (Projects 27205918 and 17201419). The second author acknowledges financial support from his Hong Kong Ph.D. Fellowship.