A three-field phase-field model for mixed-mode fracture in rock based on experimental determination of the mode II fracture toughness

In this contribution, a novel framework for simulating mixed-mode failure in rock is presented. Based on a hybrid phase-field model for mixed-mode fracture, separate phase-field variables are introduced for tensile (mode I) and shear (mode II) fracture. The resulting three-field problem features separate length scale parameters for mode I and mode II cracks. In contrast to the classic two-field mixed-mode approaches, it can thus account for different tensile and shear strength of rock. The two phase-field equations are implicitly coupled through the degradation of the material in the elastic equation, and the three fields are solved using a staggered iteration scheme. For its validation, the three-field model is calibrated for two types of rock, Solnhofen Limestone and Pfraundorfer Dolostone. To this end, double-edge notched Brazilian disk (DNBD) tests are performed to determine the mode II fracture toughness. The numerical results demonstrate that the proposed phase-field model is able to reproduce the different crack patterns observed in the DNBD tests. A final example of a uniaxial compression test on a rare drill core demonstrates that the proposed model is able to capture complex, 3D mixed-mode crack patterns when calibrated with the correct mode I and mode II fracture toughness.


Introduction
Accurate prediction of fracture in rock and rock-like materials is vital for a number of engineering applications, ranging from building processes to deep geothermal applications.In recent years, numerical methods have vastly complemented geo-mechanical testing and can provide important insights into crack initiation and propagation processes.As an alternative to discrete approaches such as XFEM [35], and cohesive zone models [38], the phase-field approach to fracture [17,23] has gained more and more popularity.Due to its elegant way of representing the crack using a smooth and continuous scalar-variable and its formulation as a minimization problem, the phase-field approach facilitates the solution of complex fracture scenarios.In contrast to the aforementioned discrete approaches, crack propagation follows directly from the solution of a partial differential equation without the need for complex remeshing procedures or ad-hoc criteria for crack initiation.Consequently, a wide range of phase-field approaches have been proposed including models for ductile fracture [2,50], heterogenenous [27] or anisotropic material [48] and specific materials such as fiberreinforced concrete [1,54] and poro-elastic media [56].For the simulation of fracture in rock it is important to account for the difference in mode I (tensile) and mode II (shear) fracture resistance, as the fracture toughness for mode II fracture is usually higher.The first and most intuitive phase-field approach to capture the mixed-mode behavior in rock was presented by Zhang [55].Here, different critical energy release rates for mode I and mode II fracture are introduced and the crack driving force is split into two separate parts which correspond to the different crack modes [55].Bryant and Sun [18] propose a modification of the mixed-mode approach with consistent kinematics based on the determination of the local crack propagation direction.The approach by Fan [21] extends the splitting method for masonry-like material [24] to account for mixed-mode behavior by introducing a split into mode I and mode II components based on the local crack direction similar to [18].In contrast to the mixed-mode model proposed in [55], the latter two methods do not suffer from an overestimation of the driving force under pure mode I loading [21].However, they require the solution of a maximization problem to determine the local crack driving direction.A further drawback of the above mentioned methods is related to the length scale parameter for the phase-field regularization.As the material strength depends on the choice of the length-scale parameter, Tanné [47] suggests to regard the length-scale as a material property and calibrate it with the material's tensile strength.This, however, can only describe the nucleation of mode I cracks.To overcome this problem, [22] propose a length insensitive multi-phase-field formulation for the simulation of mixed-mode fracture in quasi-brittle materials.Based on the ideas presented in [10], the approach uses two different phase-fields, one for cohesive tensile fracture and one for frictional shear fracture.In the present contribution, we address the length scale problematic by proposing a three-field phase-field model that uses two different length scales, one for mode I and one for mode II failure.The model can be calibrated using the respective tensile and shear strength of the material.The flexible setting is easy to implement and allows for different splits between mode I and mode II components.That way, the three-field model can be tailored for specific applications and to the available computational resources.Proper calibration of the three-field model for different rocks requires their specific and unique mechanical properties.Decisive parameters for the rock's plastic behavior, its elastic properties and tensile strength, can easily be determined by standardized tests, such as the uniaxial compression test [20,36] and indirect methods like the Brazilian disc test [5,8,46].Mixed-mode phase-field models also require the critical energy release rates to properly assess mode I and II fracture initiation and propagation.A variety of different tests have been proposed for the determination of the mode I fracture toughness [31,51] and numerous data has been collected for different types of rock.However, data on the mode II fracture toughness is limited.To obtain the true mode II fracture toughness not only the loading of the test specimen has to be in mode II, but the crack initiation has to be driven by shear.Only a limited number of tests for the determination of the mode II fracture toughness have been proposed, including the punch through shear test [6] and the shear box test [41].Recently, a Double-edge Notched Brazilian Disk (DNBD) test was suggested by Bahrami [7], which features a simple experimental setup, enables the determination of the true mode II fracture toughness, and readily allows for the observation of fracture patterns using high-speed cameras.In the present contribution, we present a full workflow based on the experimental determination of the mode II fracture toughness using DNBD tests.To the authors' knowledge, this is the first mixedmode model which is calibrated and successfully applied to reproduce both 2-and 3-dimensional, mixed-mode fracture scenarios.The paper is structured as follows.In Section 2, the three-field model and its discretization with the Finite Cell Method is introduced.The DNBD experiments including the computation of the mode II fracture toughness are presented in Section 3, followed by the numerical results in Section 4. Here, the proposed three-field is validated based on the DNBD tests and a complex application example of a uniaxial compression test after ISRM SM 1979, with determined uniaxial compressive strength and its stress-strain curve, is presented.

