Brittle anisotropic fracture propagation in quartz sandstone: insights from phase-field simulations

We developed a generalized multiphase-field modeling framework for addressing the problem of brittle fracture propagation in quartz sandstones at microscopic length scale. Within this numerical approach, the grain boundaries and crack surfaces are modeled as diffuse interfaces. The two novel aspects of the model are the formulations of (I) anisotropic crack resistance in order to account for preferential cleavage planes within each randomly oriented quartz grain and (II) reduced interfacial crack resistance for incorporating lower fracture toughness along the grain boundaries that might result in intergranular crack propagation. The presented model is capable of simulating the competition between inter- and transgranular modes of fracturing based on the nature of grain boundaries, while exhibiting preferred fracturing directions within each grain. In the full parameter space, the model can serve as a powerful tool to investigate the complicated fracturing processes in heterogeneous polycrystalline rocks comprising of grains of distinct elastic properties, cleavage planes, and grain boundary attributes. We demonstrate the performance of the model through the representative numerical examples.


Introduction
Fractures are an integral part of reservoir rocks [1] and are formed when the stresses exceed the strength of the rock material. These stresses can occur due to a wide variety of mechanical conditions causing the cracks to propagate isotropically or anisotropically depending upon the physical properties of rocks. Sandstones, at microscopic length scale, are predominantly composed of quartz grains of different shapes and sizes that are cemented together as Nishant Prajapati, Christoph Herrmann, and Michael Späth have contributed equally to this work. Nishant Prajapati nishant.prajapati@kit.edu a result of precipitation of minerals (e.g., silica SiO 2 and calcite CaCO 3 ) in the pore spaces [2][3][4]. Each quartz grain exhibits a crystallographic orientation which is different from its neighboring grains, thereby resulting in the presence of grain boundaries. The material properties vary depending upon the crystallographic directions. For grains of prismatic habit (with (1011), (0111), and (1010) facets, see Fig. 1a), experiments [5][6][7][8][9] report that the fracture toughness (and therefore, the crack resistance) varies along different crystallographic directions, leading to preferential planes for intragranular crack propagation. Furthermore, when the crack resistance along the grain boundaries is lower as compared to grains, intergranular fracturing might occur. Depending upon the geophysical conditions in the source area, rock fragments and grains of other minerals (e.g., k-feldspar) also contribute to the bulk of sandstone [10]. The crack paths are also influenced by the presence of secondary minerals with different material properties. All these grain-scale heterogeneities and anisotropies, in addition to the mechanical loading conditions, account for the fracture path, imparting complexities to the process of crack propagation in sandstones. Fractures that originate at microscale not only serve as nucleation sites for  [11] but they also control the mechanical strength of rocks and provide permeability for the flow of fluids and hydrocarbons. Therefore, a deep understanding of fracturing process is crucial for understanding the stress states and petrophysical properties in reservoir rocks.
Computational modeling and simulations are powerful tools that assist in the investigation and characterization of fractured hydrocarbon and geothermal reservoirs. They are employed for predictive fracture-induced failure analysis in a wide variety of other rock engineering applications such as rock cutting and underground excavations. Different numerical approaches have been proposed to address the problem of fracturing in geological systems. Zhang and Jeffrey [12] utilized the boundary element method (BEM) for analyzing the formation of fracture networks through fluid-driven crack growth in 2-D. However, as only boundary is discretized in this approach, the method faces difficulties in addressing heterogeneities and anisotropies of the rocks. Using the extended finite element method (XFEM), Wang et al. [13] simulated the interactions between hydraulic and natural fractures, and the formation of fracture networks in 2-D. Moreover, 3-D investigations were conducted by Virgo et al. [14][15][16] using the discrete element method (DEM) for analyzing the interactions of veins and fractures. DEM-based models treat material as accumulation of spherical particles and track each individual particle and its interactions with the neighboring particles over time. As an increase in the number of particles n leads to higher computational costs (typically scaling as order of n), these models are not well-suited for large-scale numerical simulations. For a comprehensive understanding of the computational methods in fracture mechanics for rocks, the interested readers are referred to the recent review article [17]. In the abovementioned numerical approaches, cracks are treated as sharp material discontinuities, placing these models in the category of so-called sharp interface models. Development of robust algorithms for the update of crack surfaces and their numerical implementation based on sharp interface approach is a tedious task, especially when dealing with the problem in 3-D.
In contrast to sharp interface methods, diffuse interface approaches such as the phase field method describe crack surfaces as diffuse regions and do not require explicit tracking of interfaces. Due to this reason, the method has emerged as a robust and computationally efficient alternative and has been extensively employed for modeling both brittle [18] and ductile [19][20][21][22] fracture propagation in homogeneous materials. In geological materials, approaches like coupling phase-fields models with mechanics have been developed for addressing different problems (e.g., cracking from crystallization in pores [23], fracture with pressure sensitive plasticity in geomaterials [24], mixed mode brittle fracturing in anisotropic brittle rocks [25], crystal plasticity formulation coupled with anisotropic phase-field fracture for crack propagation in single-crystal halite (rock salt) [26]). Furthermore, recent works [27][28][29][30][31][32] are coupled geomechanics based on phase-field fracture with the flow equations for simulating hydraulic fracturing in porous media. While the aforementioned works [23][24][25] focused on fracturing in homogeneous geomaterials, Na and Sun [26] simulated the elastoplastic crack propagation in an anisotropic halite bicrystal. In heterogeneous materials, comprising of multiple phases, Wang et al. [33] published one of the earliest works to address the problem of crack propagation with application to polycrystalline systems. However, the problem of intergranular and transgranular fracturing using the phase-field method, that largely remained unaddressed, was investigated for polycrystalline materials by Abdollahi and Arias [34] and Oshima et al. [35]. However, they considered isotropic fracture toughness in the bulk of the crystals, and therefore, are not appropriate to describe fracturing in polycrystalline materials exhibiting preferential cleavage planes within each randomly oriented crystal. This aspect was addressed by Liu and Juhre [36] based on the anisotropy formulation considered in Clayton and Knap [37]. Nguyen et al. [38] extended this model to further account for the inter-and transgranular cracking based on the cohesive zone model (CZM) for describing the decohesion of the grain boundaries. However, one limitation of CZM is that the cohesion surfaces can only lie along the finite element edges. Therefore, the crack paths become mesh-dependent. Furthermore, in the anisotropy formulation presented in the works [25,26,36,38], the orientation dependency is only incorporated in that term of the free energy functional which contains the gradient of the phase fields, causing different interface widths of the crack phase-field propagating in different crystals. Variable crack interfacial widths in different solid phases lead to different crack surface energies, resulting in energetic deviations from the sharp interface results.
In the present work, we developed a new and robust multiphase-field crack model that accounts for the varying fracture toughnesses using the interpolation functions of the phase fields and ensures constant crack interface widths in different crystals. Motivated from the modeling framework of Schneider et al. [39], the present work showcases novel formulations for incorporating crack resistance anisotropy and lower crack resistances along the solid-solid interfaces, well-suited for describing brittle crack propagation in sandstones that exhibit crystal symmetry and crack resistance anisotropy at grain scale. These formulations enable the model to simulate trans-and intergranular fracture propagation, while exhibiting preferential crack growth depending upon the weakest planes within each grain.
The present article is organized as follows: Section 2 elaborates the multiphase-field model for crack propagation, formulations for anisotropic, and reduced interfacial crack resistances. In Section 3, we demonstrate the capabilities of the model through the representative numerical examples that include the applications to an exemplary geological vein system. Finally, we conclude the article by recapitulating the main inferences drawn in the numerical investigation and directions for further work in Section 4.

