Multiphase-field modelling of crack propagation in geological materials and porous media with Drucker-Prager plasticity

A multiphase-field approach for elasto-plastic and anisotropic brittle crack propagation in geological systems consisting of different regions of brittle and ductile materials is presented and employed to computationally study crack propagation. Plastic deformation in elasto-plastic materials such as frictional, granular or porous materials is modelled with the pressure-sensitive Drucker-Prager plasticity model. This plasticity model is combined with a multiphase-field model fulfilling the mechanical jump conditions in diffuse solid-solid interfaces. The validity of the plasticity model with phase-inherent stress and strain fields is shown, in comparison with sharp interface finite element solutions. The proposed model is capable of simulating crack formation in heterogeneous multiphase systems comprising both purely elastic and inelastic phases. We investigate the influence of different material parameters on the crack propagation with tensile tests in single- and two-phase materials. To show the applicability of the model, crack propagation in a multiphase domain with brittle and elasto-plastic components is performed.

fractures which impact the strength and other mechanical properties of these rocks [3][4][5]. There is an extensive literature of experiments dealing with crack initiation, propagation and coalescence in rocks or rock-like materials under different loading conditions [6][7][8][9][10][11][12][13]. Depending upon the geophysical factors (e.g. pressure, temperature, mineralogy, strain rates and grain size), different rocks exhibit diverse modes of failure, ranging from linear elastic brittle fracturing [14,15] to elasto-plastic deformation followed by crack propagation [16,17]. In order to model the pressure-sensitive plastic yielding in geological materials and porous media, different plasticity models such as Mohr-Coulomb [18] and Drucker-Prager [19] are designed. Macroscopic plastic yielding occurs due to frictional sliding of particles or formation and compression of voids in porous material at microscale. Numerical works of Zreid and Kaliske [20,21] demonstrate the sound predictive capacity of the Drucker-Prager model in recovering the experimentally determined mechanical behaviour of concrete that comprises a matrix of rocky material and a binder. On the other hand, Griffith's criterion [22] is among the most popular theories to describe fracturing in linear elastic materials. It states that a crack propagates when the energy release rate exceeds a critical value necessary to create the crack surfaces. This critical value is known as crack resistance or fracture toughness of the material. Irwin's modification [23] extended the theory to ductile materials and suggested that a plastic zone is formed near the crack tip when a load is applied, which further grows with the increase in loading. Moreover, he introduced the socalled stress intensity factors in order to describe the stress and displacement fields near the crack tip. The critical values of the stress intensity factors for different modes are mathematically related to the crack resistance and are used as engineering design parameters for fracture-tolerant materials, employed for the construction purposes (e.g. buildings and bridges).
In order to computationally deal with crack propagation in geological materials and porous media, different numerical approaches may be used, e.g. boundary elements method (BEM) [24][25][26], peridynamics [27,28], discrete element method (DEM) [29,30] and extended finite element method (XFEM) [31,32]. Using BEM, Zhang et al. [24] investigated the hydraulic fractures in an existing network in 2D. Furthermore, Wu and Olson [25] and McClure et al. [26] extended and applied the BEM to 3D fracture networks. However, the method, although computationally cheaper than others, faces problems in dealing with heterogeneous and anisotropic materials, as the discretisation is only limited to boundaries. Peridynamics was utilised by Ha and Bobaru [27] to analyse dynamic crack branching in brittle materials and was coupled with fluid flow by Ouchi et al. [28] to model fluid-driven fractures in porous media.
Virgo et al. [30] employed DEM in order to investigate crack propagation in a heterogeneous vein-hosted rock structure under different boundary conditions. However, due to high computational costs, the method is not wellsuited for performing large-scale numerical investigations. Wang et al. [32] investigated hydraulic cracks in existing fracture networks using XFEM. This method, which is an extension of the finite element method, computes the displacements by modifying the shape functions in the elements exhibiting cracks. For an elaborate review of the numerical methods and advances in computational fracture mechanics, interested readers are referred to [33].
All the above-mentioned models treat crack surfaces as sharp discontinuities in materials and require explicit tracking of the interfaces. Design and development of robust algorithms based on these methodologies and their numerical implementation become tedious in the treatment of complex crack geometries (e.g. crack branching and merging, multiple cracks), especially in addressing these processes in 3D. The phase-field method, already well established in the materials science community for modelling phase transition processes involving multiple components and phases (e.g. solidification, crystal growth; for further examples, see review articles [34][35][36][37][38][39]), has emerged as a robust methodology for solving problems of crack propagation in brittle and ductile solids [40][41][42][43][44][45]. Unlike the aforementioned sharp interface approaches, the phase field models treat fracture surfaces as diffuse regions, where a scalar field describes the state of a material and varies smoothly from the fully damaged material state to the fully intact material. The intermediate diffuse regions are interpreted as partially damaged material states, which correspond to the evolution of micro-cracks and micro-voids. Furthermore, the approach supersedes the need of tracking the interfaces explicitly through any external fracturing criteria, making it a convenient and computationally efficient approach in dealing with geometrically complex crack branching and merging problems. One of the first works of Bourdin et al. [41] explored the phase-field method for describing fracture. A treatment of the process taking into account thermodynamic aspects was presented by Miehe et al. [42]. Different investigations of ductile crack propagation using the phase-field method are presented in [46][47][48][49][50][51]. Choo and Sun [50] and Kienle et al. [51] utilised a modified Drucker-Prager yield criterion to model fracture formation in frictional material. Despite their achievements, these studies were limited to failure prediction in homogeneous materials. Investigations addressing fracture propagation in heterogeneous materials were initiated by Spatschek et al. [52], and further developed by [53][54][55] for multiphase systems. Here, one topic was the impact of fracture toughness of distinct polycrystal domains and their grain boundaries on the crack path. Furthermore, the impact of anisotropy in the brittle polycrystalline solids was investigated by [56][57][58][59]. In the recent works (e.g. [60][61][62][63]), flow physics and fracture mechanics were coupled to investigate hydraulic fracturing in porous media. The phase-field solution shows a sensitivity depending on the choice of the length parameter. There have been approaches setting the length parameter as a physical material property or presenting techniques which ensure insensitive solution results [64][65][66].
In this work, we present a multiphase-field model, capable of describing fracture formation processes in heterogeneous multiphase materials, where separate regions of (i) purely brittle and anisotropic material and (ii) frictional and porous material are present. With the above-mentioned capabilities, the model may be well suited for various civil engineering and geotechnical applications (e.g. computational geomechanics, subsurface energy production). In the present work, we limit the fracture formation under pure mechanical loading conditions. However, we remark that the presented framework allows for further extensions (e.g. coupling with temperature or fluid flow equations), and for simulating the process of fracture formation under a wide range of boundary conditions. The content of this work is structured as follows. In Section 2, the multiphase-field model for crack propagation in elasto-plastic and anisotropic elastic materials is adapted from the works of Prajapati et al. [67] and Herrmann et al. [68]. Section 3 gives a description of the numerical solution scheme of the implementation. In Section 4.1, we consider the Drucker-Prager plasticity model for describing the plastic behaviour in inelastic solid phases. This plasticity model was considered owing to its relative simplicity in reflecting the characteristics of the soil and rock materials with few parameters. Through the numerical examples of twoand four-phase specimens under uniaxial mechanical loading, the model implementation is validated by comparing numerical results with interface solutions. Crack propagation in single-as well as two-phase materials is presented in Section 4.2. To demonstrate the capabilities of the model, crack propagation was simulated in multiphase structures, with the material parameters of sandstone and a brittle rock for different solid phases, in Section 4.3. Section 5 concludes the presented work and gives an outlook on possible future research directions.

