Phase-field modeling of hydraulic fracture network propagation in poroelastic rocks

Modeling of hydraulic fracturing processes is of great importance in computational geosciences. In this paper, a phase-field model is developed and applied for investigating the hydraulic fracturing propagation in saturated poroelastic rocks with pre-existing fractures. The phase-field model replaces discrete, discontinuous fractures by continuous diffused damage field, and thus is capable of simulating complex cracking phenomena such as crack branching and coalescence. Specifically, hydraulic fracturing propagation in a rock sample of a single pre-existing natural fracture or natural fracture networks is simulated using the proposed model. It is shown that distance between fractures plays a significant role in the determination of propagation direction of hydraulic fracture. While the rock permeability has a limited influence on the final crack topology induced by hydraulic fracturing, it considerably impacts the distribution of the fluid pressure in rocks. The propagation of hydraulic fractures driven by the injected fluid increases the connectivity of the natural fracture networks, which consequently enhances the effective permeability of the rocks.


Introduction
Hydraulic fracturing is a process of fracture of rock formation induced by pressurized liquids. From the last decade, hydraulic fracturing has been attracting more and more attention because of its applications in the extraction of oil and gas from unconventional reservoirs [1] and the enhancement of geothermal energy systems [2].
In a process of hydraulic fracturing, high-pressurized fluids are injected into rock layers in order to widen the existing fractures and create new artificial fractures. The resulting fracture networks consequently provide high-permeability pathways for the recovery of oil, gas, and geothermal energy. Due to their complexity, analytical solutions to hydraulic fracturing problems are very rare, particularly for these practical problems in engineering which involve complicated geometries and boundary conditions. Thereby, predictions of hydraulic fracturing processes rely heavily on the available numerical techniques.
In recent years, numerous computational methods have been developed to model hydraulic fracturing processes in rocks that can be roughly divided into discrete and continuous categories. The discrete approaches attempt to capture the exact topology of fracture either in an explicit way via remeshing techniques such that mesh edges are potential crack surfaces [3][4][5] or directly using Lagrangian particle methods or combined Eulerian and Lagrangian methods such as discrete element methods [6], smooth particle hydrodynamics [7] and material point method [8,9], or in an implicit manner such as the enrichment of shape functions that used in the extended finite element method [10][11][12]. The fluid flow process in the fracture is taken into account according to Reynolds' equation of lubrication theory [13,14]. It is noticeable that these discrete approaches either require special treatments for complex fracture topologies or sophisticated numerical schemes for discontinuities that considerably limits their applications in engineering practice. Alternatively, sharp fractures in smeared approaches are diffused via smooth transitions between intact and fully damaged material states such that discontinuities at fractures are removed. Typical continuous approaches include the peridynamic method [15][16][17], the gradient damage method [18], and the phase-field method [19][20][21], to name a few. In addition, another type of methods combined the continuous and discontinuous methods, such as the combined finitediscrete element method [22], and the numerical manifold method [23][24][25][26]. These discrete and continuous methods or the combined approaches all have been applied for modeling of hydraulic fracturing. Recently, the phase-field method for modeling of fracture is gaining an increasing popularity because it combines the advantages of the smeared and the discrete approaches and gives a precise meaning to the idea of using damage models to approximate the discontinuous fractures. The phase-field method introduces an auxiliary continuous scalar variable, referred to as phase field, to smear and regularize the sharp fracture discontinuities. Due to this feature, complex fracture phenomena for example crack branching and coalescence can be tackled without computational difficulties. Moreover, the phase-field approach was developed on the basis of variational principle for brittle fracture [27]. Regarded as a general Griffith's theory, this principle is capable of predicting both the crack initiation and the crack propagation path which further enhances the attraction of the phase-field approach. Owing to the advantages of the phasefield approach, numerous efforts have been devoted to exploring its possibility for modeling hydraulic fracturing propagation.
The first attempt to simulate hydraulic fracture using the phase-field approach was made by Bourdin et al. [20] that a term related to the fluid in cracks was included in the variational principle to consider fluid-pressure-induced fractures in impermeable materials. Phase-field models for hydraulic fracturing in saturated permeable porous media were then developed by Miehe et al. [19], Lee et al. [28], and Mikelic et al. [29,30] based on Biot's equation of poroelasticity. Xia et al. [31] and Shu et al. [32] further considered the material heterogeneity and the dynamic effects in the phase-field model for hydraulic fracture in porous media. Santillan et al. [33,34] proposed another framework for hydraulic fracture in both impermeable and permeable solids using phase-field approach with fluid flow in fracture being a lower-dimension manifold. In cases of permeable solids [34], the fluid flux between the fracture and the porous solid is employed as a coupling variable. These studies have shown that the phase-field method is a promising approach for modeling hydraulic fracturing propagations. Despite their contributions, these studies focused majorly on the theoretical and numerical development of the phase-field model for hydraulic fracturing. Many issues such as the effects of material properties on the hydraulic fracturing, interactions between the hydraulic fracture and a pre-existing fracture, and hydraulic fracturing propagation in natural rocks with complex fracture networks remain open questions and require further studies.
In this paper, a phase-field model is developed for hydraulic fracturing. Implemented using the finite element method, it is applied to study the hydraulic fracturing propagation in rocks with pre-existing natural fractures. Specifically, the hydraulic fracturing propagation in a poroelastic rock sample with two pre-existing fractures is investigated from the viewpoint of phase-field modeling. Emphases are placed on the influence of the initial distance between the fractures and the permeability of rocks. Additionally, the hydraulic fracturing propagation in a rock sample with more complex natural fracture networks is studied, in which interactions between multiple natural fractures and hydraulic fracture are concerned.
2 Phase-field model of hydraulic fracturing in permeable rocks The governing equations for the phase-field modeling of hydraulic fracturing in poroelastic media are introduced.