Multiphase-field model for crack propagation
We consider a domain Ω ⊂ R d , d ∈ {1, 2, 3} consisting of N ∈ N solid phases (e.g., quartz and secondary minerals) and a crack phase c. The presence of each phase α ∈ {1, . . . , N; c} is determined by an order parameter φ α (x, t) : Ω × R + 0 → [0, 1] at spatial position x and time t. The region of domain which is solely occupied by phase α is known as α-bulk, and is mathematically written as B α = {x ∈ Ω | φ α (x, t) = 1}. The total bulk region in the domain Ω is given by the union of bulk regions of individual phases, i.e., B = ∪ α B α . The grain boundaries (i.e., solidsolid interfaces) and the fracture surfaces (i.e., crack solid interfaces) are described as diffuse regions, collectively called the total interface region I Ω = Ω\B. A binary interface between any two phases α and β is written as where the value of the order parameter φ α varies from 1 in α-bulk to 0 in β-bulk smoothly and monotonically over the diffuse interface, and vice versa for φ β . The location of interfaces, both solid-solid and crack solid, are determined by the isosurface φ α = 0.5. At any material point x in the domain Ω, the extent of damage of the solid phase is determined by the crack phase field φ c , with φ c = 0 and φ c = 1 describing the fully intact and the fully destroyed material states, respectively. The interface region shared by a solid and crack phase where φ c ∈ (0, 1) represents the partially damaged material states. For the sake of convenience, all the order parameters are collectively represented as a set of phase fields given by , t), . . . , φ N (x, t)} represents the phase fields of all the solid phases. In a system entirely composed of for example N quartz grains, multiple solid phase fields φ α ∀ α ∈ {1, . . . , N} are required in the subset φ φ φ s (x, t) to assign distinct crystallographic orientation (and thereby direction dependent material properties) to different grains. At every material point x, the summation constraint N α=1 φ α = 1 − φ c is fulfilled. According to Griffith's theory of brittle fracture [40], crack propagation takes place when the released elastic strain energy at a material point exceeds the surface energy required to create free crack surface. Therefore, in order to model fracturing process based on the thermodynamic Griffith's criterion, we formulate the total free energy functional of the system of volume Ω as follows: In the integrand on the RHS of Eq. 1, G c (φ φ φ s , ∇ ∇ ∇φ c ) denotes the effective crack resistance. The term in curly brackets, that comprises of the sum of gradient c |∇ ∇ ∇φ c | 2 and potential energy density contribution w c (φ c )/ c , maps the energy of a propagating crack in a regularized manner. The scalar parameter c controls the crack interface width. For the present work, we chose a one side well-type potential w c (φ c ) = k w φ c 2 for the crack phase energy density contribution. The values of constants k = 0.5 and k w = 1 are computed by integrating the surface energy of the diffuse crack over the whole interface and equating it to the sharp interface solution [39]. The effective crack resistance G c (φ φ φ s , ∇ ∇ ∇φ c ) is given by the volumetric interpolation of the crack resistance of the individual solid phases G α c (∇ ∇ ∇φ c ), and reads as follows: Following the work of Moelans [41], a normalized interpolation function for the solid phases h α s (φ φ φ s ) of the form is chosen in order to satisfy the local constraint N α=1 h α s (φ φ φ s ) = 1. The treatment of crack resistance for anisotropic crack growth and formulation of weak grain boundaries for intergranular fracturing are discussed in Sections 2.2 and 2.3, respectively. The term f el (φ φ φ, ε ε ε) in Eq. 1 denotes the effective elastic strain energy density as a function of phase fields φ φ φ and the total strains ε ε ε. Under the assumption of small deformations, the total strain at any material point is given by the following: in terms of the displacement field u(x, t). Similar to effective crack resistance, the effective elastic strain energy density is also given by the interpolation of the elastic strain energy densities of the individual solid phases f α el (φ φ φ, ε ε ε α ), and reads as follows: The evolution of crack phase field is given by the following: where μ is a positive scalar relaxation parameter that controls the crack propagation velocity. When accounting for the dynamic effects in the momentum balance equation, the crack velocity can be controlled by appropriately formulating the term μ. Under the quasi-static assumption of the present work, the term μ serves as a positive constant relaxation parameter, and is chosen sufficiently large in magnitude such that the simulations are numerically stable. On substitution of the free energy from Eq. 1 to 6, we get the following: The term ∂f el (φ φ φ, ε ε ε)/∂φ c accounts for the driving force for crack propagation. In order to prevent fracture growth under compressive states that may lead to physically unrealistic crack patterns, the elastic strain energy density of each phase is split into positive and negative parts following the work of Miehe et al. [42]. Example of unrealistic crack patterns in compression has been reported in, e.g., Bourdin et al. [43]. For a detailed review of the different tension-compression split formulations (e.g., [42,44]) and the comparative analysis of their performance, interested readers are referred to Ambati et al. [45]. For the present work, the two parts of the phase-dependent elastic strain energy density are defined as follows [42]: where C C C α = λ α 1 1 1 ⊗ 1 1 1 + 2μ α I is the phase-dependent isotropic stiffness tensor, λ α and μ α are phase-inherent Lame parameters, and α i for i = 1, 2, 3 are the phasedependent principal strains. The operator ± when applied over any scalar quantity q is evaluated as q ± = (q±|q|)/2. Therefore, the driving force using the tension-compression split formulation of Eq. 8 is computed using only the positive part of the strain energy density. In the absence of solid-solid transitions, the evolution of solid phase fields is obtained as follows: The balance of linear momentum under the quasi-static assumption (i.e., neglecting inertia and body forces) given by the following: is solved for the displacement field u. The effective stresses σ σ σ are given by the following: as an interpolation of the phase-dependent stresses σ σ σ α which are computed as follows: Unlike the crack driving force, the phase-dependent stresses are computed by the degradation of the total strain energy density. This retains the linearity of the momentum balance equation (Eq. 11). Therefore, the present formulation is rendered as a hybrid model after Ambati et al. [45]. Within the staggered algorithmic implementation (as also utilized in the present work), and from the computational point of view, the hybrid models have shown to lead to similar results (in terms of fracture patterns as well as load-displacement behavior) compared to those where tension-compression split enters the calculations of the crack phase field as well as the stresses, but at much lower computational costs than the latter (about one order lower in magnitude [45]). Owing to these reasons, we chose this hybrid formulation for our present work. The phase-inherent stresses in the interface region are determined based on the mechanical jump conditions discussed in Schneider et al. [46] and Herrmann et al. [47].
The balance of momentum (Eq. 11) and the phasefield evolution equations (Eqs. 7 and 10) are solved in a staggered manner. The balance of momentum is solved implicitly for the mechanical fields. While, the evolution equation for the crack phase (Eq. 7) is solved using explicit Euler scheme for time derivative and central differences for the spatial derivatives. A rotated staggered grid is employed for the discretization, where the phase fields and stresses are evaluated at the center of each cell, while the displacements are stored at the cell corners. Therefore, the present grid is analogous to a regular finite element (FE) mesh with equal-sized linear elements and full integration. For each time step, the mechanical and phase fields serve as input for each other. For each time increment Δt, the conditionφ c ≥ 0 is enforced. Furthermore, for the purpose of reducing the computational costs, a critical value of φ critical c = 0.9 was set, above which a material point is considered fully cracked (i.e., φ c = 1.0). This numerical treatment ensures that the values of the crack phase field remains in the range φ c ∈ [0, 1]. It is noteworthy that the incremental time update of the solid phase fields (in Eq. 10) is performed as a volumetric interpolation of the incremental time update of the crack phase field. This essentially implies that, if the computational domain is initialized such that it satisfies the summation constraint ( N α=1 φ α = 1−φ c ) locally, the constraint will be sustained at each successive time step. Within each mechanical load increment, the balance of linear momentum is solved for each relaxation time step for the evolution of crack phase field (Eq. 7) until it reaches a steady state. Thus, the next increment is applied as soon as ∂φ c /∂t < 0.0001 is fulfilled. This numerical treatment leads to negligible influence of the time step Δt and the scalar relaxation parameter μ on the simulation results. The numerical and model parameters are listed in Table 1. The model is implemented in a multiphysics simulation framework called PACE3D [48], which is written in C programming language. The solver is highly parallelized based on the message passing interface (MPI) standard that includes domain decomposition and dynamic redistribution schemes, facilitating large-scale simulations for resolving, e.g., microstructures with large number of grains. For the present work, simulations were performed on linux high performance cluster using multiple central processing units (CPUs).