Model formulation
To simulate elasto-plastic crack formation in geological materials and porous media, the model for anisotropic crack propagation in brittle grains introduced by Prajapati et al. [67] is adapted to conduct simulation studies of elasto-plastic crack propagation in a multiphase system comprising different brittle and ductile regions. Schneider et al. [70][71][72] introduced an approach to fulfil the mechanical jump conditions in the diffuse interface to obtain correct stress and strain fields in multiphase domains within infinitesimal and finite deformations. Herrmann et al. [69] extended this work to include the inelastic material behaviour based on J2 plasticity (see e.g. [73]). Their homogenisation approach [69] is, in principle, independent of the utilised plasticity model. Therefore, we adapt this approach with the Drucker-Prager plasticity model for describing the pressure-sensitive flow behaviour of geological materials and porous media (e.g. soil or rock) [19]. For a detailed discussion of the mechanical jump conditions in the diffuse interface, for small as well as large deformations, interested readers are referred to the works of [69][70][71][72]. The crack model originates from Schneider et al. [55].

Multiphase-field model
We consider a set of N + 1 phase-field order parameters where each phase-field order parameter φ α (x, t) ∈ [0, 1] describes the volume-fraction of a particular solid phase or crack phase α ∈ {1, . . . , N, c} at position x and time t. The interface between different order parameters is characterised by a diffuse region of finite width. The magnitude of order parameter φ α (x, t) varies smoothly from 1 inside the phase α to 0 outside through the diffuse interface region. Solid phases correspond to different, physically separated regions with diverse material properties that may characterise the different rock types. The crack phase, represented by φ c , describes the material degradation, i.e. φ c = 0 for fully intact material and φ c = 1 for a fully damaged material. At each point, the sum of all occurring phase fields is constrained by: Based on Griffith's criterion [22,40], the total free energy functional F of a domain of volume V is given by: where G c represents the effective crack resistance, f el denotes the effective strain energy density and f pl is the effective plastic dissipation, in which ε el is the elastic strain tensor,ε pl is the accumulated plastic strain and ∇φ c is the gradient of φ c . The crack resistance maps the energy of an emerging crack by using the gradient term |∇φ c | 2 and a potential term ω c (φ c ) in a regularised manner.
In the present work, we chose a single obstacle potential ω c (φ c ) = k ω φ c over the conventional well potential, as it allows much smaller interface widths, resulting in lower computational costs [74]. The constants k = 3π 2 /64 and k ω = 8/π 2 2 are calculated on the basis of the regularisation of the crack resistance in an integral manner along interface normal direction and equating it with the sharp interface solution [74]. The length-scale parameter c controls the crack interface width. The temporal and spatial evolution of the crack phase-field is described by the Allen-Cahn equation [75] in terms of a variational derivative of the free energy functional which is given as: where M is a constant positive mobility. Substituting the expression of the free energy functional (2) in the above (3) results in: As solid-solid transformations are not considered in the present work, the evolution of solid phases is formulated as: with the interpolation function h α s for the solid phases given as the normalised phase fraction: This formulation ensures the condition N α h α s (φ) = 1. As soon as the crack phase field φ c reaches a given critical value φ c,crit , a fully damaged material state is considered. In this case, the crack phase field is set to 1 and all solid phase-fields are set to 0. Moreover, effects of crack healing are neglected and only the growth is considered; therefore, φ c ≥ 0 is enforced [76]. The effective strain energy density f el and the effective plastic dissipation density f pl are defined as the volumetric averages of the phase-inherent strain energy densities f α el and the phase-inherent plastic dissipation densities f α pl , respectively, reading: Detailed expressions of the phase-inherent strain energy densities f α el and the phase-inherent plastic dissipation density f α pl are discussed in Section 2.2. Furthermore, the effective crack resistance is computed as a volume average of the phase-inherent crack resistances: For investigations of anisotropic brittle crack propagation in polycrystalline quartz sandstones, Prajapati et al. [67] introduced an anisotropy of the crack resistance which is given by: The crack resistance G α c,0 is reduced by the anisotropy factor F α a,x = F α a,y = F α a,z in the x-direction of the rotated coordinate system and represents an ellipsoid in 3D space. The rotated normal vector of the diffuse crack interface n α c (∇φ c ) is calculated as the rotation of the normal vector to the crack interface ∇φ c /|∇φ c | using the rotation matrix Q α : The rotation matrix describes the anisotropy orientation of a solid phase. For materials exhibiting isotropic crack resistance follows F α a,x = F α a,y = F α a,z .

