## 1 Introduction

Subsurface systems are replete with fractures of various kinds. Examples include joints and faults, which are tensile and shear fractures, respectively, from the viewpoint of fracture mechanics (Pollard and Fletcher, 2005; Schultz, 2019). Slip surfaces can also be viewed as shear discontinuities, of which propagation can be described by fracture mechanics (Palmer and Rice, 1973; Puzrin and Germanovich, 2005). The deformation and growth of these fractures govern the stability and performance of a wide range of geotechnical systems from slopes and earth retaining systems to tunnels and caverns. Also, modern subsurface energy technologies such as unconventional resource recovery and enhanced geothermal systems purposely generate fractures in deep underground deposits. Unlike other types of fractures, the geologic fractures in these problems often involve frictional contact and show significant pressure dependence. Also, fractures in geomaterials are quasi-brittle instead of purely brittle, and their critical fracture energies in tensile and shear usually differ by an order of magnitude. Computational modeling is essential to addressing such complex behavior of geologic fractures.

A variety of approaches have been developed for computational modeling of geologic fractures. Following the classical approaches in computational fracture mechanics, earlier works treated geologic fractures as sharp discontinuities across which displacement and/or strain fields are discontinuous. These approaches honor the sharpness of crack geometry, being consistent with fracture mechanics principles derived from sharp cracks. However, it is challenging to represent sharp crack geometry in numerical methods. For example, the standard finite element method can only model cracks aligned with element boundaries. To overcome this limitation, embedded finite elements (e.g., Regueiro and Borja, 1999; Foster et al., 2007) and extended/generalized finite elements (e.g., Liu and Borja, 2009; Sanborn and Prévost, 2011) have been proposed for modeling geologic fractures that pass through inside elements. However, these non-standard finite element methods require not only additional shape functions for enrichment but also sophisticated algorithms (e.g., the level-set method) to track crack surfaces.

Over the past decade, the phase-field method has become a popular approach for computational modeling of fractures (e.g., Bourdin et al., 2008; Miehe et al., 2010b; Borden et al., 2012). The phase-field method approximates a sharp crack surface as a diffuse interface and describes its evolution by a partial differential equation formulated from fracture mechanics concepts. In this way, the phase-field method can simulate complex cracking processes such as kinking, branching, and coalescence without any surface-tracking algorithms, being solvable by the standard finite elements widely used in practice. While a phase-field model demands a significant computational cost, it has been increasingly affordable and widely applied to many kinds of fracture problems in various materials. Nevertheless, the standard phase-field model of fracture and its variants do not incorporate several distinct features of geologic fractures: frictional contact, pressure-dependence, quasi-brittleness, mode-mixity, and roughness, among others. Therefore, although these phase-field models have been employed in many geomechanical applications (e.g., Lee et al., 2016; Zhang et al., 2017; Choo and Sun, 2018a; b; Ha et al., 2018; Zhou et al., 2018), there is a need for new phase-field models tailored to geologic fractures.

Recently, the author and his coworkers have developed a class of novel phase-field models for geologic fractures, which incorporate frictional contact (Fei and Choo, 2020a), quasi-brittle shear fracture with friction and pressure effects (Fei and Choo, 2020b), mixed-mode fracture (Fei and Choo, 2021), roughness effects (Fei et al., 2022), and inertia effects with rate- and state-dependent friction Fei et al. (2023). This report reviews these phase-field models with particular attention to the double-phase-field model (Fei and Choo, 2021) which allows one to simulate the combination of cohesive tensile fracture and frictional shear fracture without any algorithms for surface tracking and contact constraints. The double-phase-field model has not only been well-validated against various experimental data (Fei et al., 2021) but has also been applied to gain insights into geologic fractures that are challenging to investigate experimentally (Choo et al., 2023; Sun et al., 2024). The remainder of this paper introduces the formulation, validation, and applications of the phase-field model.

## 2 Phase-field model

### 2.1 Geometry approximation