Anisotropic crack resistance formulation: preferrential growth directions
For quartz grains exhibiting prismatic habit (see Fig. 1a), experiments report that the fracture toughness varies along different crystallographic directions [5][6][7][8][9]. The fracture toughness values along the cand a-axis (see Table 2) have been reported in the works of Norton & Atkinson [7] and Ferguson et al. [8], respectively. Based on the available data, we assume an anisotropy of the general form in order to describe the spatial variation of crack resistance within each quartz grain. Here, G α c,0 denotes the isotropic crack resistance of the solid phase α, and n α c,x (∇ ∇ ∇φ c ), n α c,y (∇ ∇ ∇φ c ), and n α c,z (∇ ∇ ∇φ c ) are the x-, yand z-components, respectively, of the rotated unit vector n α c (∇ ∇ ∇φ c ) normal to the crack interface (shared with the solid phase α), given by the following: The phase-dependent rotation matrix Q α describes the orientation of the solid phase anisotropy, in accordance with the crystallographic orientation of each quartz grain. The phase-dependent anisotropy strength factors f α d , ∈ [0, 1] for d = x, y, and z reduce the crack resistance in the x-, y-, and z-directions, respectively, of the rotated coordinate system. As the considered 3-D growth habit of quartz (Fig. 1a) exhibits three a-axes lying in the plane perpendicular to the c-axis, we assume an isotropic crack resistance in this plane. In order to model this transversely anisotropic behavior, we chose f α x = f α y = f α ∈ (0, 1) and f α z = 1. Such a selection isotropically reduces the crack resistance along the x-y-plane (or x-direction in case of 2-D) that corresponds to the plane containing the a-axes (or a-axis direction in case of 2-D) of the randomly oriented quartz grain. The anisotropy strength factor f α represents the ratio of the crack resistances along the a-and c-axis. Therefore, if the crack resistance along c-axis is G α,c-axis c = G α c,0 , then the corresponding value along the plane containing a-axes is given by G α,a-axis c = G α c,0 f α , giving rise to fracture propagation along preferential plane.  [38] and Liu and Juhre [36], in the present formulation, the G c (φ φ φ s , ∇ ∇ ∇φ c ) term with the orientation dependency is multiplied to both potential and gradient terms (Eq. 1) of the diffuse crack surface energy, leading to invariant interface widths in different anisotropic solid phases. In the case of well-type potential where the diffuse interface is present in the whole computational domain, if the interface width is defined by the region lying between φ c = 0.15 and φ c = 0.85, the present formulation ensures constant crack interface widths.