Elasto-plastic multiphase-field model with Drucker-Prager plasticity
The mechanical calculations are performed within the theory of small deformations for the sake of simplicity and to limit the computational costs. However, the presented model is not restricted to small deformations and can be extended to account for finite strains [71]. The total strain ε is given as ε = 1 2 (H + H T ), where H is the displacement gradient. Furthermore, the elastic ε el and plastic strain ε pl contributions are additively retrieved as ε = ε el + ε pl .
Herrmann et al. [69] introduced an elasto-plastic multiphase model according to the mechanical jump conditions, where the inelastic deformations were calculated based on J2 plasticity theory (see e.g. [73]). The homogenisation approach for computing the phase-inherent stresses σ α and plastic strains ε α pl in the interface regions results in well-defined stress and strain fields between brittle and elasto-plastic phases. We incorporate the Drucker-Prager plasticity theory [19] into the homogenisation approach, which describes plastic deformation in geological materials and porous media.
The pressure-sensitive yield function of the Drucker-Prager plasticity model is given by: where c α and β α are the phase-inherent cohesion and friction angle, respectively. It is used to calculate the phaseinherent plastic strains ε α pl , depending on the phase-inherent stresses σ α . The deviatoric part of the stress tensor σ α is given by (σ α ) = σ α − p α (σ α )I in terms of the pressure p α (σ α ) = 1/3(σ α 11 + σ α 22 + σ α 33 ) and the second-order unit tensor I . The second invariant of the deviatoric stress is given as The cohesion c α of phase α is dependent on the accumulated plastic strainε α pl . For the case of linear hardening, the cohesion is given as c α = c α 0 + H αεα pl , where c α 0 is the initial cohesion and H α is the hardening modulus. In the principal stress space, the Drucker-Prager flow surface forms a cone where the apex lies on the positive hydrostatic axis, and its position depends upon the chosen values for the cohesion and the friction angle. Figure 1 illustrates the variation of the yield surface for different values of cohesion c α and friction angle β α . If the stress state at a material point is outside the cone, the inequality (12) is violated and plastic flow occurs.
The evolution of plastic strains is then given by the associative flow rule that reads: where N α determines the direction of plastic flow and is perpendicular to the yield surface. A predictor-corrector return-mapping algorithm is used to compute the plastic strains. When the predicted stress state falls inside the complementary cone of the yield surface (green region in Fig. 1(a), the stress state is mapped to the apex of the cone (turquoise in Fig. 1(a)). Otherwise, it is mapped to the cone surface. Due to the pressure dependency of the yield criterion, plastic deformation can occur even within a purely hydrostatic stress state. The elasto-plastic multiphase-field model fulfils the jump conditions of both, stress vector t = σ [n] as well as the displacement gradient H , given by the sharp interface theory [77]. With the conservation of linear momentum at a material singular surface, the continuity of the stress vector t follows. Here, n αβ is the singular surface normal vector between a α and β phase. Hadamard's compatibility condition [77] is given for small-strain theory by It is a kinematic compatibility condition which ensures that the displacement gradient is continuous in the tangential directions of a material singular surfaces. Here, a αβ is an arbitrary vector and ⊗ is the dyadic product. In order to determine the jumps of the phase-inherent displacement gradients (H α − H β ) at the α-β transition along the n αβdirection in a diffuse multiphase region with N occurring solid phases, N −1 unknown vectors, i.e. a 12 , a 13 , . . . , a 1N , have to be evaluated. These unknown vectors are set with respect to the phase with the largest volume fraction (phase 1 in this case) and are arranged in an N − 1 tupleâ given bŷ The effective displacement gradientH (φ) is expressed as the volumetric interpolation of the phase-inherent Therefore, the displacement gradient of phase 1 is written as Deformation gradients for the remaining phases can be expressed as Utilising Hooke's law σ α (â) pl ] and the phase-inherent displacement gradients (in Eq. 18) the unknown vectors (in Eq. 16) can be calculated by solving the following system of equations based on the balance of linear momentum on material singular surfaces (14): Here, C α denotes the phase-inherent stiffness. Due to the plastic strains, the system of equations is highly non-linear, and is solved by using the Newton-Raphson scheme: The phase-inherent strain energy density for elastic and isotropic materials is additively decomposed into positive and negative parts: where λ α and μ α are the phase-inherent Lamé coefficients, ε α i is the phase-inherent elastic principal strain and ε α i ± = (ε α i ± |ε α i |)/2. This spectral decomposition was introduced by Miehe et al. [43] in order to avoid non-physical crack propagation in compressive stress states. The crack interpolation function h c (φ c ) degrades the strain energy density of solid phases and is given by [55]: In contrast to [68], only one interpolation function is used for the brittle and the elasto-plastic material because crack propagation in the elasto-plastic material is possible irrespective of the occurrence of plastic deformation. The static balance of linear momentum without body forces, ∇ ·σ = 0, is solved for the displacement field u with the degraded volume-averaged stresses: With the solution forâ from Eq. 19 at hand, the phaseinherent stresses can be calculated as follows: Unlike the strain energy formulation in Eq. 21, the full stress tensor is degraded by the interpolation function h c (φ c ), similar to [68]. Therefore, this formulation is considered as a hybrid model [47]. The plastic dissipation density for each phase α [78] is given by: in terms of the phase-inherent hardening modulus H α and the accumulated plastic strainε α pl .