The departure point of the phase-field modeling of fracture is to diffusely approximate the sharp geometry of a crack (e.g., a line crack in a 2D domain, and a surface crack in a 3D domain) by a phase field,

Mixed-mode fractures—the combination of tensile and shear fractures—are common in geomaterials. Standard phase-field models approximate both tensile and shear fractures with the same phase-field variable. This approach works well when the characteristics of tensile and shear fractures are similar, but it is undesirable when the two types of fractures exhibit contrasting characteristics. The latter is the case for geologic fractures: Shear fractures in geomaterials have much higher fracture energies than tensile fractures and involve frictional effects. Therefore, Fei and Choo (2021) have proposed a double-phase-field model that approximates tensile and shear fractures separately by two distinct phase fields. Figure 1 illustrates the double-phase-field approximation of mixed-mode fractures.

Figure 1

**Figure 1**. Phase-field approximation of mixed-mode (tensile and shear) fractures. Source: Fei and Choo (2021).

### 2.2 Formulation

Phase-field modeling describes the cracking process as the evolution of the phase-field variable over time. The governing equation for the phase-field evolution can be derived either from a variational principle or microforce theory. The latter approach allows one to easily incorporate complex features such as quasi-brittleness, contact dependence, and friction effects. From microforce theory, the governing equations for tensile and shear fractures, respectively, are derived as.

$\nabla \cdot \left(\frac{\partial \psi \left(\mathit{\u03f5},{d}_{I},\nabla {d}_{I},{d}_{II},\nabla {d}_{II}\right)}{\partial \nabla {d}_{I}}\right)-\frac{\partial \psi \left(\mathit{\u03f5},{d}_{I},\nabla {d}_{I},{d}_{II},\nabla {d}_{II}\right)}{\partial {d}_{I}}=0,\phantom{\rule{0ex}{0ex}}(1)$

$\nabla \cdot \left(\frac{\partial \psi \left(\mathit{\u03f5},{d}_{I},\nabla {d}_{I},{d}_{II},\nabla {d}_{II}\right)}{\partial \nabla {d}_{II}}\right)-\frac{\partial \psi \left(\mathit{\u03f5},{d}_{I},\nabla {d}_{I},{d}_{II},\nabla {d}_{II}\right)}{\partial {d}_{II}}=0,\phantom{\rule{0ex}{0ex}}(2)$

where

$\psi ={\psi}^{\text{e}}+{\psi}^{\text{f}}+{\psi}^{\text{d}},\phantom{\rule{0ex}{0ex}}(3)$

where

$-{g}_{I}^{\prime}\left({d}_{I}\right){\mathcal{H}}_{I}-\frac{{\mathcal{G}}_{I}}{\pi L}\left(2{L}^{2}\nabla \cdot \nabla {d}_{I}-2+2{d}_{I}\right)=0,\phantom{\rule{0ex}{0ex}}(4)$

$-{g}_{II}^{\prime}\left({d}_{II}\right){\mathcal{H}}_{II}-\frac{{\mathcal{G}}_{II}}{\pi L}\left(2{L}^{2}\nabla \cdot \nabla {d}_{II}-2+2{d}_{II}\right)=0.\phantom{\rule{0ex}{0ex}}(5)$

Here,

$\theta =\mathrm{arg}{\mathrm{max}}_{\theta}\left[\mathcal{F}\left(\theta \right)\right]{|}_{\mathit{\u03f5}},\phantom{\rule{0ex}{0ex}}\phantom{\rule{0ex}{0ex}}\text{where}\phantom{\rule{0ex}{0ex}}\phantom{\rule{0ex}{0ex}}\mathcal{F}\left(\theta \right)\u2254\frac{{\mathcal{H}}_{I}\left(\mathit{\u03f5},\theta \right)}{{\mathcal{G}}_{I}}+\frac{{\mathcal{H}}_{II}\left(\mathit{\u03f5},\theta \right)}{{\mathcal{G}}_{II}},\phantom{\rule{0ex}{0ex}}(6)$

where

The degradation functions,