Governing equations for poroelastic media
Assuming that the solid skeleton is elastic and the fluid within the porous is incompressible, the deformation of a saturated porous rock under quasi-static loading is described by Biot's theory of poroelasticity [35,36].
The equilibrium equation for the saturated media is where b is the body force and σ is the total stress decomposed as with σ eff , p and α being the Biot's effective stress, the pore fluid pressure, and Biot's effective stress coefficient, respectively, and 1 = [1 1 0] T in two-dimensional cases. The effective stress and the linear strain ε is linked via the elastic constitutive equation where D is the elastic constitutive matrix which, in planestrain cases, is written as in which E is Young's modules and υ is Poisson's ratio. According to Biot's theory [35,36], the continuity equation for fluid flow is where q is the fluid flux, Q is the source/sink term, and θ is the fluid content defined as In the above equation, u is the displacement field of the solid skeleton and M is the Biot's modulus. The migration of fluids in porous media is described by Darcy's law that where k is the vector consisting of permeability coefficients of the material. Substitution of Eqs. (6) and (7) into the fluid continuity Eq. (5) results iṅ Complemented with the boundary conditions, the governing equations for poroelastic media are summarized as

Phase-field evolution equations for poroelsatic media
The theory of the phase-field model is constructed on the basis of the variational principle proposed by Francfort and Marigo [37] which interprets fracture as a competition between the surface energy for crack formation and the bulk energy stored in materials. In this section, the phase-field evolution equations for rocks are derived following [19] by first defining the total energy Ψ total as a sum of the bulk energy Ψ bulk stored in the rock and the energy dissipated in the fracturing process Ψ fracture that For a saturated poroelastic rock with a domain Ω bounded by its surface ∂Ω, the bulk energy is in the form of where the density of bulk energy is In this work, we only consider the quasi-static loading condition and ignore the dynamic responses in the hydraulic fracturing process. Acording to Miehe et al. [19], the elastic energy density and the contribution from the fluid are calculated as The total stress and the fluid pressure are accordingly the derivatives of φ bulk with respect to the strain and the fluid content that are and Substituting Eqs. (12) and (14) into Eqs. (15) and (16) leads to and which are apparently in consistence with the corresponding definitions in the classical Biot's theory, for instance Eqs. (2) and (6). According to Francfort and Marigo [37], the fracture energy Ψ fracture is defined as where G c is the Griffith critical energy release rate and Γ is the crack set. To overcome the numerical difficulties induced by the discontinuity, the crack is regularized by introducing an auxiliary scalar variable, referred to as the phase field s(x, t) ∈ [0, 1]. The scalar variable, s, indicates the damage degree of the material with s = 1 and 0 indicating the fully damaged and the intact states, respectively. The sharp crack thereby is smeared out as a diffusive fracture. Specifically, the discontinuous fracture surface Γ(s) [19,21] is replaced by with l 0 being the length scale, and consequently, the energy required to generate a new crack is rewritten as In order to describe the stiffness degradation due to the damage process, a degradation function g(s) = (1 − s) 2 was introduced by Miehe et al. [38], which gives g(0) = 1, g(1) = 0 for representing the constraints of intact when s = 0 and fully broken when s = 1, respectively. In addition, the derivation of the degradation function equals to 0, i.e., g ′ (s) = 0, which ensures the derivation of the energy density with respect to the phase field, ∂ s ψ, converges to a final value if the damage converges to the fully broken state s = 1. Here The stress-strain relationship considering damages is expressed as which recovers the classic linear elastic constitutive relationship when material are intact (e.g., s = 0). Using Eqs. (14), (21), and (22), the total energy functional in Eq. (10) is written as Sequentially, the total energy density function is in the form The phase-field evolution equations are thus derived via the derivative of the total energy density function with respect to the phase field, for example ∂φ total ∂s ¼ 0, leading to As shown, the phase-field evolution is driven by the elastic strain energy function of the material, φ e , which however results in three unrealistic fracture phenomena: (1) the energy density φ e is not a monotonously increasing function meaning that the fracture is recovered if φ e decreases; (2) the energy density induced by pure compression may also lead to fracture; and (3) a material may be more-or-less damaged as long as φ e is greater than zero.
To overcome the first two issues, Miehe et al. [19] decomposed the elastic strain energy function into tension and compression parts as by compression, based on the two ramp functions 〈x〉 + = (x + |x|)/2 and 〈x〉 − = (x − |x|)/2 and the decomposition of the strain tensor with λ and μ are lame parameters of the solid skeleton. A history was then suggested [19] to drive the evolution of the phase field that Since H + is monotonously increasing and contributed by only the tension part, it prevents from the fracture reversibility and the compression-induced fracture.
In the damage evaluation Eq. (29), the maximum energy density H + is a monotonously increasing function of the strain ε. In order to further prevent from the unnecessary material stiffness reduction at very low stress states, Miehe et al. [19] introduced a damage evolution criterion with a threshold by using the fracture energy density with σ c being the critical effective stress of matrix. The energy criterion with threshold can be rewritten as The total energy density decribed by Eq (25) can be rewritten as The evolution equation is then expressed as [19,39].
where the driving energy is defined as H max = max(φ + − ψ c ) + which poses a threshold and thus resolves the third issue. The evolution equations for the phase field with the corresponding boundary conditions are summarized as follows 2.3 Deformation-dependent permeability for the fracture The migration of fluids in the porous media is assumed to be a type of Poiseuille fluid flow and the phase field is used to determine the deformation-dependent permeability at the fracture. The permeability tensor k in Eq. (8) is decomposed into two parts [19], expressed as where ϵ is an parameter used to localized the permeability in the fracture, k m = k/η1 is the permeability tensor of matrix (k is the intrinsic permeability, η is the fluid dynamic viscosity), and k f is the anisotropic permeability depending on the fracture aperture defined as where 1 is the identity 2-order matrix, n s ¼ ∇ s ∇ s k k is the normal vector defined by the phase field, and w = h ∇ u • n s is the fracture width representing the displacement jump with h being a length parameter (approximated by the mesh size in practice [19,40]). The normal displacement jump along the fracture interfaces can be expressed as w n = w • n s according to [31].