Numerical treatment
In the present work, we utilise an equidistant, orthogonal grid for the numerical calculations. Using this grid, the curvatures can be sufficiently described through the diffuse interfaces; thus, no complex discretisation is necessary. The values of the phase fields are evaluated at the central positions of the grid points, while the mechanical fields are computed in a staggered manner implicitly. The phase-field evolution (3) is numerically solved using an explicit Euler scheme for the time derivative and second-order accurate central difference scheme for the spatial derivatives [55].
The rotated staggered grid, where the stress and strain fields are calculated, is equivalent to finite element (FE) discretisation with linear elements using full integration. The phase fields and the mechanical fields serve as input for each other. Based on elastic assumption, the initial predictor displacement fields are computed by solving the static conservation of linear momentum ∇ ·σ = 0 (23), using the values of the plastic fields from the previous time step. The plastic strains are computed using a returnmapping algorithm for Drucker-Prager plasticity (see [73] for algorithmic equations), where each phase is treated separately. After the evaluation of plastic strains, the jump conditions in the diffuse interface region might be violated, and thus are adjusted by solving (19). The calculations of plastic strains and adjustment of the displacement fields to satisfy the jump conditions are performed iteratively until a steady state is reached. For the calculation of ∇ ·σ = 0, a biconjugate gradient stabilized (BICGSTAB) algorithm with a Jacobi preconditioner is utilised, which shows a slightly better convergence than the standard conjugate gradient (CG) algorithm. The mechanical fields serve as input for the evolution of the crack phase field (3). The crack phase field φ c is updated such that it fulfilsφ c ≥ 0 after each time step. After the crack phase field reaches a steady state, the next increment in the mechanical load is applied. A material point is considered fully cracked above a critical value φ c,crit , while a steady-state crack phase field is assumed to be reached when ∂φ c (x, t)/∂t < [∂φ c (x, t)/∂t] steady state . As a result of this procedure, the

Representative numerical examples
Fracture formation in heterogeneous multiphase materials, including anisotropies and inelastic effects, is a highly complex phenomenon. Therefore, a systematic study is conducted, where each model ingredient is validated in a step-by-step procedure. In Section 4.1, we validate the implemented Drucker-Prager plasticity model coupled with the homogenisation approach by comparing the simulation results for the loading of two-phase and fourphase specimens obtained from PACE3D with the sharp interface results of the commercial FEM solver ABAQUS. In order to ensure a direct comparison, the material parameters, boundary conditions and the integration scheme are kept identical for both the solvers. In Section 4.2, crack propagation in ductile material and impact of various parameters on the material response is investigated. Finally, in Section 4.3, we demonstrate the capabilities of the model through simulation of crack propagation in multiphase geological structures comprising of brittle and ductile regions. Furthermore, we discuss the factors influencing the crack pathways in such systems.

Model validation: Drucker-Prager plasticity with homogenisation approach
In ductile sandstones comprising of different regions with varying plastic properties, stresses and plastic strains have to be treated in accordance with mechanical jump conditions in the diffuse interface and junctions. As a starting point, we compare our plasticity implementation in multiphase systems with the results of ABAQUS, through the following examples (in Sections 4.1.1 and 4.1.2). The simulations with PACE3D and ABAQUS are performed on an equidistant orthogonal grid, under the same boundary conditions. The mechanical load is applied quasi-statically in one time-step. Linear elements with full Gaussian integration are chosen in ABAQUS, which is equivalent to the finite difference discretisation with a rotated staggered grid in PACE3D.

Benchmark case 1: serial and parallel combinations of a two-phase specimen
We consider a computational domain comprising of two ductile phases (see Fig. 2(a, b)). The computational domain was generated by the following preprocessing procedure: (i) filling the domain with two different solid phases, followed by (ii) a step to create diffuse solid-solid interface.
To study the influence of the diffuse solid-solid interface width, the length-scale parameter, s , was varied in order to generate different initial simulation settings. We refer to them as phases 1 and 2, where the material parameters are represented with the indices 1 and 2. The phase-field parameter φ 1 determines the presence of the bulk and the diffuse interface regions. For both phases, Young's modulus E 1 = E 2 = E 0 = 25 GPa and Poisson's ratio of ν 1 = ν 2 = 0.33 is assigned. In order to obtain different plastic strains, the cohesion in phase 2 is chosen to be c 2 = 1.2 c 0 , whereas the cohesion of phase 1 is c 1 = c 0 = E 0 /1000. Moreover, the linear hardening modulus H 1 = H 2 = E 0 /10 and the friction angle β = 15 • are chosen. Based on the applied loading conditions, two simulation setups are considered: serial and parallel combinations (see Fig. 2(a, b)). In a single step, orthogonal displacements u 1 and u 2 = u 1 /5 are applied in the serial and parallel combinations, respectively, which amounts to the same total strain in the case of a homogeneous and purely elastic material for the two loadings. On the remaining boundaries, the orthogonal displacements are set to 0. Simulations are performed for different interface widths of the solid phases, i.e. s = {0, 3Δx, 5Δx}, with Δx being the grid spacing in the x-direction. Here, s = 0 implies a sharp interface between the two solid phases, without utilising the second step of the above-mentioned preprocessing procedure where only a diffuse interface between the crack and corresponding solid phase is present in one computational cell. For the numerical examples, where multiple solid phases and the crack phase are present, the choice s = 0 essentially implies that, at any given computational cell, the crack phase field and only one solid  Figure 2 (c, d) depicts the variation of the normalised von Mises equivalent stress σ vM /σ vM,max along the red-dashed line (in Fig. 2(a, b)) for different interface widths and the sharp interface solution from ABAQUS. Due to the distinct values of cohesion of the two phases, a different amount of plastification occurs in each phase, resulting in different stress distributions. In the sharp interface solution of ABAQUS (blue curve), the normalised von Mises equivalent stress sharply jumps at the two-phase boundary. For the diffuse interface solution (dark green and red lines from PACE3D), a smooth and continuous increase over the whole interface width is observed, following the chosen interpolation function h α s . A higher value of the length parameter s leads to widening of the diffuse interface, and therefore, results in an expanded transition zone of the von Mises stress (see zoomed inset picture in Fig. 2(c)), without influencing the values in the bulk regions. Similarly, the normalised accumulated plastic strainε pl /ε pl,max along the red dashed line (in Fig. 2(e, f)), for the respective serial and parallel combinations, exhibits a continuous and smooth decrease ofε pl in the diffuse interface solution, and a jump at the phase boundary in the sharp interface solution. It is worthy to mention that the von Mises stresses and the accumulated plastic strains in the bulk phases are invariant for the different widths of the diffuse interface and also match with the sharp interface solution.