To make the stress-strain response of a phase-field model independent of the length parameter, one must employ a particular degradation function derived to provide a length-insensitive stress-strain response. The double-phase-field model employs two such degradation functions, one derived for cohesive tensile fracture (Wu, 2017) and the other derived for frictional shear fracture (Fei and Choo, 2020b). In this way, it is possible to directly use the tensile and shear strength properties measured from experiments, without concerning the chosen value of

Last but not least, it should be emphasized that the double-phase-field formulation can be solved well with the standard finite element method. (See Fei and Choo (2021) for the full finite element formulation.) Specifically, one can use a three-field finite element discretization in which the displacement field and the two phase fields are the unknown variables, with a material (integration) point update algorithm to calculate contact-dependent stresses. The displacement fields and phase fields can be solved in a staggered manner for robustness (Miehe et al., 2010a). In doing so, no algorithm is necessary for tracking the crack geometry or imposing the contact constraints. This ease of implementation with the standard finite element method is a distinct advantage of the phase-field method.

### 2.3 Material parameters

The material parameters of the double-phase-field model are summarized as follows.

It is noted that when fracture test data are available, the tensile and shear fracture energies can also be calibrated to match the fracture data.

### 2.4 Validation

The double-phase-field model has been validated with several experimental data on mixed-mode fracture in rock specimens with preexisting flaws under compression. Fei and Choo (2021) and Fei et al. (2021) present qualitative and quantitative validations with data from uniaxial compression tests on rock specimens with various flaw configurations. These papers have shown that the double-phase-field model can simulate complex mixed-mode cracking patterns and stress-strain responses similar to the experimental data. Yet these validations were restricted to 2D (planar) fracture under uniaxial compression. More recently, in Sun et al. (2024), the double-phase-field model has been validated with fractures in two types of more general conditions, namely, 1) 2D fracture under true triaxial compression, and 2) 3D fracture under uniaxial compression. It is noted that although 3D fractures under true triaxial compression are the most general condition, they cannot be characterized experimentally through the existing techniques. In the following, the validation of 2D fractures under true triaxial compression is reported briefly.

To validate the double-phase-field method for cracking behavior under true triaxial compression, it was applied to simulate a true triaxial compression test on a fine sandstone specimen with two coplanar, fully-penetrating 2D flaws (Zhou et al., 2021). The width, length, and ligament length of each flaw were 1mm wide, 16mm, and 8mm, respectively, and the flaw inclination angle was

Figure 2 compares the simulation and experimental results of mixed-mode fracture from two coplanar flaws under true triaxial compression. One can see that the crack type, crack geometry, and crack coalescence pattern of the simulation and experimental results are quite similar. From the outer tips of the flaws, primary and secondary tensile cracks emerged, and then oblique shear cracks developed from both the inner and outer tips of the flaws. Also, from the shear crack fronts, anti-wing cracks grew toward the direction of the major principal stress. Eventually, the two flaws coalesced into a shear crack, and coplanar shear cracks nucleated from the outer tips of the flaws and grew along the flaw planes. Apart from this qualitative agreement in terms of the cracking pattern, the simulation and experimental results show good quantitative agreement in terms of the stress-strain curves. It is noted that the same conclusion has been reached from other validation studies (Fei and Choo, 2021; Fei et al., 2022; Sun et al., 2024).

Figure 2

**Figure 2**. Validation of the double-phase-field model with experimental data on cracking from two coplanar flaws under true triaxial compression (Zhou et al., 2021). **(A)** Cracking patterns from the simulation and experiments. **(B)** 3D mixed-mode processes in the simulation (Points A–E correspond to those in the stress-strain curve below). **(C)** Difference between the major principal stress *versus* strain from the simulation and experiments. (Note that the experimental stress-strain curve is redrawn without the initial nonlinear portion related to the closure of preexisting defects.) Source: Sun et al. (2024).

## 3 Applications of phase-field model