Finite element discretization and staggered solution algorithm
In this section, the governing equations of the coupled problem are discretized using the standard Galerkin finite element method. The quadrilateral bilinear elements are applied to approximate the displacement field u, the fluid pressure field p and the phase field s, namely where b u, b p and b s are the vectors for the displacement, pressure, and phase field at element nodes and N u , N p and N s are the shape functions.
Applying the above approximations to Eq. (9), we have where the coefficients in the matrix are In Eq. (39), D ′ is the modified elastic constitutive matrix, written as The matrices B u and B p are the derivatives of shape functions expressed as and The backward Euler finite difference scheme is adopted for where Δt is the time increment, u n + 1 and p n + 1 denote the displacement and pressure at time t n + 1 , and u n and p n denote the displacement and pressure at time t n . Substituting Eqs. (48) and (49) into Eq.(38) arrives in Similarly, the phase field equation can be discretized as: where and The matrix B s is the derivative of shape functions for the phase field that is In summary, the discretized governing equations include Eqs. (38) and (51) that can be solved using a staggered scheme. In detail, in a typical time increment the displacement and pressure of poroelasticity problem is solved with a fixed phase field based on Eq. (50). Then, the phase field is resolved using Eq. (51) with the updated displacement field and pressure field. A detailed solution algorithm is summarized in Fig. 1. The staggered solving scheme is a kind of explicit approach; the smaller time steps imply the more stable the solution process; however,