Benchmark case 2: four-phase specimen
In this section, we consider a simulation domain comprising four different ductile phases (see Fig. 3(a)). Such a specimen is expected to exhibit more complex stress states compared with the previous cases. The values of material parameters (i.e. E and ν) are kept identical in all the phases and are given in  Figure 3 (b-e) depict the plots of normalised von Mises equivalent stress σ vM /σ vM,max and the normalised accumulated plastic strainε pl /ε pl,max along the dashed lines -, highlighted in Fig. 3(a). The scaling factors (maximum values) along all lines are distinct. For the lines intersecting the binary interfaces, i.e. and , the values of von Mises stress and accumulated plastic strain in the bulk regions, obtained from the PACE3D simulations, are invariant with respect to the interface width and further exhibit good quantitative agreement with the results of ABAQUS. Moreover, in the diffuse interface region, a continuous and smooth transition is seen (see zoomed inset pictures in Fig. 3(b)), in accordance with the results for the These results advocate the accuracy of the implemented plasticity model with the homogenisation approach in precisely computing the stresses and strains in the bulk phases, while exhibiting good agreement in the interface regions.

Elasto-plastic crack propagation
We investigate the elasto-plastic crack propagation in ductile geological materials and porous media based on Drucker-Prager plasticity in Section 4.2.1. Using a twophase specimen comprising brittle and ductile phases, the brittle to ductile transition at the interface during crack propagation is discussed in Section 4.2.2.

Single-phase ductile specimen with pre-existing crack
We consider a square computational domain (151Δx × 151Δx) having a pre-existing crack (length = 68Δx) which shares a diffuse interface with a ductile solid phase (see Fig. 4(a)). At the lower boundary, the displacements in all directions are set to 0 (u = 0), while the right and left boundaries are rendered stress free (σ [n] = 0). On the upper boundary of the domain, an orthogonal displacement u is applied in an incremental manner. As discussed in Section 3, the load increment is applied only after the evolution of the crack phase reaches a steady state. If not specified differently, the material parameters are chosen as follows: E 0 = 25 GPa, c 0 = E 0 /2500, H 0 = E 0 /5, β = 40 • , G c,0 = 1 J/m 2 and c = 5Δx. We remark that the chosen value of the hardening modulus is higher in orders of magnitude compared with engineering alloys. However, given the small-strain assumption of the present model, and non-availability of precise data of the above parameters for the geo-materials along with experimental benchmarks, the present numerical examples are illustrated only for the purpose of demonstrating the model capabilities, and should not be considered precisely quantitative in nature. Thus, in order to generate a sufficiently high driving force for a macroscopic fracture propagation, a high hardening modulus is chosen. Figure 4 (f-i) depict the load-displacement response for varying values of cohesion c, hardening modulus H , friction angle β and the crack interface length parameter c , respectively, while keeping the other parameters fixed. Force and displacement are plotted as normalised quantities, i.e. F /F max and u/u max , scaled by their maximum absolute value. After the loading is sufficiently high, and the yield condition is violated, inelastic deformation starts. The plastic strains evolve in such a way that a so-called active plastic zone is formed around the crack tip (see Fig. 4(b)).
As soon as a cell is fully cracked, the plastic fields are set to 0, resulting in a plasticity field around the crack boundaries. This zone is described in literature as a relic plastic zone [4]. For the chosen range of material parameters, only a small region around the crack tip plastifies. As no precise material properties and loading conditions of experimental samples are available, we remark that reproduction and comparison of plastic zones in simulations and experiments are not straightforward. Nonetheless, the occurrence of an active and a relic plastic zone in the simulations is in qualitative agreement with experimental observations [4].
A sensitivity analysis of the load-displacement response for different values of ∂φ c (x, t)/∂t and φ c,crit is depicted in Fig. 4(c and d). It is observed that, when the chosen value ∂φ c (x x x, t)/∂t is smaller than 1 × 10 −4 , the evolution of the crack phase-field is sufficiently relaxed and no visible difference in the material response is present anymore (see Fig. 4(c)). A higher value of φ c,crit results in the retardation of the drop-off of the load (Fig. 4(d)), while the computational time for macroscopic crack formation increases significantly (i.e. approx. 2.5 times higher for φ c,crit = 0.95 than for φ c,crit = 0.9 in Fig. 4(e)). Therefore, we chose φ c,crit = 0.9 as a trade-off between sufficiently accurate results and computational costs for the further simulations. Due to the small plastic region and the high hardening modulus, the non-linear inelastic regime in all the load-displacement curves is small compared to the elastic regime. The result of a purely elastic crack propagation, in which the cohesion c = 100c 0 is chosen, is given in Fig. 4(f) as the brown line in the plot. As the cohesion decreases, the plastic strains at the crack tip increase, resulting in the lower elastic strain and driving force for the propagation. Therefore, a higher displacement has to be applied before further crack propagation, as reflected in the plot in Fig. 4(f). As the hardening modulus decreases, the amount of plastic strains increases due to the slower expansion of the flow cone of Drucker-Prager plasticity (discussed in Section 2.2). As a result, a higher displacement load is required for lower values of hardening modulus (see Fig. 4(g)). Figure 4(h) illustrates the influence of varying friction angle on the load-displacement response. As the friction angle decreases, the apex of the flow cone moves to the right on the positive hydrostatic axis (illustrated in Fig. 1(b)), resulting in smaller plastic strains compared with higher friction angles at the same stress state. Therefore, for smaller friction angles, the crack starts propagating at lower displacements. Furthermore, in the present diffuse interface framework, there is an influence of the crack interface width (related to c ), as depicted in Fig. 4(i). The sensitivity of the material response with the interface width is attributed to the evolution equation of the crack phase (see Eq. 4). In order to describe a smooth transition in the interface with a sufficient number of grid points (approx. 10 see e.g. [80]) and to reduce computational time, we chose c = 5Δx for the simulations in the forthcoming sections. The length-scale parameter c can be associated to a physical interface width and can be set as a material parameter, provided the availability of precise information about the physical material properties, real microstructures and experimental loading conditions (see e.g. [65]).