This section introduces two examples of how the phase-field model has been applied to investigate mixed-mode fracture in rocks which are extremely challenging to characterize by the existing experimental methods. The two examples are: 1) size effect on mixed-mode fracture in rocks with preexisting flaws (Choo et al., 2023), and 2) intermediate principal stress effect on the 3D cracking behavior of rocks under true triaxial compression (Sun et al., 2024).

### 3.1 Size effect on mixed-mode fracture in rocks

Laboratory specimens with preexisting flaws, such as the specimen simulated in the validation example, have been conventionally used as small-scale analogs of rock masses, and their failure behavior under compression has been extensively studied through experimental and numerical methods. However, no study has been concerned with the energetic size effect—determined by the relative size between the fracture process zone (FPZ) and the structure size—on the failure behavior of rock specimens with preexisting flaws. In quasi-brittle materials like rocks, the size of FPZ is significant in relatively small structures like laboratory specimens, and it gives rise to rather ductile failure behavior that deviates from the description of linear elastic fracture mechanics (LEFM). However, in large structures like field-scale rock masses, the size of FPZ is negligible compared with the structural size; in this case, the structure fails like a purely brittle material described by LEFM. This energetic size effect was first revealed in Bažant (1984) in the context of tensile failure in quasi-brittle materials and has been investigated for various types of quasi-brittle failures. Nevertheless, the energetic size effect on mixed-mode fracture in rocks under compression has not been investigated systematically. One primary reason may be that it is highly challenging to examine the size effect on mixed-mode fracture in flawed rocks under compression with the existing experimental methods. For example, it is extremely difficult to prepare field-scale (meter-scale) rock specimens with preexisting flaws and characterize the mixed-mode fracture under compression.

The double-phase-field model can be an ideal tool to investigate the energetic size effect on mixed-mode fracture in flawed rocks. The reasons are: 1) it can simulate complex crack patterns without geometry tracking algorithms, 2) it can distinguish between tensile and shear cracks naturally, and 3) it combines phase-field models of cohesive tensile and shear fractures in quasi-brittle materials. As such, in Choo et al. (2023), the double-phase-field model is leveraged to perform the first systematic investigation of mixed-mode fracture in rocks under compression. It should be noted that unlike tensile fracture where size affects the strength but not crack morphology, mixed-mode fracture in rocks under compression may be subjected to two types of size effects in terms of crack morphology as well as strength. Indeed, experimental studies have found different cracking patterns in rocks and rock-like materials with the same specimen and flaw geometry (e.g., Lee and Jeon, 2011), and the difference may be attributed to the different brittleness/ductility of the materials. For this reason, it is intriguing to explore whether structure size also affects the mixed-mode cracking pattern in flawed rocks under compression.

To investigate the size effect on the mixed-mode cracking behavior of rock masses under compression, the double-phase-field model is used to simulate a series of uniaxial compression tests on single- and double-flawed rock specimens of seven different sizes with geometrical similarity. The sizes of the seven specimens were determined such that they range from millimeter-scale specimens to meter-scale specimens. Specifically, setting a 76.2mm wide and 152.4mm tall laboratory specimen as a reference, seven values of the scaling factor