Reduced interfacial crack resistance formulation: weak grain boundaries
Intergranular fracture propagation might occur in sandstones if the crack resistance along the grain boundaries is lower than that of the grains. In a binary interface shared between two solid phases α and β (Fig. 2a), the crack resistance according to the normal interpolation using Eq. 2 reads as follows: where G α c (∇ ∇ ∇φ c ) and G β c (∇ ∇ ∇φ c ) are the crack resistances of the two phases. A sinusoidal profile was employed for the spatial interpolation of the solid phase fields in the interface region (see Fig. 2b). The reduced interfacial crack resistance is formulated as follows:  which is a third-order polynomial function of φ α or φ β . Figure 2c depicts the plot of normalized crack resistance G red c (φ α , φ β , ∇ ∇ ∇φ c )/G α c (∇ ∇ ∇φ c ) along the lines I-II (of Fig. 2a) according to Eq. 17 for different values of the reduction factor ζ , when G β c /G α c = 1. In case of quartz, where the maximum difference in the crack resistance between the two grains can be of a factor of 0.795 (as discussed in Section 2.2), the reduced interfacial crack resistance interpolation for different values of reduction factor ζ is shown in Fig. 2d. The variation of normally interpolated crack resistance according to Eq. 16 is also depicted in black color in Fig. 2d. The diffuse interface description of the grain boundaries enables the reduction of the crack resistance in a continuous and smooth manner, highlighting the utility of multiphase-field approach (with different phase fields for each grain) in the present case.