Crack propagation in a two-phase specimen: brittle to elasto-plastic transition
In this section, we investigate the fracturing behaviour in a two-phase specimen consisting of a purely elastic and an elasto-plastic phase. The simulation setup along with the initial crack is displayed in Fig. 5(a). The initial tip of the crack lies in the elastic phase. Moreover, the purely elastic solid is on the left (described by the phase-field φ 1 ), while the elasto-plastic material (with phase-field φ 2 ) obeying Drucker-Prager flow criterion is present on the right. A diffuse solid-solid interface is present between the two phases, where the diffuse region lies in the transition enclosed by the dashed lines in Fig. 5(a). The boundary conditions and the incremental displacement loading are the same as in the previous Section 4.2.1. Moreover, the material parameters in both phases are identical and their chosen values are also the same as in Section 4.2.1. In contrast to the cohesion of the ductile phase, c 2 = c 0 , a higher value of c 1 = 100c 0 is chosen for the purely elastic phase, so that no plastification can occur.
We analyse results at three different positions of the crack tip: In the bulk region of the brittle phase, Just after passing the brittle-ductile diffuse interface and In the bulk region of the ductile phase, as shown in Fig. 5(b). After an initial linear elastic behaviour, when the displacement is sufficiently high, crack propagation begins in the brittle phase (at position in Fig. 5(b)), resulting in a drop in the reaction force until the crack reaches the diffuse interface of the sample (see in Fig. 5(d)). Plastic strains evolve at the crack tip accompanied by a reduction of the strain energy density in the diffuse interface region. Consequently, further displacement increments are applied to propagate the crack. Due to the employed homogenisation approach, a continuous relic plastic zone is formed over the diffuse interface region between the brittle and ductile phases, as illustrated in the inset picture of Fig. 5(c). After passing the diffuse interface, the crack reaches the bulk region of the ductile phase (at position in Fig. 5(b)), and crack propagates without further displacement increments. The reaction force drops to 0 until complete fracture of the specimen. In order to analyse the stress distribution near the crack tip, we plot the normalised stress in the loading direction, i.e. σ yy /σ yy,max along the white dashed-dotted line at the three positions , and (see Fig. 5(e)). The chosen value of the scaling factor for the normalised stresses corresponds to the stress peak at position . Due to the incorporated linear hardening, the stress increases with the plastification of phase φ 2 . Furthermore, the stress peaks in the plastic regions no longer remain close to the crack tip, as illustrated in Fig. 5(e). Moreover, simulations are performed for different values of cohesion c 2 of the ductile phase (see Fig. 5(d) for the load-displacement responses). It is observed that all curves show the same behaviour in the brittle phase until the crack reaches the diffuse interface. In the diffuse interface, lower cohesion of the second phase results in a higher amount of plastic strains at the crack tip, therefore requiring higher displacement loads for crack propagation. Finally, in order to analyse the influence of the interface width of the solid-solid interface s , additional simulations are conducted for s = 3Δx and s = 8Δx. The curves (in Fig. 5(f)) reveal an insignificant impact of the solid-solid interface width on the resulting loaddisplacement responses. For the values of s = {3Δx, 8Δx}, a maximal deviation of about 1.3% is obtained that is arising due to the homogenisation approach with phase-inherent plastic strains.