Figure 3 presents the simulated mixed-mode cracking patterns in double-flawed specimens of seven sizes (

Figure 3

**Figure 3**. Mixed-mode cracking patterns in double-flawed specimens with different scaling factors:

### 3.2 Intermediate principal stress effect on the 3D cracking behavior of rocks

So far, mixed-mode fracture in rocks has mostly been studied under uniaxial compression. However, uniaxial compression is far from *in-situ* stress conditions of rocks. Underground rocks are usually under true triaxial stress conditions in which all three (major, intermediate, and minor) principal stresses are compressive and distinct. Even at surfaces such as excavation boundaries, the stress condition is biaxial in which only the minor principal stress is zero. Therefore, the intermediate principal stress, *in-situ* stress states, and hence its effect should be considered properly.

Besides the stress condition, the 3D geometrical feature of the preexisting flaw(s) is a critical factor of cracking behavior under compression. Most existing studies have studied cracks emanating from 2D (planar, penetrating) preexisting flaws. However, a few studies have shown that when the preexisting flaw is 3D (internally embedded), new types of cracking patterns such as petal cracks and crack wrapping emerge (e.g., Dyskin et al., 2003; Yin et al., 2014; Lu et al., 2015). Yet such 3D cracking behavior has only been studied under uniaxial and biaxial compression regimes. Therefore, the 3D cracking behavior of flawed rocks under a wide range of true triaxial stress conditions remains elusive.

Indeed, it is virtually impossible to experimentally characterize 3D cracking processes in rocks under true triaxial stress conditions. This is because while a high-speed imaging system is necessary to identify a mixed-mode fracture process faithfully, it cannot be applied to a rock specimen under a true triaxial cell. Alternatively, high-fidelity numerical simulations based on sound physical principles can be employed to investigate 3D cracking processes under true triaxial compression. To this end, the double-phase-field method can again serve as an ideal tool because of its ability to simulate complex 3D tensile and shear fractures individually.

In Sun et al. (2024), a series of numerical true triaxial compression tests were performed on single-flawed and double-flawed cubic specimens. The flaws were internally embedded such that fully 3D cracking patterns could develop. The test protocol was designed to investigate how the orientation of

Figure 4

**Figure 4**. 3D cracking patterns in a double-flawed specimen (inclination angle: 45**(A)** The intermediate principal stress **(B)** The minor principal stress **(C)** The major principal stress

Remarkably, based on the observations made from the numerical true triaxial compression tests, three mechanisms of the cracking behavior of 3D flawed rocks under true triaxial conditions were proposed. First, the normal stress on the flaw controls the tensile fracture. Second, the Coulomb stress on the flaw controls the shear fracture. Third, the Coulomb-to-normal stress ratio affects the mixed-mode cracking pattern which controls the peak stress. These mechanisms would not have been uncovered without the use of the double-phase-field model.

## 4 Closure

This paper has introduced a novel phase-field approach for computational modeling of geologic fractures and its application to investigating the mechanics of complex rock fractures. Particular emphasis has been placed on the double-phase-field model developed for mixed-mode fracture in rocks and similar quasi-brittle materials. The double-phase-field model has two standout features: 1) it can simulate complex fractures without sophisticated algorithms for crack geometry tracking and contact constraints, and 2) it can naturally distinguish between tensile and shear fractures. The combination of these two features makes the model an ideal tool for studying the mechanics of complex mixed-mode fractures that cannot be well characterized by the existing experimental methods. The model is however not without drawbacks. The most critical drawback is that phase-field modeling entails a significant computational cost because an extremely fine discretization is necessary around the crack. Another drawback is that it is less straightforward to extract discrete quantities (e.g., aperture) from phase-field approximated fractures. Work is underway to overcome these drawbacks through the combined use of novel formulations, efficient algorithms, machine learning methods, and modern computing platforms.

## Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

## Author contributions

JC: Writing–original draft, Conceptualization, Funding acquisition, Methodology.

## Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (No. RS-2023-00209799).

## Acknowledgments

The author wishes to express his deep thanks to his former students Fan Fei and Yuan Sun for their collaboration on the works reported in this paper.

## Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

Bažant, Z. P. (1984). Size effect in blunt fracture: concrete, rock, metal. *J. Eng. Mech.* 110, 518–535. doi:10.1061/(asce)0733-9399(1984)110:4(518)

Bobet, A., and Einstein, H. (1998). Fracture coalescence in rock-type materials under uniaxial and biaxial compression. *Int. J. Rock Mech. Min. Sci.* 35, 863–888. doi:10.1016/s0148-9062(98)00005-9

Borden, M. J., Verhoosel, C. V., Scott, M. A., Hughes, T. J., and Landis, C. M. (2012). A phase-field description of dynamic brittle fracture. *Comput. Methods Appl. Mech. Eng.* 217, 77–95. doi:10.1016/j.cma.2012.01.008

Bourdin, B., Francfort, G. A., and Marigo, J.-J. (2008). The variational approach to fracture. *J. Elast.* 91, 5–148. doi:10.1007/s10659-007-9107-3