A three-field phase-field formulation
In this chapter, the theoretical and numerical background of the phase-field formulation are introduced.The proposed phase-field model is a three-field problem based on the ideas presented in [10] and [55] with two separate phase-field variables associated to mode I and mode II fracture, respectively.

Mathematical formulation
Let Ω ⊂ R d , d = 2, 3 be an open bounded domain which is cut by a set of discrete cracks, as shown in Figure 1 A point in Ω is denoted by x and u(x), ε(x) and σ(x) ∈ R d are the displacement, strain and stress fields, respectively.An isotropic and linear elastic material with small deformations and quasi-static conditions is assumed.In this case, the strain tensor given as ε = 1 2 (∇u + ∇ T u) and the elastic strain density as Ψ(ε) = 1 2 λ tr 2 (ε) + µ tr(ε 2 ), where λ and µ are the Lamé constants.

Background
The phase-field approach to fracture is based on the variational formulation by Francfort [23] and its subsequent regularization by Bourdin [15,16].Here, the discrete crack is approximated using a scalar variable s, the so-called phase-field, which smears the crack over a regularization width l 0 .The phase-field parameter attains a value of zero on the crack and is one if the material is undamaged.Crack propagation is considered as a minimization problem of the associated functional Here, G c is the critical fracture energy, g(s) is the degradation function which models the loss of stiffness in the material due to damage, w(s) is the energy dissipation function, and c w a scaling parameter.Formulation (1) suffers from interpenetration of crack surfaces and non-physical crack patterns in compression.Thus, commonly an additive split of the elastic strain energy density Ψ(ε) = Ψ + (ε) + Ψ − (ε) is used.Based on the spectral split proposed in [34], a hybrid formulation was introduced in [3].Here, the split is only accounted for in the phase-field equation, which results in a linear elastic problem.This reduces the computational effort while providing comparable results [3].
Following variational theory the Euler-Lagrange equations of the functional Eq. ( 1) can be derived, which yields the following coupled system of equations for the hybrid formulation Here, the AT-2 model with w(s) = 1 − s 2 and c w = 1/2 ([47]) is used.A quadratic degradation function g(s) = (1 − η) s 2 + η is chosen, where η is a small numerical parameter which prevents full degradation of the material.The coupled system (2) is subject to the boundary conditions Fig. 1: Sharp crack topology (left) and phase-field crack surface (right) adapted from [25].The discrete cracks Γ c I and Γ c II are represented by two distinct phase-field variables s I and s II ranging from 0 to 1.For visualization purposes, these are summarized in a variable s ∈ [−1, 1], which takes a value of 0 in the undamaged region, 1 on a tensile crack and −1 on a shear crack.
The history variable H, as introduced by Miehe [34], ensures irreversibility of the phase-field and facilitates the use of a staggered solution scheme.It is defined as where • denotes the Macaulay brackets with • + = x − | x | and ε + is the positive part of the strain tensor resulting from a spectral split.In the coupled system (2), the ratio H/G c drives the evolution of the phase-field.Although the critical fracture energies for mode I and mode II fracture can vary considerably, this is not accounted for in the formulation above.To overcome this limitation Zhang [55] proposed a mixed-mode modification of (2), where two critical fracture energies G c I and G c II are introduced.The driving force is replaced with a weighted average of mode I and mode II driving forces weighted by their respective critical fracture energies.The phase-field equation (2a) is modified according to where the mode I and mode II driving forces H I and H II are defined as The mixed-mode formulation by Zhang [55] uses a single length-scale parameter for both tensile and shear fracture, and thus, is not able to account for different tensile and shear strength of the material.
To overcome this limitation we extend the formulation above to a three-field problem, which will be introduced in the next section.