Crack propagation in multiphase system: application to exemplary non-homogeneous geological structures
In the above numerical examples, we present a systematic study of modelling elasto-plastic crack propagation in single-and two-phase domains. Through the benchmark examples (in Section 4.1), we validate the implemented plasticity model in the diffuse interface framework (using the homogenisation approach) for binary and higher order interfaces, and show that a sound agreement is achieved when compared with the sharp interface results of the commercial FE solver (ABAQUS). Moreover, the simulations of crack propagation in a single-phase ductile material elucidate the influence of different material parameters on the load-displacement response, while successfully recreating the (a) active plastic zone near the crack tip and the (b) relic plastic zone around the crack boundaries (see e.g. [4]). The present approach further captures the brittle-to-ductile transition during fracturing of a two-phase heterogeneous specimen with a pre-existing crack and reveals that a continuous relic plastic zone can be expected for a smooth transition between brittle and ductile regions. All these examples advocate that the present approach is capable of computationally mimicking the elasto-plastic fracturing in multiphase systems.
In this section, we apply the model to multiphase geological structures consisting of spatially separated regions of brittle and ductile rocks, which are observed in the earth's crust. For the sake of demonstrating the numerical performance, we chose the elastic material parameters (i.e. Young's modulus E and Poisson's ratio ν) for the ductile phase corresponding to sandstone, while for the brittle phase, the values corresponding to quartz are chosen. The elastic material parameters for sandstone vary [81][82][83][84][85] depending upon various geophysical factors such as water content, average grain size and porosity. We consider the values for sandstone from [81] and [82], as listed in Table 1. The values of the plastic parameters, i.e. cohesion c and friction angle β, are taken from [81]. The elastic material parameters for the brittle phase are chosen from [86]. Furthermore, the brittle phase is assumed to exhibit isotropic material behaviour (i.e. isotropic stiffness), consistent with the implemented tension-compression split (see Eq. 21). As the elastic material parameters, the fracture toughness for sandstone is dependent upon the geophysical factors as well and is reported in the range of K Ic = 0.3-1.5 MPa √ m [83,85]. The energy release rate of an elasto-plastic fracture includes both the specific fracture energy and the plastic dissipation, which can be hardly separated in experimental measurements of fracture toughness. Therefore, the measured fracture toughness includes both elastic and dissipative effects [87]. In the present work, the dissipative effects are already considered in the strain energy density. Therefore, we assess a reduced G c compared with experimental results (G c = 1.5 J/m 2 ) for both sandstone and the brittle phases. Due to the incorporated isotropic Drucker-Prager plasticity, the fracture toughness for sandstone is embraced to be isotropic. However, in order to show the influence of the anisotropy (Eq. 10), the fracture toughness of the brittle phase is assumed to be anisotropic with a factor of F a,x = 0.7 and F a,y = F a,z = 1.
Crack propagation in a heterogeneous system is a complex phenomenon and is dependent on mechanical properties, material response and structure of the rocks. In a rock composed of brittle and ductile phases, the crack path is primarily decided by an interplay of the elastic strain energy distribution, the fracture resistance of the rock phases and their associated anisotropies (see [67]). In order to illustrate this interplay, we consider two different geological structures comprising of brittle and ductile phases (see Figs. 6(a) and 7(a)). The simulation domain size for both the setups is 151Δx × 151Δx, where Δx represents a physical length of 100 μm. The phase structures for the two simulation domains are constructed using a Voronoi algorithm. The blue and orange phases are assigned the brittle material parameters, where the ellipse determines the orientation of the fracture toughness. The yellow phase is assigned the ductile material properties of sandstone and isotropic fracture toughness as highlighted by the circle in Figs. 6(a) and 7(a). The ductile and brittle material parameters are given in Table 1. An initial crack is set into the blue brittle phase on the left-hand side of both the setups, and an incremental displacement loading is applied in the orthogonal direction of the upper domain edge. The maximal displacement for structure 1 is u max,1 = 1.8 μm and for structure 2 u max,2 = 2.06 μm. We analyse the results at two different stages of crack propagation for both the structures.
For structure 1, Fig. 6(b-d) depicts the crack phase, the accumulated plastic strain and the strain energy density at an intermediate stage after 2000 dimensionless time steps. The inset picture of Fig. 6(b) displays the contour plot of the strain energy density when the crack propagation is about to begin at position . Due to uniaxial loading, the elastic driving force is observed to be the highest in positive x-direction in the cell next to the crack tip. On the other hand, the crack resistance is lowest along the direction highlighted in the same inset picture. As a result of the interplay of the two factors, the initially horizontal crack takes an intermediate path towards the direction of the lowest crack resistance, as shown in Fig. 6(b). As soon as it reaches the interface of the brittle-ductile transition, the plastic strains begin to accumulate near the crack tip forming plastic zones, as depicted in the zoomed inset picture in Fig. 6(c). In the bulk region of isotropic ductile phase, the direction of crack path is determined solely by the elastic strain energy distribution, which is observed to be high in the region adjacent to the tip due orthogonal loading (see Fig. 6(d)). Therefore, the crack propagates in a straight path (in positive x-direction) perpendicular to the direction of the loading, leaving behind a continuous relic plastic zone in the ductile phase, as can be seen in the zoomed inset picture of Fig. 6(f). As the crack moves into the brittle material, it again aligns itself towards the direction of lowest crack resistance until complete fracturing (after 8800 time steps, in Fig. 6(e, f)). For structure 2, similar to the former case, the crack deflects in the brittle phase according to the previously discussed interplay (Fig. 7(b)) and enters the ductile phase, with plastification near the tip (Fig. 7(c)). However, although the crack resistance is isotropic in the ductile phase, the crack prefers to stay in this phase passing through the regions of higher strain energy density (illustrated in Fig. 7(d)) which is energetically more favourable, and thereby, following a curved path (see Fig. 7(e)). Towards the end, when the crack again enters the brittle phase, a slight deflection is observed in the direction of lowest crack resistance as expected. In the ductile regions, both fractured structures exhibit a relic plastic zone, as illustrated in Figs. 6(f) and 7(f). Fig. 6 a Simulation setup for structure 1 depicting a domain with brittle (orange and blue) and ductile (yellow) phases, and pre-existing crack phase along with the domain size and the direction of the incremental displacement loading. b The crack path in the brittle phase is determined by the interplay of the strain energy density and the direction of lowest crack resistance. The contour plots of c illustrate the accumulated plastic strains and d the strain energy density when the crack reaches the diffuse interface. In the interface regions, the plastification results in the formation of a plastic zone. e The crack is continuously deflected in the direction of the reduced crack resistance. f A continuous relic zone is left in the ductile regions after complete fracturing Fig. 7 a Simulation setup for structure 2 depicting a domain with brittle (orange and blue) and ductile (yellow) phases, and pre-existing crack phase along with the domain size and the direction of the incremental displacement loading. b The crack deflects itself in the direction decided by the interplay of the strain energy density and direction of the lowest crack resistance. The contour plots of c refer to the accumulated plastic strains and d the strain energy density when the crack reaches the diffuse interface. e The crack takes a curved path in order to stay in the ductile material with regions of high strain energy density. f A continuous relic zone is left in the ductile regions after complete fracturing