Representative numerical examples
Quartzarenite sandstone, based on the grain composition, comprises primarily of quartz grains along with a small amount (usually less than 5%) of secondary minerals (kfeldspar). Crack propagation in such rocks is controlled by the material properties (i.e., Young's modulus E, Poisson's ratio ν, and fracture toughness K I c ) of the grains and the grain boundaries. As discussed previously, quartz grains have been found to exhibit different fracture toughness along different crystallographic directions that gives rise to preferential direction for intragranular crack propagation. Furthermore, as the fracture tip reaches a grain boundary (quartz-quartz or quartz secondary mineral), depending upon the local stress state and the crack resistance in the vicinity of the tip, three different modes of crack propagation might occur, which are as follows: 1. Intergranular mode: Propagation along grain boundaries 2. Transgranular mode: Propagation through the grains 3.

Mixed mode: Combination of inter-and transgranular modes
In the upcoming sections, we demonstrate the capabilities of the present multiphase-field model in simulating crack propagation under different modes in quartzarenite sandstone through the representative numerical examples. The material parameters of quartz used for performing the simulations are given in Table 2, unless otherwise mentioned. We remark that although quartz is known to exhibit elastic anisotropy [50,51], we utilized the values of the quasiisotropic elastic constants for quartz (Table 2) calculated by Heyliger et al. [50], in order to account for the tensioncompression split [42] that was proposed for materials exhibiting isotropic elastic properties. For determining the phase-dependent isotropic stiffness tensor C C C α , the values of Lame parameters λ α and μ α for each phase were computed The present analysis is restricted to 2-D, to reduce the complexity of geometrical parameters, and to limit the computational effort. All the simulations were performed under the plane strain assumption. Thus, the crack resistance G c (used in the model formulation) is calculated as G c = K 2 I c (1 − ν 2 )/E. In the numerical setup of all the examples, the phase fields (both solid and crack) were initialized to ensure that the summation constraint N α=1 φ α = 1 − φ c is locally sustained at each spatial point x.