The three-field phase-field model
Based on the formulation by Zhang [55], we propose a three-field phase-field formulation which introduces different scalar variables for tensile and shear failure.The scalar variables s I , s II ∈ [0, 1] represent mode I and mode II fracture, respectively, with s I = 0 on a tensile crack and s II = 0 on a shear crack as depicted in Figure 1.The driving force is split up following 8.Here it is assumed that the tensile phase-field s I is driven by H I , while the shear field s II is driven by H II .Adopting a phase-field evolution according to (2b) for each of the damage variables the three-field problem is obtained as Here, a degradation function g(s I , s II ) is defined as which accounts for the damage of mode I and mode II cracks.Different length scale parameters l 0,I and l 0,II for shear and tensile fracture, respectively, are introduced.In the case of pure mode I or pure mode II failure, the formulation falls back to the original mixed-mode formulation (7).It should be noted, that the split of the elastic strain energy density in tensile and shear components as proposed by [55] suffers from an overestimation of the force response under pure mode I loading [55].
Alternatively, different approaches based on the split by Amor [4] or directional dependent splits based on local crack coordinates as proposed by Strobl [45] or Steinke [43] can be integrated in the proposed three-field formulation.An overview of existing splitting methods can be found in [21].In the following, we

Discretization
The numerical framework is based on the approach presented in [28,37], which combines the phasefield approach with an embedded domain technique, the finite cell method [39], and multi-level hp-adaptive refinement [53].

The Finite Cell Method
The finite cell method (FCM) is based on an implicit representation of the geometry.Instead of generating a boundary conforming mesh, the actual geometry is recovered during integration with the help of an indicator function.As depicted in Figure 2, the physical domain Ω phy is embedded into a larger domain of simple shape Ω ∪ = Ω phy ∪ Ω f ict which can easily be meshed.To account for the actual geometry, an indicator function α(x) is defined which takes a value close to zero in the surrounding, so-called fictitious domain and is equal to one in the physical region: Here, α is a small numerical parameter greater than but unequal to zero to ensure stability [39].
The weak form is multiplied by α(x) eliminating contributions from the fictitious domain.Advanced integration schemes such as quad-and octree-subdivision approaches are needed for a sufficiently accurate integration of the cells cut by the domain boundary [19].For further elaboration on the FCM and its combination with multi-level hp-adaptive refinement the reader is referred to [19] and [53].

Weak Form
Let the trial spaces for the displacement solution S u , the mode I phase-field solution S s I and the mode II phase-field solution S s II be defined as Fig. 2: Embedding concept of the Finite Cell Method following [39].
where H 1 refers to the Sobolev space of degree one.Furthermore, let the spaces for the test functions be defined as The weak formulation of the coupled three-field problem for the FCM states: Find u ∈ S u , s I ∈ S s I and s II ∈ S s II such that (σ, ∇w) , the penalty method is used to apply Dirichlet boundary conditions for the elastic problem in a weak sense, where β is the penalty parameter.

Solution of the Coupled Three-Field Problem
The system (19) is discretized in a finite element setting using integrated Legendre polynomials as basis functions for the finite test and trial spaces as explained in [37].In each displacement step of the quasi-static simulation, a staggered solution scheme is used to solve the discretized equations.First Eq. (19b) is solved for the tensile field, then Eq. (19c) is solved for the shear field and finally Eq. (19a) is solved for the displacements u, which are used to update the history variables H I and H II .In the next staggered step, the phase-field equations are solved using the updated history variables and the staggered iterative scheme continues until convergence.As a stopping criterion for the staggered procedure, the residua of the three solution fields are compared against a certain threshold.The iterations are terminated after staggered step i if where In contrast to the classic two-field mixed-mode approaches and similar to Fei [22], the history variables H I and H II are updated depending on the dominating crack mode.If H I ≥ H II or s I (x) < 0.5 we assume that a mode I is present and only perform an update H I .Accordingly, if H I < H II or s II (x) < 0.5, we assume that a mode II crack is present and update H II .
3 Determination of the mode II fracture toughness In this section, the mode II fracture toughness is determined for two types of rock, namely Solnhofen Limestone (SPK) and Pfraundorfer Dolostone (PFD).The expected crack path is shown in Figure 3b), and consists of a straight shear crack connecting the two external notches.The two investigated rocks are analogue geothermal carbonate reservoir rocks.The Solnhofer limestone is very homgenous and very fine grained (0.1 -0.055 mm) and the Pfraundorfer Dolostone consists of 99% dolomite with Fig. 3: Schematic setup for determination of the mode II fracture thoughness (left), and expected crack pattern (right) following [7] .
small vugs [49].The experimental setup is based on the double-edge notched Brazilian disk (DNBD) tests presented in [7].In contrast to conventional mode II tests, the DNBD test features not only shear-based crack tip loading, but also the material failure is shear-induced.