Conclusion and outlook
In this work, we present an approach for modelling crack propagation in heterogeneous rocks consisting of brittle and elasto-plastic regions. The impact of anisotropic fracture toughness of the brittle phase and the plastic flow of the ductile phase on the crack propagation paths is illustrated. At first, the implemented plasticity model (based on the Drucker-Prager yield criterion) and its extension to the multiphase systems (based on the homogenisation approach) are validated through the following numerical experiments: 1. The simulations of two-phase specimens loaded according to parallel and serial combinations are able to recover the values of the von Mises equivalent stresses and accumulated plastic strains in the bulk regions, regardless of the chosen interface width. Moreover, in the bulk regions, the diffuse interface results exhibit sound agreement when compared with the sharp interface solutions of the commercial solver ABAQUS. In the diffuse interface region, a smooth and continuous transition of both the fields is obtained in accordance with the employed interpolation function h α s . 2. In a four-phase specimen, similar to the above case, the values of the von Mises equivalent stress and accumulated plastic strains in the bulk regions are independent of the chosen interface width, and further match (within 0.03%) with those obtained from ABAQUS. Near the domain edges and the quadruple junction, a relatively higher deviation of 4.9% is observed, which is attributed to the different solution schemes or chosen residuals in both the approaches.
After validating the plasticity model and the homogenisation approach for multiphase systems, we investigate the formulation of crack evolution based on Griffith's criterion. The simulations of crack propagation in elasto-plastic specimens successfully recreated the active plastic zone near the crack tip and a relic plastic zone around the crack boundaries, consistent with experimental observations [4]. Furthermore, the following aspects about elasto-plastic crack propagation are illustrated: 1. A homogeneous ductile specimen obeying Drucker-Prager flow criterion is analysed with a pre-existing crack and undergoing incremental displacement-driven loading. The amount of plastic strain near the crack tip increases with decreasing cohesion c and hardening modulus H values and increasing friction angle β values. The variation of parameters yields that higher displacements are required for crack propagation. Moreover, due to the influence of c in the evolution equation of the crack phase-field, the load-displacement behaviour exhibits a dependency on the crack interface width. 2. Simulations of crack propagation in a two-phase specimen comprising of a purely elastic and an elastoplastic phase further suggest that a continuous relic zone is formed if the transition region between the brittle and ductile regions is smooth and continuous. Moreover, after the occurrence of plastification in the ductile regions, the earlier vicinal stress peak shifts away from the crack tip.
Finally, in order to demonstrate the capabilities of the presented modelling, crack propagation is simulated in exemplary multiphase geological structures. To evaluate the influence of anisotropy in fracture toughness on the crack path, the brittle phases are assigned with an anisotropy of G c . The simulation results of two different structures reveal the following: 1. In a brittle phase exhibiting anisotropy in the fracture toughness, the crack path is determined by an interplay of the strain energy density and the direction of the lowest crack resistance (see also [67]). 2. In a ductile phase exhibiting isotropic plasticity and fracture toughness, the crack path is solely determined by the strain energy density distribution.
The presented numerical model is capable of simulating crack propagation under a wide variety of boundary conditions and heterogeneities of multiphase systems, on the basis of plasticity, anisotropic fracture toughness and formulation for crack growth in homogeneous and heterogeneous systems with purely elastic and inelastic phases. Given the strongly non-linear nature of the multiphase-field model, due to the incorporation of the Drucker-Prager plasticity, brittle phase anisotropy and the mechanical jump conditions, a staggered scheme for the update of mechanical and phase fields is employed. The mechanical fields were updated implicitly, while the phase fields are updated in an explicit manner. For this coupled multiphase-field modelling of fracture, from the algorithmic point of view, a fully implicit solution scheme should be implemented, while assessing the convergence characteristics. With the aforementioned model ingredients (i.e. plasticity, anisotropic elasticity, jump conditions and coupling with fracture), the present work serves as an initial treatment and paves the way for more advanced and quantitative ductile fracture models for heterogeneous materials at large strains. The constraint on the choice of hardening modulus can be further relaxed by appropriately formulating the crack resistance and/or crack driving force, such that the accumulated plastic strains enter their formulations, i.e. G c (φ, ∇φ c ,ε pl ) and f α el (ε α el , φ c ,ε pl ) [47]. Another approach would be the extension of the implemented rate-independent plasticity model to a viscoplastic formulation [50]. Model extensions based on these topics are currently a work in progress, and will be addressed in the future. The presented modelling framework can be extended to include other plasticity models (e.g. Mohr-Coulomb). Although the present work is a two-dimensional treatment of the fracturing process, the extension to three dimensions is straightforward, although with higher computational costs. The present work, when coupled with fluid flow and temperature equations, would enable further investigations on the formation of fracture networks (by hydraulic pressure and mechanical deformation) and alternating crack-sealing processes in geothermal/hydrocarbon reservoirs, along with their impact on the flow behaviour of rocks. When coupled with the crystal growth model of Ankit et al. [88], this work can also be employed to investigate different vein formation mechanisms.