Choo, J., Sohail, A., Fei, F., and Wong, T.-f. (2021). Shear fracture energies of stiff clays and shales. *Acta Geotech.* 16, 2291–2299. doi:10.1007/s11440-021-01145-5

Choo, J., and Sun, W. (2018a). Coupled phase-field and plasticity modeling of geological materials: from brittle fracture to ductile flow. *Comput. Methods Appl. Mech. Eng.* 330, 1–32. doi:10.1016/j.cma.2017.10.009

Choo, J., and Sun, W. (2018b). Cracking and damage from crystallization in pores: coupled chemo-hydro-mechanics and phase-field modeling. *Comput. Methods Appl. Mech. Eng.* 335, 347–379. doi:10.1016/j.cma.2018.01.044

Choo, J., Sun, Y., and Fei, F. (2023). Size effects on the strength and cracking behavior of flawed rocks under uniaxial compression: from laboratory scale to field scale. *Acta Geotech.* 18, 3451–3468. doi:10.1007/s11440-023-01806-7

Dyskin, A., Sahouryeh, E., Jewell, R., Joer, H., and Ustinov, K. (2003). Influence of shape and locations of initial 3-D cracks on their growth in uniaxial compression. *Eng. Fract. Mech.* 70, 2115–2136. doi:10.1016/s0013-7944(02)00240-0

Fei, F., and Choo, J. (2020a). A phase-field method for modeling cracks with frictional contact. *Int. J. Numer. Methods Eng.* 121, 740–762. doi:10.1002/nme.6242

Fei, F., and Choo, J. (2020b). A phase-field model of frictional shear fracture in geologic materials. *Comput. Methods Appl. Mech. Eng.* 369, 113265. doi:10.1016/j.cma.2020.113265

Fei, F., and Choo, J. (2021). Double-phase-field formulation for mixed-mode fracture in rocks. *Comput. Methods Appl. Mech. Eng.* 376, 113655. doi:10.1016/j.cma.2020.113655

Fei, F., Choo, J., Liu, C., and White, J. A. (2022). Phase-field modeling of rock fractures with roughness. *Int. J. Numer. Anal. Methods Geomechanics* 46, 841–868. doi:10.1002/nag.3317

Fei, F., Mia, M. S., Elbanna, A. E., and Choo, J. (2023). A phase-field model for quasi-dynamic nucleation, growth, and propagation of rate-and-state faults. *Int. J. Numer. Anal. Methods Geomechanics* 47, 187–211. doi:10.1002/nag.3465

Fei, F., Sun, Y., Wong, L. N. Y., and Choo, J. (2021). “Phase-field modeling of mixed-mode fracture in rocks with discontinuities: from laboratory scale to field scale,” in *55th US rock mechanics/geomechanics symposium. ARMA–21–1223*.

Foster, C. D., Borja, R. I., and Regueiro, R. A. (2007). Embedded strong discontinuity finite elements for fractured geomaterials with variable friction. *Int. J. Numer. Methods Eng.* 72, 549–581. doi:10.1002/nme.2020

Ha, S. J., Choo, J., and Yun, T. S. (2018). Liquid CO_{2} fracturing: effect of fluid permeation on the breakdown pressure and cracking behavior. *Rock Mech. Rock Eng.* 51, 3407–3420. doi:10.1007/s00603-018-1542-x

Ingraffea, A. R., and Heuze, F. E. (1980). Finite element models for rock fracture mechanics. *Int. J. Numer. Anal. Methods Geomechanics* 4, 25–43. doi:10.1002/nag.1610040103

Lee, H., and Jeon, S. (2011). An experimental and numerical study of fracture coalescence in pre-cracked specimens under uniaxial compression. *Int. J. Solids Struct.* 48, 979–999. doi:10.1016/j.ijsolstr.2010.12.001

Lee, S., Wheeler, M. F., and Wick, T. (2016). Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model. *Comput. Methods Appl. Mech. Eng.* 305, 111–132. doi:10.1016/j.cma.2016.02.037