Double-edge Notched Brazilian Disk (DNBD) tests
The experimental setup is schematically depicted in Figure 3 a).A Brazilian disk specimen containing two external notches of length a and thickness w is diametrically compressed at an angle α.Wooden plates are attached to the top and bottom of the specimen to induce the load.The use of flexible materials allows for an even distribution of the load and prevents local concentration of the stress, which can lead to fracturing at the loading points.At the same time it introduces an additional non-linear displacement.Increasing the load angle α, the mode I intensity factor K I increases, which results in smaller shear stresses and consequently a higher failure load.Bahrami [7] suggest to use an angle α in the range of 10°to 20°to prevent failure at the loading points, which can occur for very small as well as too large loading angles.Moreover, they state that the size of the ligament l which defines the distance between the two notches (cf. Figure 3) should be chosen depending on the disk radius such that 0.2 < l/R < 0.35.

Experimental results
Brazilian Disk samples of SPK and PFD with a radius of 47 mm and a thickness of 20 mm were prepared by cutting two external notches of length 37 mm.This results in a ligament of 20 mm and a ligament to radius ratio of l/R = 0.21, which lies within the suggested optimal range [7].Two different geometries were tested, one with notches of width 2.2 mm which were inserted by hand using a saw and a second geometry with notches of width 1.1 mm cut with a water jet cutter.While the water jet cutter resulted in perfectly aligned notches, the samples prepared by hand showed slight deviations in the geometry including misaligned notches and different notch length.For all experiments, a loading angle of α = 15°is chosen.The specimen are loaded until fracture with displacement controlled steps of 0.5 mm/min.The fracture process is recorded using a Photron Mini UX100 high-speed 10000 fps camera with a resolution of 1280×480 pixels.A spray pattern is applied to the specimen to perform a digital image correlation (DIC) analysis using the open-source 2D-DIC software Ncorr [9].Tab.2: Experimental results for K II .

Calculation of fracture toughness
The mode II fracture toughness is defined as the critical mode II stress intensity factor at failure.Following Bahrami [7], the mode II fracture toughness can be computed based on the failure load F , the geometrical dimensions of the specimen and the mode II stress intensity factor K * II following where t is the thickness of the disk, a is the notch length, and R is the radius of the disk.Bahrami [7] performed finite element analyses to determine the mode I and II stress intensity factors for different DNBD geometries.The relevant values for the selected geometry are summarized in Table 1.Following Equation 21, the mode II fracture toughness is calculated and the average of all experiments is taken.The results are summarized in Table 2.A mode II fracture toughness of 4.79 is calculated for SPK, while a lower K II of 3.99 is obtained for PFD.The mode I fracture toughness of SPK is taken from [42].For PFD it is calculated based on an empirical relationship stated in [52] and the material parameters determined in [44].As stated in Table 2, this yields a K I /K II ratio of 2.78 for PFD, and a ratio of 4.74 for SPK.
Type SPK PFD

Crack pattern and DIC analysis
In addition to the simple shear crack between the two notches as observed in [7], two additional crack patterns are obtained during the DNBD tests on SPK and PFD.As shown in Figure 4, two mixed-mode patterns are observed.Type A features two tensile wing cracks which initiate at the top and bottom notch, and emerge symmetrically.Further displacement contributes to the formation of a shear band connecting the notches, which leads to abrupt rupture of the specimen.Type B features similar tensile cracks.However, small deviations in the geometry lead to shear cracks emerging from the tips of the tensile cracks connecting to the notches.Type C is the pure shear crack observed by Bahrami [7].In Figure 4, right, the number of crack types observed for SPK and PFD are listed.While patterns A and B occur for SPK as well as PFD, pattern C is only observed for PFD.Consequently, both types of rock facilitate the formation of tensile wing cracks preceding the shear failure.For SPK, in all experiments two tensile wing cracks preceding the shear crack could be observed.In the case of PFD, in two out of ten experiments no tensile cracks initiated, and in four experiments only one tensile wing crack emerged.This observation fits very well to the K I /K II ratios determined.SPK has a higher K I /K II ratio than PFD, which explains the higher number of tensile cracks observed experimentally.