Experimental tests and discussions
The applications of the developed phase-field model for modeling hydraulic fracturing in saturated poroelastic rock will be carried out in this section. Three examples with increasing complexity are concerned which are the propagation of a single hydraulic fracture, the hydraulic fracturing propagation interacted with a natural fracture, and the hydraulic fracturing in natural fracture networks.

Validation
Firstly, we validate the phase-field model and verify the developed computational codes by comparing numerical simulation result with the classical analytical solution by Sneddon [41]. In this validation example, we consider a single fracture with length 10 m in a square domain of dimensions 100 m × 100 m as shown in Fig. 2.
The domain is assumed to be linear elastic, homogeneous, and isotropic. The analytical solution exists for the fracture aperture under constant internal pressure, given by [34].
In this validation example, the displacement is fixed in all directions on the external boundaries, i.e., u = 0, v = 0. T h e m e c h a n i c a l p r o pe r t i e s a r e se t a s Yo u n g 's modulus E = 16.06GPa and Poisson's ratio ν = 0.1812. The pressure applied in the fracture is p =1 × 10 5 Pa. To keep the same condition as the assumption of the analytical solution, we set the intrinsic permeability of the matrix as k = 0. Biot's coefficient is assumed as α = 0 in the matrix and α = 1 in the fracture, which assumes that the fracture propagation is only contributed by the pressure applied in the fracture. The domain is discretized using four-node quadrilateral finite elements with a mesh size being 0.25 m. Following [19], the length scale of the phase-field model is set to be l 0 = 3h, where h is the mesh size.
The simulated vertical displacement field for the case of mesh size 0.25 m and comparison of fracture aperture between the analytical solution and numerical result are presented in Fig. 3, showing that the vertical displacement field is symmetric along the fracture. The numerical simulated fracture aperture for the case of mesh size 0.25 m matches very well with the analytical result (the red curve). The convergence of the simulation results with respect to mesh sizes is also evaluated by analyzing the problem using meshes of different sizes (e.g., h = 0. 25 m, 0.4 m, and 0.5 m). With the decrease of mesh sizes from 0.5 to 0.25 m, the simulated fracture aperture gradually converges to the analytical result. This validation result shows that the phase-field model developed in this study is accurate and the numerical results simulated by using the relatively small mesh size, i.e., h = 0.25 m presented in this work are reliable.