Liu, F., and Borja, R. I. (2009). An extended finite element framework for slow-rate frictional faulting with bulk plasticity and variable friction. *Int. J. Numer. Anal. Methods Geomechanics* 33, 1535–1560. doi:10.1002/nag.777

Lu, Y., Wang, L., and Elsworth, D. (2015). Uniaxial strength and failure in sandstone containing a pre-existing 3-D surface flaw. *Int. J. Fract.* 194, 59–79. doi:10.1007/s10704-015-0032-3

Miehe, C., Hofacker, M., and Welschinger, F. (2010a). A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator splits. *Comput. Methods Appl. Mech. Eng.* 199, 2765–2778. doi:10.1016/j.cma.2010.04.011

Miehe, C., Welschinger, F., and Hofacker, M. (2010b). Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations. *Int. J. Numer. methods Eng.* 83, 1273–1311. doi:10.1002/nme.2861

Palmer, A. C., and Rice, J. R. (1973). The growth of slip surfaces in the progressive failure of over-consolidated clay. *Proc. R. Soc. Lond. A. Math. Phys. Sci.* 332, 527–548. doi:10.1098/rspa.1973.0040

Pollard, D. D., and Fletcher, R. C. (2005). *Fundamentals of structural geology*. Cambridge University Press.

Puzrin, A. M., and Germanovich, L. (2005). The growth of shear bands in the catastrophic failure of soils. *Proc. R. Soc. A Math. Phys. Eng. Sci.* 461, 1199–1228. doi:10.1098/rspa.2004.1378

Regueiro, R. A., and Borja, R. I. (1999). A finite element model of localized deformation in frictional materials taking a strong discontinuity approach. *Finite Elem. Analysis Des.* 33, 283–315. doi:10.1016/s0168-874x(99)00050-5

Sanborn, S. E., and Prévost, J. H. (2011). Frictional slip plane growth by localization detection and the extended finite element method (xfem). *Int. J. Numer. Anal. Methods Geomechanics* 35, 1278–1298. doi:10.1002/nag.958

Shen, B., and Stephansson, O. (1994). Modification of the G-criterion for crack propagation subjected to compression. *Eng. Fract. Mech.* 47, 177–189. doi:10.1016/0013-7944(94)90219-4

Sun, Y., Fei, F., Wong, L. N. Y., and Choo, J. (2024). Intermediate principal stress effects on the 3D cracking behavior of flawed rocks under true triaxial compression. *Rock Mech. Rock Eng.* 57, 4607–4634. doi:10.1007/s00603-024-03777-x

Wong, L. N. Y., and Einstein, H. H. (2009). Crack coalescence in molded gypsum and carrara marble: part 2—microscopic observations and interpretation. *Rock Mech. Rock Eng.* 42, 513–545. doi:10.1007/s00603-008-0003-3

Wong, N. Y. (2008). *Crack coalescence in molded gypsum and Carrara marble*. Ph.D. thesis, Massachusetts Institute of Technology.

Wu, J.-Y. (2017). A unified phase-field theory for the mechanics of damage and quasi-brittle failure. *J. Mech. Phys. Solids* 103, 72–99. doi:10.1016/j.jmps.2017.03.015

Yin, P., Wong, R., and Chau, K. T. (2014). Coalescence of two parallel pre-existing surface cracks in granite. *Int. J. Rock Mech. Min. Sci.* 68, 66–84. doi:10.1016/j.ijrmms.2014.02.011

Zhang, X., Sloan, S. W., Vignes, C., and Sheng, D. (2017). A modification of the phase-field model for mixed mode crack propagation in rock-like materials. *Comput. Methods Appl. Mech. Eng.* 322, 123–136. doi:10.1016/j.cma.2017.04.028

Zhou, S., Zhuang, X., Zhu, H., and Rabczuk, T. (2018). Phase field modelling of crack propagation, branching and coalescence in rocks. *Theor. Appl. Fract. Mech.* 96, 174–192. doi:10.1016/j.tafmec.2018.04.011