Numerical results
The phase-field model presented in Chapter 2 is calibrated for each type of rock using the computed mode II fracture toughness.In Section 4.1 the DNBD tests are analysed.A first study (Section 4.1.1) shows that the proposed model can reproduce the three crack patterns observed experimentally.To prove the correctness of the computed material parameters, specific geometries of the DNBD tests are simulated and compared with the experimental results in Section 4.1.2.As an outlook, a uniaxial compression test on a rare drill core is presented in Section 4.2 .

DNBD Experiments
The 2D setup for the mixed-mode simulation is shown in Figure 5, and the simulation parameters for the different crack types are listed in Table 3.The computational domain is initially discretized with 25 × 25 elements and refined towards the specimen boundary with a refinement depth of k = 1 and towards the notches with a refinement depth of k = 3.An adaptive refinement strategy based on the value of the phase-field is used, which refines the mesh in each staggered step based on the criterion s < 0.7 up to a depth of k = 4 for SPK, and up to a depth of k = 3 for PFD.To resolve the geometry, a quad-tree subdivision approach with a partitioning depth of s = 3 is used.Dirichlet boundary conditions are applied on two circular arcs of length 10 mm representing the contact zone between the wood and the sample on the top and, respectively, the bottom of the specimen (see Figure 5).While the displacements on the lower arc are fixed, negative displacements in the ydirection are applied on the top arc.An adaptive load stepping scheme is used based on the ideas presented in [26].Here, we use a step-size controller relating the step size change to the ratio of tolerance ε and current error ε i which computes the load step size in iteration i + 1 based on where κ = 1.1, u init = 5 • 10 −3 and u min = 5 • 10 −4 .The material parameters for both types of rock are listed in Table 4.While G c I is obtained directly from literature in the case of Solnhofen Limestone [42], for Pfraundorfer Dolostone it is computed from K I using the relation ( [32]).The mode II fracture energy G c II is obtained for both rocks by inserting the experimentally determined K II into the analogue expression Following the argumentation in [47], l 0 is treated as a material parameter.The length scale l 0,I associated to mode I failure is determined using the tensile strength σ t of the material, while the length scale l 0,II is computed based on the shear strength τ , following Here, we assume that the shear strength is equal to the maximum strength measured in the DNBD experiments.Consequently, we obtain different length-scale parameters for tensile and shear failure, namely l 0,I SPK = 0.259 and l 0,II SPK = 0.682 for SPK, and l 0,I PFD = 0.916 and l 0,II PFD = 0.656 for PFD.

Crack patterns
As explained in Section 3.2.2,three different crack patterns could be observed in the DNBD tests: the mixed tensile-shear crack (A), the circular mixed tensile-shear crack (B) and the pure shear crack (C).To assess the possibilities of the proposed model, geometrical and material parameters are varied to see if all crack types can be reproduced.For crack types A and C, a notch width of 1.1 mm is used, which corresponds to the geometry prepared by a water jet cutter.For Type B, a notch width of 2.0 mm is used and the notches are shifted perpendicular to their connecting line using an offset δ of size 1.0 mm.The latter corresponds to the geometry prepared using a hand saw.
The results of the phase-field simulations for both types of rock, SPK and PFD, are shown in Figure 6.
Here, for visualisation purposes a combined scale with a value of s = 1 on a tensile crack, and a value of s = −1 on a shear crack is chosen.As can be seen, the proposed model is able to correctly capture all three crack types.The mixed tensile-shear crack (A) features two tensile wing cracks emerging from the tips of the notches followed by a shear crack connecting the notch tips.This crack pattern can be reproduced for both SPK and PFD using the geometric parameters of the specimens prepared by a water jet cutter.Due to the higher K I /K II ratio of SPK, the tensile wing cracks propagate further before shear failure occurs compared to PFD.The circular mixed tensile-shear crack (B) exhibits similar tensile wing cracks, however, shear cracks develop between the tips of the tensile cracks and the respective, closest notch tips.Similar to Type A, it can be obtained using both SPK and PFD material parameters.Interestingly, the numerical results confirm that Type B occurs due to the imperfectly captured geometry when cutting the notches by hand saw.Due to the misaligned notches, the shear cracks start to initiate from the tips of the tensile wing cracks and evolve towards the notch tips, instead of the shear band forming in the center of the specimen as in the case of crack Type A. Here, the proposed model gives valuable insights concerning the formation and type of cracks, while experimental determination would require advanced experimental techniques and expensive equipment.The pure shear crack (C) shows a shear band connecting the two notches.In contrast to Type A, no or only minor tensile wing cracks have formed before shear failure occurs.Whereas crack types A and B are triggered by geometric differences in the specimen, the pure shear crack is obtained numerically by decreasing the K I /K II ratio and occurs for both geometries.For SPK, a K I /K II ratio of 2.23 is chosen, while for PFD a ratio of 2.12 yields a pure shear crack.Whereas the numerical result for PFD represents a pure shear crack, the result obtained for SPK shows small tensile cracks which initiate before shear failure occurs.In summary, the proposed model can reproduce all three experimentally observed crack patterns.A detailed study on the influence of both material and geometrical parameters on the resulting crack pattern, including stochastic analysis, might generate new valuable insights and is part of future research.