Propagation of a single fracture
After validation, we further test the feasibility of phase-field model for hydraulic fracturing of a single fracture in porous media. The domain geometry and the fracture location is the Fig. 4 Phase-field evolution with fluid injection at 0.5 s, 1.5 s, 3 s, and 8 s  Fig. 2. To drive the fracture, fluids are injected into the pre-existing fracture (10 m in length) with a specific constant flowrate 0.005 m 2 /s. The properties of the rock adopted in the numerical simulation are from [42] and listed in Table 1. The length scale of the phase-field model is set to be 0.75 m according to [19]. The domain is discretized using four-node quadrilateral finite elements with a mesh size being 0.25 m. The time step in the simulation is Δt = 0.1 s. The displacements of all external boundaries are fully fixed (e.g., u = 0, v = 0), and the initial porous fluid pressure is set to be zero representing drainage boundaries. Figure 4 illustrates the snapshots of the phase field at four different times, i.e., t = 0.5 s, 1.5 s, 3 s, and 8 s and the distributions of the corresponding porous fluid pressure are shown in Fig. 4. As indicated, the hydraulic fracture propagates symmetrically in both directions due to the continuous injection of fluids (Fig. 4). The fluid pressure in the fracture increases at the beginning (Fig. 5a, b) because of the supply of fluids which then drops considerably when the crack starts to propagates (Fig. 5b-d. It is notable that the entire domain is saturated in the simulation which means the fluid lag is not considered in this model. However, the fluid pressure at tips of fracture is very low and even negative as shown in Fig. 5 indicating that a mathematical fluid lag exists. This phenomenon was also observed in the modeling of hydraulic fracture in saturated porous media using the phase-field model [29], the cohesive fracture model [43], and extended finite element method [44]. This variation can be observed precisely by showing the curve of the fluid pressure at the center of the sample against time (Fig. 6). A linear increase of the fluid pressure is observed at the initial stage, for example when t < 0.9 s. After reaching the maximum value (around 6 × 10 5 Pa) at 1.2 s, the fluid pressure decreases gradually to the asymptotic value at around 4.5 × 10 5 Pa. The simulation results presented are in line with these in [31] which verifies the developed phase-field approach for modeling hydraulic fracturing.
The evolutions of the fluid pressure at the center of the fracture from the simulation result are shown in Fig. 5. The pressure increases dramatically at the initial stage when t < 1 s and then gradually decreases until it approaches to 4 × 10 5 Pa. The decrease of pressure is caused by fracture propagation.

Interaction of hydraulic fracture with a pre-existing natural fracture
In the second example, we focus on the hydraulic fracture propagation in a poroelastic rock interacted with a pre-existing natural fracture. To this end, the same sample considered in section 4.1 is re-analyzed except that a pre-existing natural fracture is further considered in the sample (Fig. 7). Material parameters are the same as these listed in Table 1 if not otherwise specified. The influences of the initial distance between the hydraulic fracture and the natural fracture (e.g., the distance L in Fig. 7), the permeability of the rock material, and the injection rate on the fracture propagation process are investigated.   Figure 8 shows the evolution process of the phase field at different time points (i.e., t = 0.5 s, 1.5 s, 3 s, and 8 s from the left to the right columns) for different cases (i.e., L = 0 m and 10 m from the top to the bottom rows). In all cases, the injection of the fluid first leads to the coalescence of the two fractures followed by the further hydraulic fracturing propagation from the two tips of the pre-existing natural fracture. For L = 0 m, since the natural fracture is already connected with the fracture at the center, the fracture propagates from the tips of natural fracture when fluids are injected. It is noticeable that, for all cases, the crack propagation from the tips of the natural fracture is not directed horizontally but with an inclined angle of θ with respect to the horizontal direction (see also Fig. 7b).
Figures 9 and 10 present distribution of phase field and fluid pressure after injection of total specific volume for different distances between the hydraulic fracture and the natural fracture, different values of rock permeability, and different cases of injection rates, respectively. Comparing the phasefield distributions for different distances between the hydraulic fracture and the natural fracture, the fracture propagation lengths are generally shorter than the cases when the distance is larger. The inclined angle θ slightly decreases with the increase of distance, which is caused by the different pressure values when the natural fracture start to propagate. For the smaller rock permeability value, i.e., k = 2 × 10 −14 m 2 , the fracture propagation lengths are generally larger than that of Fig. 8 From left to right column is phase-field distribution at t = 0.5 s, 1.5 s, 3 s, and 8 s, and from top to bottom row is phase-field evolution with different distances between hydraulic fracture and natural fracture Fig. 7 a Geometry and boundary conditions of pre-existing hydraulic fracture and natural fracture. b Interaction between hydraulic fracture and natural fracture by fluid injection the higher rock permeability, i.e., k = 2 × 10 −13 m 2 . It indicates that the rock permeability does influence the speed of hydraulic fracturing propagation-the less the permeability is the faster the crack propagates. This can be explained by showing the distribution of the pore fluid pressure. A higher permeability leads to a more diffused zone with fluid migrations whereas the zone of high fluid pressure is concentrated at the fracture when the permeability is low. The maximum pressure in the fractures for cases of high permeability is thus smaller than that for cases of low permeability, given that the same injection rate is applied. The fluid pressure in the fractures is generally much higher for the cases with higher injection rates (Fig. 10). When L = 0 m, the injection rate has limited impact on the inclined angle θ for the natural fracture propagation. However, when L = 10 m, the cases with higher injection rate, i.e., q = 0.002 m 2 /s, lead to fracture branches in tips of the natural fracture (Fig. 9), due to relatively high elastic energy during the injection.
The variations of the fluid pressure at the center of the sample for all cases are shown in Fig. 11 which echoes the explanation for the fracture propagation shown in Fig. 9. Noteworthy, in case of L = 0, the fluid pressure in both cases of injection rate decreases gradually after the peak value in Fig. 11 whereas an apparent increase of the fluid pressure is observed after it drops from the peak for the case of L = 10 m. Such a fluctuation for L = 10 m stems from the crack k=2e-13 Distance L=10 Fig. 9 The phase field at the same specific total volume, with constant flowrates at 0.0005 m 2 /s at 12 s and 0.002 m 2 /s at 3 s in different cases of rock intrinsic permeability, i.e., k = 2 × 10 −14 m 2 (the top two rows) and k = 2 × 10 −13 m 2 (the bottom two rows) for different distances, i.e., L = 0 m and 10 m coalescence that results in a sudden change of fluid pressure distribution in the hydraulic fracture and the pre-existing fracture. For L = 0, the two cracks are already connected at the beginning, and thus, no fluctuation emerges. Additionally, the initial space to store the fluid for L = 0 is larger; thereby, the increase rate of the fluid pressure in the fracture is much smaller at the initial stage (e.g., t < 1.8 s when q = 0.0005 m 2 /s and t < 0.6 s when q = 0.002 m 2 /s).

Interaction of hydraulic fractures with natural fracture networks
In natural rocks, fractures are commonly presented in the form of networks of more complex topology. To be in line with actual situations, a rock sample with 13 randomly distributed natural fractures shown in Fig. 12 is analyzed as the third example. Material parameters are the same as those in Table 1 if not otherwise specified. Fluids are injected into a borehole (e.g., F0 in Fig. 12 located at the center of the sample with a constant flowrate at 0.001 m 2 /s. Similar to previous numerical examples, the displacements are fixed for all external boundaries, and the hydraulic pressure there is set to be zero. Figure 13 presents the distribution of the phase field at different times, i.e., t = 0.9 s, 2.3 s, 4.9 s, 7.1 s, 9.1 s, and 14.1 s, respectively. The result shows that the perforation fractures F0 continuously propagate, when fluids are injected, and gradually connect to the nearest natural fractures F3, F9, and F1 (Fig. 13a, b). Afterward, the connected natural fractures  Fig. 10 The phase field at the same specific total volume , with constant flowrates at 0.0005 m 2 /s at 12 s and 0.002 m 2 /s at 3 s in different cases of rock intrinsic permeability, i.e., k = 2 × 10 −14 m 2 (the top two rows) and k = 2 × 10 −13 m 2 (the bottom two rows) for different distances, i.e., L = 0 m and 10 m F3, F9, and F1 start to propagate to further link to other nearby fractures. For instance, fracture F3 propagates and connects to fracture F2 (Fig. 13c, d), and fracture F9 firstly connect to fracture F3 with its other side being connected to fracture F8 later (Fig. 13d, f). At the time of 14.1 s, all fractures except fracture F6 are connected. It demonstrates that hydraulic fracturing by injecting fluid through the borehole fracture F0 plays an important role in the formation of the fracture networks. The effective connectivity of fracture networks is critically enhanced due to the continuous hydraulic fracturing propagation which consequently enhance the effective permeability of the rock sample. Figure 14 illustrates the fluid pressure distributions at different injection times. The distribution of fluid pressure changes dramatically throughout the process which is owing to the pressure redistributions along the propagating fractures. In particular, the fluid pressure at fractures increases much faster than the surrounding undamaged rock matrix. The relatively higher pressure zones are concentrated in the fractures that are connected to the injection borehole, i.e., the perforation fractures F0. The evolution of pressure distribution implies that, in the hydraulic fracturing process in fractured rocks, the pressure distribution may contain significant uncertainties, which primarily depends on the damage degree of the rock and the connectivity of fracture networks.
The evolution of fluid pressure at the center of fractures near the injection borehole, i.e., the fracture F0, F1, F2, F3, F6, and F9, is presented in Fig. 15, showing complex evolution patterns. The pressure in the perforation fractures F0 (at the injection borehole) increases at the initial stage of injection (e.g., when t < 1 s), and then drops down considerably when the perforation fractures connect to the adjacent natural fractures. For the pre-existing fractures, generally the pressure values increase dramatically once they are connected to the propagating hydraulic fractures, e.g., the fractures F1, F3, and F9. Due to the deformation, the zone at front of the propagating crack tips (see also Fig. 14) may pose negative fluid pressure which has also been observed in former contribution [29,31]. The negative fluid pressure is also encountered at some natural fractures, for example, the fluid pressures at F1 and F9 drop down to negative values at the initial stage as shown in Fig. 15, which however increases gradually due to the migration of flow in the fractures and rock matrix. The simulation results indicate that the evolution of pressure is significantly affected by the propagation of hydraulic fractures. The Fig. 11 The evolution of fluid pressure at the center of the hydraulic fracture for different cases of the rock permeability, i.e., k = 2 × 10 −14 m 2 and k = 2 × 10 −13 m 2 , and different distances near the natural fracture, i.e., L = 0 m and 10 m  To investigate the influence of rock permeability on the hydraulic fracturing propagation in fracture networks, the problem is re-analyzed with permeability k = 2 × 10 −14 m 2 , 1 × 10 −13 m 2 , and 5 × 10 −13 m 2 , respectively. The distributions of the corresponding phase field and fluid pressure at t = 14.1 s are illustrated in Fig. 16 and Fig. 17. As seen in Fig. 16, minor discrepancies are observed for the phase-field distribution; however, the distribution of fluid pressure differs a lot from each other as shown in Fig. 17. A higher permeability leads to a more diffuse fluid pressure distribution, which is similar to the observations presented in section 4.2. Figure 18 demonstrates the variation of fluid pressure at the center of fracture network sample with different intrinsic permeability. The fluid pressure has similar trends for all cases; however, a smaller permeability of rocks results in a relevantly larger fluid pressure at the borehole.

Conclusions
In this paper, a phase-field model is developed and implemented for investigating the hydraulic fracturing process in saturated poroelastic rocks. Particularly, the hydraulic fracture interacted with a pre-existing natural fracture in a rock is studied with emphasis on the influence of the fracture distance and the permeability of rocks on the crack propagation process. Moreover, interactions between fractures in fracture networks, which are of more complex topology and highly in line with the actual situations, are concerned. The main conclusions drawn from this study are as follows: & The phase-field approach alleviates the difficulties in modeling hydraulic fracturing caused by discontinuities across fractures and thus is capable of tackling problems of complex fracture networks. & The distance between the hydraulic fracture and the near natural fracture has a significant impact on the propagation direction of the natural fracture once it is connected to the fracture where fluids are injected. & The rock permeability has a limited impact on the propagation directions. However, it affects the propagation speed considerably. The smaller the rock permeability, the faster the hydraulic fracture propagates. & The injection rate has significant impact on the fracture propagation process, which may cause fracture branches in the propagation of natural fracture once it is connected to the hydraulic fracture. & The hydraulic fracturing process will significantly increase the connectivity of the natural fracture networks due to the propagation of hydraulic fractures, which consequently increases the effective permeability of the rocks. & The evolution of pressure in the fracture networks is affected by the propagation of hydraulic fractures. The evolution of the pressure in the injection fracture could be  Although two-dimensional cases are concerned in this study, the extension of the proposed phase-field model to three-dimensional cases is straightforward. Additionally, the hydraulic fracture propagation, fluid pressure evolution, and interaction with natural fracture networks are studied in the presented numerical examples. Many factors, such as the in situ stress boundary conditions, the more realistic rock and fluid properties, the operational parameters, the dynamic responses and mixed fracture mode under the combined action of shear, and compression stresses [45] in the hydraulic fracturing process, have not been considered in this study yet and require efforts for detailed investigation in the future.
Funding information This study received support by the Australian Research Council Discovery Project funding scheme (Project Number DP150104257).
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://creativecommons.org/licenses/by/4.0/. Fig. 17 The pressure field at the injection time t = 14.1 s in different cases of rock permeability, i.e., k = 2 × 10 −14 m 2 (a), k = 1 × 10 −13 m 2 (b), and k = 5 × 10 −13 m 2 (c) Fig. 18 The evolution of fluid pressure at the center of the fracture network sample for different cases of the rock permeability, i.e. k = 2 × 10 −14 m 2 , k = 1 × 10 −13 m 2 , and k = 5 × 10 −13 m 2