Two-phase specimen: isotropic crack resistance
We consider a two-phase specimen with a preexisiting fracture under tensile load (see Fig. 3a). The fracture is represented by a diffuse crack phase field φ c as depicted in Fig. 3b. For the solid phases α and β, equal values of the elastic constants (i.e., E and ν) are chosen which are given in Table 2. The crack resistance of phase α is chosen to be G α c = (K c-axis In all the cases, the crack propagates along a straight line perpendicular to the direction of displacement increments (see Fig. 4a). Figure 4b-f depict the displacement, crack driving force, and stress fields for the case of λ = 3 when the crack tip is at position (III). Due to isotropic crack resistance, all the fields exhibit symmetricity with respect to the white dash-dotted line passing through the crack path. In all the simulations, at any stage of crack growth, the following inferences about the displacement and stress fields are drawn: 1. The displacement field u 22 (Fig. 4b) on the right side of the crack tip (i.e., unfractured side) varies smoothly in the direction of loading, while the fractured side on the left exhibits sharp displacement transition between the separated parts as expected. 2. Symmetric crack driving force field next to crack tip (Fig. 4c), resulting in a horizontal crack path. 3. The normal stress components σ 11 and σ 22 (Fig. 4d and e) exhibit peak values next to the crack tip. 4. The magnitude of shear stress component σ 12 is maximum right below and above the crack tip (see Fig. 4f).
When the crack tip is at its initial position ((I) in Fig. 4a) and the crack is about to propagate, we analyze the variation of normal components of stress (i.e., σ 11 and σ 22 ) along the white dash-dotted line (Fig. 4a) and shear component (i.e., σ 12 ) along the green dash-dotted line (in Fig. 4a), as shown in Fig. 4g-i. It is observed that the position and magnitude of the peaks of all components remain invariant for different values of λ for a given G α c . However, values of the normal stress components in phase β increase as the crack resistance increases. The variation of shear component along the green dash-dotted line exhibits negligible dependency on the crack resistance of phase β. We plot the load-displacement response for different cases (see Fig. 4j). It is inferred that: 1. For λ ≤ 1 (yellow regime): The load monotonically increases with displacement, reaches a maximum ((I) in Fig. 4j), while the crack tip is still at its initial position ((I) in Fig. 4a) just before the start of crack propagation. As the value of λ decreases, a softening behavior is observed which is attributed to the crack resistancedependent normal stress components (in Fig. 4g and h). As the crack begins to propagate, the load drops to zero without further displacement increments. 2. For λ > 1 (gray regime): Similar to the above case, a monotonic load-displacement behavior is observed until maximum (I), after which the crack starts to propagate in phase α, resulting in lowering of the load until the crack tip reaches the interface. At this position, due to higher crack resistance of phase β, the load begins to increase monotonically with displacement Fig. 4 a Crack propagates along a straight path from phase α to phase β with isotropic crack resistances G α c and G β c , respectively. Contour plots of b displacement u 22 in y-direction, c driving force field and stress components, d normal stress in y-direction σ 22 , e normal stress in x-direction σ 11 , and f tangential stress σ 12 , when the crack tip is at position (III), for the crack resistance ratio λ = 3. Plots of normalized stress components g σ 22 /σ 22,max , h σ 11 /σ 11,max along the white dash-dotted line and i σ 12 /σ 12,max along the green dash-dotted line (highlighted in Fig. 4a)  and reaches a second maximum ((II) in Fig. 4j). As soon as the crack growth starts in phase β, the load drops sharply to zero without additional displacement increment. The magnitude of load at second peak (II) is directly proportional to the value of λ.

Single-phase specimen: anisotropic crack resistance
We demonstrate the performance of the crack resistance anisotropy formulation (presented in Section 2.2) which allows the fracture to propagate along preferred planes within a quartz grain, and discuss the influence of the phaseinherent anisotropy strength parameter f α on the crack path. The numerical setup of a single-phase specimen with a preexisting crack is shown in Fig. 5a and b. The elastic constants (E and ν) for the specimen correspond to quartz and are given in Table 2. The crystallographic orientation of the grain with the a-axis inclination of − 45 • from the horizontal is indicated by the arrow in Fig. 5a. The directions parallel to the a-axis possess the lowest crack resistance. The value of crack resistance along the c-axis is G c-axis c = 57.71 × 10 −6 J/mm 2 . Figure 5c shows the crack growth along the preferred plane at an intermediate stage for For different values of the anisotropy strength parameter f α , simulations were performed while keeping other parameters identical. For each simulation, 17 CPUs on the high-perfomance cluster were employed. Figure 6 depicts the driving force field just before the outset of crack propagation (first row). The asymmetricity of the driving force field increases as the anisotropy strength increases (for decreasing values of f α ) resulting in higher deflection of the crack towards the direction of lowest crack resistance, illustrated by the simulated crack paths in second row of Fig. 6. Third row of Fig. 6 Based on the simulated crack propagation for quartz (Fig. 6b), crack deflections within each quartz grain in sandstone are expected to be small.

Two-phase specimen: reduced interfacial crack resistance
Intergranular fracture propagation might occur in sandstones if the crack resistance along the grain boundaries is lower compared to that of quartz grains. Figure 7a depicts the numerical setup of a two-phase specimen with a preexisting crack and two phases α and β possessing equal isotropic crack resistance G α c = G β c =Ḡ c = 57.71 × 10 −6 J/mm 2 in the bulk regions. The crack resistance along the grain boundary was lowered by incorporating the reduced interfacial crack resistance formulation (discussed in Section 2.3) using the reduction factor ζ , schematically illustrated by the crack resistance field in Fig. 7b and c. The crack is described by the crack phase field in Fig. 7b. The values of elastic constants (E and ν) are taken from Table 2. The reduced crack resistance in the interface region leads to asymmetricity in the crack driving force field pointing along the grain boundary. As the ζ decreases, this asymmetricity increases, resulting in different modes of crack propagation. For ζ = 0.9 (Fig. 8a), no visible deflection is observed in the fracture path as the crack driving force field is nearly symmetric, resulting in transgranular mode of fracturing. As ζ decreases, a transition to mixed mode occurs for ζ = 0.5 and 0.25 ( Fig. 8b and c) due to the increasing asymmetricity in the driving force field, leading to partial deflection of crack along the grain boundary. For ζ = 0.15 (Fig. 8d), the crack deflects completely along the grain boundary exhibiting intergranular fracture propagation. It is noteworthy that when the crack is passing through the grain boundary, the crack interface width (calculated between the φ c = 0.15 and φ c = 0.85) of 9-10 cells is obtained in all the cases. Each simulation in this section was performed using 17 CPUs of the high-performance cluster.

Crack propagation in exemplary multiphase geological structures
We consider a multiphase system analogous to a geological vein structure with 10 grains (data set of digital morphology available from Prajapati et al. [52]) with the material properties of quartz given in Table 2. In order to demonstrate the performance of the discussed model formulations in Sections 2.2 and 2.3, the following four cases for the crack resistance field in the computational domain were considered: -Case I: Equal and isotropic crack resistance in all the grains and along the grain boundaries (see Fig. 9a). The The boundary conditions for all the cases are identical and same as those of previous simulations, but a larger computational domain (238 × 238 × 1 cells) is considered. An initial crack phase is present in all the structures, shown in gray scale color map in Fig. 9. Figure 10 depicts different stages of crack propagation for the four cases along with the crack driving force field next to crack tip in the outlet pictures. For case I with homogeneous material properties and crack resistance, a straight crack path orthogonal to the direction of displacement loading is observed (Fig. 10a), due to symmetric crack driving force field near the crack tip at all the stages. For case II with lower crack resistance along the grain boundaries, the driving force field becomes asymmetric in the interface region,  leading to partial as well as complete deflections of crack in the grain boundaries at different points (see Fig. 10b). For case III with anisotropic crack resistance (Fig. 10c), the crack path is deflected according to the crystallographic orientation of each grain in the direction of lowest crack resistance. Due to a weak anisotropy strength, crack deflections in different grains are small (discussed in Section 3.2) without any influence of grain boundaries. It is worth mentioning that the crack interface width (between the isolines φ c = 0.15 and φ c = 0.85) is nearly same (16-18 cells) in different grains, as highlighted in the outset picture in the last stage of Fig. 10c. For case IV with anisotropic crack resistance of the bulk phases and reduced crack resistance values (ζ = 0.2) in the grain boundaries, the crack driving force field and thereby the crack path is governed by the interplay of both the factors.

Concluding remarks
The present work showcases a rigorous and generalized multiphase-field modeling framework developed for addressing the problem of crack propagation in multigrain systems at the grain scale, based on the thermodynamic Griffith's criterion [40]. With the novel formulations of anisotropic and reduced interfacial crack resistance, the model is well-suited for describing inter-and transgranular fracturing in brittle anisotropic sedimentary rocks such as sandstone. The introduced anisotropy in the crack resistance ensures that the crack interface width remains nearly constant in different grains. Utilizing the full space of the anisotropy strength parameters in the anisotropic crack resistance formulation, more complex material behavior based upon the fracture toughness can be modeled. The form of interfacial crack resistance formulation provides a general framework to assign lower or higher (in case of stiffer grain boundaries) values of crack resistance along the grain boundary varying smoothly between bulk-interfacebulk regions. For the purpose of demonstrating the model performance, a uniaxial loading along the direction perpendicular to the initial crack was applied in all the presented numerical examples. From the simulation results, the following conclusions about the numerical modeling and the process of fracture propagation can be drawn: 1. Crack path is controlled by the crack driving force field near the crack tip which is interrelated to the stress state and the crack resistance at a material point due to the coupling of set of equations for phase-field evolution and mechanical equilibrium. 2. For a material exhibiting anisotropic crack resistance, this crack driving force field becomes asymmetric exhibiting a tendency to point in the direction of lowest crack resistance and providing preferential growth directions for crack propagation. 3. When the crack resistance along a grain boundary is lower than that of the bulk, the driving force field points along the grain boundary due to lower crack resistance, facilitating intergranular fracturing.
Simulation results obtained by varying the values of different parameters show consistent fracture patterns, elucidated as follows: 1. As the strength of the crack resistance anisotropy in a given direction increases, the magnitude of deflection of the fracture in the direction of lowest crack resistance also increases, as expected. 2. As the interfacial crack resistance of the grain boundary reduces, the tendency of the crack to propagate along the grain boundaries increases.
In the present work, the values of the material parameters corresponds to quartz and are taken from literature [7,8,50]. However, we remark that the model parameters are chosen keeping in mind the numerical stability, computational costs, and the total time required for a simulation, and may not necessarily correspond to real geological systems. In the future works, the model will be extended to account for crack propagation under different loads (e.g., modes I, II, III) and more complex mechanical loading conditions (e.g., hydraulic fracturing due to fluid pressurization). In the loading scenarios where large parts of the solid are expected to exhibit compressive states, appropriate prescription for the condition of crack closure needs to be incorporated. Provided the exact material properties at grain scale, digital rock data, mechanical loading conditions, and the resulting fracture patterns in experiments, the model parameters can be calibrated to accurately mimic fracture growth patterns in polycrystalline rocks, and enhance the understanding of the formation of fracture networks. Although the presented set of simulations is limited to 2-D, the extension of the model to 3-D is straightforward, albeit with additional computational costs. The presented modeling framework, in its current form, can be applied to different polycrystalline systems consisting of multiple phases of varying elastic properties, anisotropies in stiffness as well as fracture toughness of the bulk phases and the interfacial crack resistance along the grain boundaries. Although the present work illustrates the application of the model to address the problem of brittle microfracturing in sandstones, the method is also applicable to multiphase/multigrain systems at larger scale. At macroscale, a coupling with macroscopic plasticity models (e.g., von Mises and Drucker-Prager) can be considered to describe plastic flow in certain inelastic phases of a heterogenous system. Therefore, the present work paves the way for future investigations pertaining to elastoplastic fracture propagation in multiphase systems.
Funding information Open Access funding provided by Projekt DEAL. This study is financially supported by the Helmholtz Association through the program "Renewable Energies (RE)," efficient use of geothermal energy 35.14.01. Furthermore, contributions of modeling aspects have been incorporated from the progress within the projects NE822/34-1 funded by the German Research Foundation (DFG) and "Kooperatives Promotionskolleg" funded by the Ministry of Baden-Wuerttemberg.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.