Validation
In this section, a detailed analysis of the mixed tensile-shear crack (Type A) is presented.To this end, only experiments yielding this specific crack pattern are considered.This corresponds to 5 SPK and 5 PFD experiments of the geometry prepared using a water jet cutter, which yielded mostly Type A crack patterns.In Figure 7, left, experimental and numerical crack patterns are compared for both types of rock.On the left, five different phases of crack propagation are evaluated for the experiments based on the high-speed camera recordings and compared with the simulation results.Here, the first column shows the displacement in the x-direction computed with DIC, while the second column shows the corresponding crack path.The computed phase-field crack paths are depicted in the third column (SPK) and the fourth column (PFD).In phase , a wing crack starts  to initiate at the notch where the displacement is applied.This behavior can be observed both in the experiments and in the numerical simulation.However, in the phase-field analysis, the second tensile crack initiates much earlier and not only after the first tensile crack has almost fully developed.This difference in behavior likely stems from the different boundary conditions applied in experiments and simulations.In the experiment, the disk can compress and sink into the soft wood, resulting in a change of the contact area over time.This behavior is not captured by the boundary conditions set for the numerical simulation.In phase , the propagation of the tensile cracks continues with increasing displacement.The growth of the tensile wing cracks steadily decelerates as soon as the wing cracks propagate up to the height of the opposite notch tip (phase ).Next, a shear band starts to develop in the center of the disk.In contrast to the tensile cracks, which initiate locally and grow from the crack tip with increasing load, the shear failure follows a different pattern.As can be seen, the shear crack initiates at the center of the disk and the associated damage covers a wider area.The damage increases gradually along the connection line between the notches (phase ) until failure occurs abruptly in phase .In summary, the numerical behavior, including the localisation of the shear band in the center of the specimen, agrees very well with the experimental observations.The fully developed crack patterns for SPK and PFD are shown in Figure 7, right.A direct comparison of the crack patterns of the two types of rock shows that the tensile wing cracks propagate further in the case of SPK.Due to the higher K I /K II ratio, the initiation of shear cracks occurs later, which enables the wing cracks to propagate beyond the notch tip of the opposite notch.This difference is also visible in the experimental crack paths.The computed load-displacement curves are shown in Figure 8 for SPK, top, and PFD, bottom).Due to the different boundary conditions in the experiment no direct comparison of experimental and simulated load-displacement curves is presented here.As the plastic deformations of the wooden plates which are used to transfer the load to the specimen can not be captured by the numerical model, in the following, the computed failure loads are compared against the averaged experimental failure load.For both types of rock, the crack phases -are marked.First, the force increases linearly until the tensile wing cracks start to initiate ( ).This results in a sudden drop in force which occurs at 10.18 kN in the case of SPK and at 14.37 kN in the case of PFD.Due to the different length scales for tensile cracks (l 0,I SPK = 0.259 and l 0,I PFD = 0.916), the loss in force is higher for PFD.The tensile wing cracks are not yet fully developed, when the force starts increasing again.The slope is almost linear until , when the shear band starts to develop.Failure occurs at the maximum bearable force of F max, SPK = 20.98 kN and F max, PFD = 16.8 kN, respectively.At this point ( ), the shear damage is already clearly visible.Once the shear band has fully developed ( ), the force drops to zero.As can be seen, the computed failure loads are in very good agreement with the experimental values.For both types of rock, the computed values lie within the range of values observed experimentally.The numerical failure loads show a relative deviation of 11.48% for SPK and of 2.98% for PFD from the respective averaged experimental failure loads.

Rare drill core
In this section, the proposed three-field model is applied to a complex, 3-dimensional crack scenario.Within the framework of the Geothermal-Alliance Bavaria, various laboratory experiments could be carried out on rare drill cores of the exploration well Moosburg SC 4 (MSC-4) ( [13,14,40]).Drilled in 1990 to a total vertical depth of 1585 m, the MSC-4 well is unique for being fully cored over the entire reservoir section of the Upper Jurassic carbonates (Malm aquifer) with a thickness of 453 m ( [11,12,33]).In the following, a uniaxial compression test on a rare drill core from a dolomitic part of the reservoir is presented.The low-porosity dolostone sample shows crystalline  sizes of 0.25 mm to 1.0 mm as well as small and large vugs.To obtain the drill core's exact geometry, the cylindrical sample with a height of 98.3 mm and a diameter of 49.7 mm was CT-scanned with a resolution of 0.11 mm × 0.11 mm × 0.11 mm.A uniaxial compression test was performed using a displacement-controlled test speed using a displacement rate of 0.06 mm/min.To be able to analyze the experimental crack pattern, the compression test was recorded using a Phantom Flex4K highspeed camera with 2000 fps in full HD resolution.The experimental setup and observed crack pattern is shown in Figure 9, left.Upon loading, a vertical crack initiates on the upper side of the large, central pore on the front side of the specimen (Figure 9, ).In addition, cracks emerge which connect to the smaller pore on the right side of the specimen as well as the lower side of the specimen (Figure 9, ).The vertical crack continues to grow upwards until it reaches the top side of the drill core.The premature damage on the left side of the large pore leads to the development of smaller cracks which will connect to a pore on the back of the core sample, as seen in Figure 9, .Here, failure occurs and the right half of the specimen is blasted off.The failure pattern is dominated by the pores inside the rock, which trigger the initiation of the cracks.

Simulation Setup
The geometry and boundary conditions for the numerical simulation are depicted in Figure 9, right, and the simulation parameters are listed in Table 5.The boundary conditions are set as follows: y-displacements are fixed on the top surface, while a positive displacement is applied in y-direction.
Steps of size 0.02 mm are applied until a total displacement of 0.2 mm, at which the step size is decreased to 0.005 mm.The geometry is represented using FCM.To this end, the core sample is embedded into a Cartesian mesh with 44 × 92 × 44 elements and integration is performed using an octree sub-division approach with partitioning depth s = 3.This results in a total number of 2 168 105 DOFs for the three-field system.Due to the limited access to the material of the rare drill core, the parameters for this dolostone could not be determined experimentally.Consequently, it is assumed that the MSC-4 dolostone behaves similar to the PFD and material parameters are taken from Table 4.However, the mode II fracture toughness needs to be chosen differently, as the value determined for PFD results in premature cracking and a crack pattern which does not relate to the one observed experimentally.Based on a parameter study we choose G cII = 5 • 10 −3 kN/mm, for which we obtain good agreement in both the failure load as well as the observed crack pattern.The higher value of G cII can be attributed to the three dimensional setting, in which not only mode I and mode II, but also mode III cracks occur.However, the extension of the proposed two-field problem to account for mode III cracks is the subject of future research and is beyond the scope of this work.The simulation is performed on the SuperMUC-NG cluster at the Leibniz Supercomputing Centre using 24 nodes with 48 cores and 96 GB memory per node.The hybrid MPI and OpenMP parallelisation of the three-field phase-field problem is based on the framework presented in [29,30].

Results
In Figure 10  of the central pore leading to complete failure of the sample.For a detailed comparison of the final crack pattern, the experimental cracks are highlighted in Figure 11 and contrasted with the numerical result.A common feature is the vertical tensile crack which initiates on top of the internal pore and propagates upwards.In contrast to the experiment, where a nearly straight vertical crack is observed, the computed crack tends to lean towards the outer surface of the specimen.This behavior is related to the boundary conditions.Firstly, as can be seen in the experimental setup (cf. Figure 9), the top and bottom surface of the drill core are not completely parallel.This results in a real displacement which differs from the pure in-axis displacement applied in the simulation.Moreover, as a consequence of the phase-field formulation, the crack is not able to penetrate the Dirichlet boundary.Instead, it isrepelled from the top and bottom surface of the core sample.Therefore, the part on the right side of the core sample that falls off in Figure 9, , is smaller in the simulation than in the experiment and shaped differently.The vertical cracks connecting the larger pores to the bottom side of the specimen can not be reproduced in the simulation.Similar to the experiment, shear and combined tensile-shear cracks connecting the different pores are visible outside the specimen.Since the fracture pattern could only be recorded from one side, it is difficult to judge if the final failure occurred due to a shear crack.However, the visible experimental crack pattern is captured remarkably well, and the numerical simulation generates interesting insights, including the fact that the vertical tensile crack initiates on top of the internal pore, as highlighted in Figure 9, step v 1 .In Figure 12

Conclusion
In this contribution, a three-field phase-field model for the simulation of mixed-mode fracture in rock is presented.Separate scalar phase-field variables associated to mode I and mode II failure are introduced, and the two phase-field equations are implicitly coupled through the degradation of the simulation, p=2 Fig. 12: Comparison of the experimental and computed load displacement curves.While the experimentally measured fracture force is 147.0 kN, the numerical simulation predicts failure at 142.31 kN with corresponds to a deviation of 3.2 %. material in the elastic equation.By introducing separate length scales for the mode I and the mode II problems similar to [22], the major strength of the model lies in its ability to account for different tensile and shear strengths of the material.The framework is easy to implement and flexible, as it allows the choice of different splits and degradation functions.By clearly distinguishing between tensile and shear cracks it facilitates the analysis of complex fracture patterns.
To validate the three-field approach, the model was calibrated for two types of rock, Solnhofen Limestone and Pfraundorfer Dolostone.The mode II fracture toughness for each type of rock was determined experimentally using double-edge notched Brazilian disk tests.The simulations of the DNBD tests demonstrate that the proposed model can reproduce the three crack patterns observed experimentally: a mixed tensile-shear crack, a circular mixed tensile-shear crack and a pure shear crack.Moreover, the computed failure loads agree very well with the averaged experimental results with a deviation of 11.48 % for SPK and a deviation of 2.98 % for PFD.To test the applicability of the model for realistic 3D fractures of complex shaped specimen the three-field model was applied to a uniaxial compression test on a rare drill core.For a detailed analysis of the crack patterns, the experiment was recorded using a high-speed camera.The exact geometry of the dolostone sample was extracted from a CT-scan.The computed crack pattern captures the most characteristic fractures observed experimentally.The recorded load-displacement curve can be reproduced with good agreement showing a deviation in the failure load of 3.2 %.The deviations can be explained by uncertainties in the boundary conditions as well as the diverging material parameters of the rare drill core.The example demonstrates the ability of the model to reproduce complex, 3-dimensional crack patterns in rock and its potential to generate valuable insights in the field of mixed-mode fracture.
, left.The set of discrete cracks Γ c is split into a set of tensile cracks Γ c I associated to mode I failure and a set of shear cracks Γ c II associated to mode II failure with Γ c = Γ c I ∪ Γ c II and Γ c I ∩ Γ c II = ∅ .The domain boundary ∂Ω consists of two non-overlapping parts Γ D and Γ N on which Dirichlet and Neumann boundary conditions are prescribed.

Fig. 4 :
Fig. 4: Different crack patterns observed in the experiments: mixed tensile-shear crack (A), circular mixed tensile-shear crack (B) and pure shear crack (C) (left), and frequency of occurrence for both types of rock (right).

Fig. 5 :
Fig. 5: Setup of the computational domain (left), and mesh with initial refinement and boundary conditions (right).

Fig. 6 :
Fig. 6: Comparison of experimental and simulated crack paths for the three crack types obtained for SPK and PFD: the mixed tensile-shear crack (A), the circular mixed tensile-shear crack (B) and the pure shear crack (C).

Fig. 7 :Fig. 8 :
Fig. 7: Comparison of experimental and numerical results for crack type A. From left to right: displacement in x-direction obtained by DIC, experimental crack pattern, computed phasefield for SPK and computed phase-field for PFD (left).Computation domain with mesh refinement for the fully cracked specimen for SPK (i) and PFD (ii) (right).

Fig. 9 :
Fig. 9: Experimental setup and recorded crack pattern (left), and simulation setup (right) based on the CT-scanned drill core.

Fig. 10 :
Fig. 10: Simulated crack pattern for different displacement steps.For visualisation purposes the cracks are shown as iso-volumes of the phase-field, i.e. s I ≤ 0.03 for the tensile cracks depicted in blue, and s II ≤ 0.03 for the shear cracks depicted in red.Two different views are shown: a front view corresponding to the viewing angle of the camera in the experiments (top row), and a side view (bottom row).The red and orange line in the first column illustrate the orientation of the different views.

Fig. 11 :
Fig. 11: Comparison of the experimental crack pattern with manually highlighted cracks (left), and computed tensile and shear cracks obtained with the three-field model (right).
, the experimental load-displacement curve is shown and compared with the computed result.The experimental curve clearly shows a brittle fracture behavior predicting failure of the sample at a force of 147.0 kN and a displacement of 0.24 mm.The computed curve closely follows the slope of the experimental curve until a displacement of 0.224 mm when the maximum load of 142.31 kN is reached.At this point, the tensile crack starts to develop which results in a drop in force.Once the tensile crack has stabilised, the curve slowly starts to ascend again.Shear damage accumulates which results in complete failure at a displacement of 0.256 mm.The deviation of the computed failure load corresponds to 3.2 % of the force measured experimentally.