On the numerical simulation of crack interaction in hydraulic fracturing

In this paper, we apply the enhanced local pressure (ELP) model to study crack interaction in hydraulic fracturing. The method is based on the extended finite element method (X-FEM) where the pressure and the displacement fields are assumed to be discontinuous over the fracture exploiting the partition of unity property of finite element shape functions. The material is fully saturated and Darcy’s law describes the fluid flow in the material. The fracture process is described by a cohesive traction-separation law, whereas the pressure in the fracture is included by an additional degree of freedom. Interaction of a hydraulic fracture with a natural fracture is considered by assuming multiple discontinuities in the domain. The model is able to capture several processes, such as fracture arrest on the natural fracture, or hydraulic fractures that cross the natural fracture. Fluid is able to flow from the hydraulic fracture into the natural fracture. Two numerical criteria are implemented to determine whether or not the fracture is crossing or if fluid diversion occurs. Computational results showing the performance of the model and the effectiveness of the two criteria are presented. The influence of the angle between a hydraulic fracture and a natural fracture on the interaction behaviour is compared with experimental results and theoretical data.


Introduction
The technique of hydraulic fracturing can be used to stimulate low permeable gas and oil reservoirs. Generating a large fracture network by injecting fluid under high pressure into the hydraulic fractures enhances the permeability greatly and thus increases production rates [1]. Additionally, this technique can also be used for enhanced heat recovery purposes [2]. In order to achieve an optimal fracture network in realistic situations, it is necessary to correctly describe the interaction between hydraulic fractures and the pre-existing natural fracture network. Tectonic stress rotations tilt the natural fracture network at the time of the fracture formation which may result in tilt fractures that are not aligned with the maximum horizontal stress. Natural fractures may therefore intersect with hydraulic fractures, altering the propagation path leading to complex fracture geometries [3].
Lee et al. [4] demonstrated experimentally the influence of calcite filled natural fractures on the propagation path of a mechanically induced fracture in shale rocks. Blanton et al. [5] showed that an induced hydraulic fracture may either cross a natural fracture, arrest onto the fracture, or fluid diverts into the natural fracture. These different hydraulic pathways interacting with a natural fracture are shown Fig. 1. The experiments show that the inclination angle between the natural fracture and the hydraulic fracture in relation to the in situ stress differences is important in determining which of the phenomena occurs. Zhou et al. [6] performed similar experiments and demonstrated the influence of the friction in the natural fracture.
In a theoretical study, Renshaw and Pollard [7] studied a propagating fracture across an orthogonal interface and derived an analytical criterion, which was validated by experiments. The criterion is based on linear elastic fracture mechanics (LEFM) under the assumption that no slip occurs Fig. 1 Possible pathways due to a hydraulic fracture (HF) interacting with a natural fracture (NF) along the interface before the two fractures merge. The possibility of the fracture to divert into the natural fracture is not considered in this study. The criterion has been validated and extended by means of experiments for non-orthogonal angles by Gu et al. [8] by solving the criterion numerically.
Another mechanism that causes crack tip branching occurs when the fracture propagation speed exceeds the Rayleigh wave speed of the material [9,10]. Valko and Economides [11] showed that hydraulic fractures propagate more slowly than the Rayleigh wave speed criterion. Therefore, branching of a single hydraulic fracture does not appear and does not have to be considered.
The extended finite element method (X-FEM) is a proven concept in numerical modelling of crack propagation and has the advantage that there is no need to remesh the domain [12], i.e. the topology of the original mesh does not need to be modified as the crack evolves. By exploiting the partition of unity property of finite element shape functions, a crack is represented as a discontinuity in the displacement field [13]. The discontinuity can simply be placed in arbitrary locations in the finite element mesh by adding additional degrees of freedom to existing nodes [14]. Daux et al. [15] showed that multiple discontinuities can be included in a similar manner by stacking up one additional set of degrees of freedom per discontinuity.
Dahi-Taleghani and Olsen [16] developed an X-FEM based model for the interaction between hydraulic and natural fractures. The influence of a diverted hydraulic fracture was compared to the fracture wing that did not interact. Khoei et al. [17] used a similar approach but included frictional contact based on plasticity theory of friction using a penalty method. Similar to [16], fracture crossing was not considered yet in their work. In addition, these X-FEM models did not consider the porosity of the bulk material. Nevertheless, the method was applied successfully for fracturing in porous materials, see e.g. [18][19][20], also in combination with a single hydraulic fracture [21]. The X-FEM has also been successfully applied to model fluid flow through existing fracture networks in a porous media [22,23]. Other methods that considered hydraulic fracturing in porous media are, e.g. based on interface elements with a cohesive zone [24], remeshing techniques [25] and phase field approaches [26]. The phase field method incorporates complex fracture patterns and 3D fractures directly into a model based on the variational approach [27]. X-FEM requires additional implementation aspects to include these into a model but it does allow the use of more coarse meshes.
The enhanced local pressure (ELP) model, based on the X-FEM approach, has been developed specifically for hydraulic fracturing in very low permeable rocks [28]. The pressure is assumed to be discontinuous over the fracture to prevent the necessity to resolve the very steep pressure gradient near the fracture surface. By including an additional degree of freedom to account for the pressure in the fracture, it is ensured that all the injected fluid goes into the fracture. Fluid leakage is included by an analytical solution based on Terzaghi's one-dimensional consolidation equation.
In this study, we extend the ELP model to account for multiple, interacting fractures. These fractures are included by adding a set of additional degrees of freedom for each fracture. Upon interaction of a hydraulic fracture with a natural fracture, the hydraulic fracture can either cross the natural fracture or fluid can divert into the natural fracture. The former is based on an average stress criterion and the latter on the opening displacement of the natural fracture. Both can happen simultaneously and are based on the numerical results so no additional theoretical criterion is needed. Compared to previously mentioned methods, we include fracture crossing and fluid diversion simultaneously. The advantage of X-FEM, i.e. fracture growth irrespective of the underlying finite element mesh, is exploited with the ELP model. With the proposed model, it is possible to predict the interaction mechanics between a hydraulic fracture and a natural fracture network in a poro-elastic material. We demonstrate the performance of the model by investigating four numerical examples.

Model background
In the following section, we describe the background of the numerical model. The section is subdivided into two parts, the numerical formulation (paragraph 2.1-2.5) and the numerical implementation (paragraph 2.6). The kinematic relations for the displacement and the pressure are described in a mathematical framework in paragraph 2.1. The balance equations in the porous material, at the discontinuity, and the weak form of the problem are given in the paragraphs 2.2-2.4, respectively. In paragraph 2.5 the weak form is discretized with standard finite element shape functions based on the partition of unity property of the shape functions [13]. Numerical implementation aspects of the X-FEM hydraulic fracturing model are discussed in Section 2.6.

Kinematic relations
We follow the kinematic relations as described in [28]. Consider the body , which contains m discontinuities, see Fig. 2a. Each discontinuity i separates the domain in a two parts, + i and − i . We can write the total displacement field at any time t, following a discrete modelling approach [12,14,29], as a continuous displacement field u(x, t) and m additional displacement fieldsũ i (x, t). The total displacement field can be written as [15] u( where x is the position of a material point and H d i is the Heaviside step function defined across discontinuity i as The pressure field is decomposed in a similar fashion in a continuous fieldp(x, t) and m discontinuous pressure fields At the discontinuity, inside the crack, the pressure is equal to variable p d i (Fig. 2b).
The displacement jump across discontinuity i is written as where the notations + and − are used for the same location but located compared to the positive and the negative side of discontinuity, respectively, see Fig. 3. This can be rewritten as The first term in this equation is the jump caused by discontinuity i. The second term is the jump caused by the other discontinuities (j ) in the domain.

Balance equations
The porous material is considered to be fully saturated with fluid. Assuming that the contribution of gravity, inertia and convection can be neglected, the momentum balance is written as  Here σ is the total stress, which is decomposed in Terzaghi's effective stress σ e and the hydrostatic pressure p according to [30]: In this equation, α is the Biot coefficient and I is the unit matrix. The Biot coefficient is defined as with K and K s being the bulk moduli of the porous material and the solid constituent, respectively. Neglecting mass transfer between the two constituents results in the following equation for the mass balance [18] α∇ · v s + ∇ · q + 1 where v s is the deformation velocity of the solid skeleton, q is the seepage flux, and M is the compressibility modulus defined as 1 The porosity of the porous material is defined by φ and K f is the bulk modulus of the fluid. The continuity equation for fluid flow within discontinuity i is given by [31] Here v d i is defined as a natural extension of the v s field in the porous material and the seepage flux q d i . We assume that the fluid flow can be described as a flow between two parallel plates. The hydraulic fracturing fluid is considered to be a Newtonian fluid. Under these assumptions the fluid flow in the discontinuity is defined as [32] where k d i is the permeability in discontinuity i, defined by the normal opening of the discontinuity and u n i and by the dynamic viscosity μ [32]. We refer to [28] for more details about the balance equation and the corresponding boundary conditions.

Constitutive law at the fracture and the interface
The constitutive mechanical behaviour at a fracture is described by a relationship between the traction at the interface and the displacement jump u d across the fracture [31]: where κ is a history parameter that is equal to the largest displacement jump reached. It is necessary to perform a linearisation on Eq. 14 in order to use the tangential stiffness matrix in an incremental iterative solution: The relation between the traction t d and the displacement jump u d can be any traction-separation relation and is referred to as the cohesive law. We assume that hydraulic fractures open in mode-I due to the internal pressurization. The shear tractions are therefore neglected and an exponential cohesive law is used that is only a function of the normal opening u n [33] Here is τ ult the ultimate strength of the material and G c the fracture thoughness. In contact, we assume a penalty constraint in both normal and shear direction. The linear relation between the traction and the opening displacement is defined by a stiffness parameter D n and D s for the normal and shear direction, respectively. In contact this gives the following penalty relations In the remainder of this study, natural fractures are described by these contact relations. We do not consider additional cohesion in the natural fractures due to filling materials. However, the framework allows this to be added later.

Weak form
We derive the weak form for multiple discontinuities by multiplying the balance equations with admissible test functions in the same form as the displacement field and the pressure field as where η and ζ are the admissible displacement and pressure fields, respectively. Multiplying the momentum balance in Eq. 7 with test function η, including boundary conditions and using Gauss's theorem, gives the weak form of the momentum balance: Here, we use the notations for the Heaviside function of variational field k integrated over discontinuity j for the positive and negative side, respectively. The traction at the interface is equal to The external applied traction on the body is given by t p (see Fig. 2a). Separation of the tractions into one part where the variational displacementη k is acting on discontinuity d d k and into one part for the remainder of discontinuities, i.e. j = k gives Multiplying the mass balance in Eq. 10 with test function ζ gives where f f is the external applied fluid flux (see Fig. 2). These equation must be satisfied for all the variations of η and ζ . Separating the balance equations in a continuous The same can be done for the m discontinuous equations, with (η = 0 andζ = 0), giving which must hold for k = 1...m. The last balance equation is conservation of mass for the fluid flow in the discontinuity given in Eq. 12. In a weak form, multiplied by a test function ψ, and integrated over the fracture this is written as The first two terms represent the analytical calculated fluid leakage. The third term is the tangential fluid flow based on lubrication theory. Term four and five are representing the opening rate terms in normal and shear direction, respectively, the sixth term is the compressibility of the fluid within the fracture, and the last term is the fluid injection within the fracture. The vector s d j represents the tangential vector at discontinuity surface j . In the remainder of this study, we neglect the fifth term representing the volume change due to elongation of the discontinuity surface in tangential direction. We assume that this contribution is small compared to the volume of the opening of the fracture in mode I.

Discretization
The spatial discretization of the balance equations is based on the partition of unity property of finite element shape functions as described in the work of Babuška and Melenk [13]. The displacement field, the pressure field, the pressure in the fracture and the variational forms are discretized similarly following the Bubnov-Galerkin approach for a single element by: where N, H , and V are matrices containing the standard shape functions for respectively, the nodal displacement, the pressure, and the pressure in the fracture for all nodes that support the element. Note that the shape functions for the nodal displacement and the pressure are two-dimensional functions while the pressure in the fracture is described in a one-dimensional domain (Fig. 4). The termsû andp contain the degrees of freedom of the continuous displacement and pressure fields, respectively. Whileũ andp contain the values of the enhanced nodes. The term p d contains the nodal values of the pressure in the discontinuity. We use the underline to distinguish between the field variable and the discrete values. The discretized strain in the bulk can be derived by differentiation where B contains the spatial derivative of the standard shape functions. Filling in the discretized form in the balance Eqs. 23-26 give the continuous equations as and k = 1...m discontinuous equations with the vector m in the two-dimensional situations being defined as m = ( 1, 1, 0 ) T . Finally the discretized form of the mass balance in the discontinuity is given as These final equations can be linearized in a standard way and solved using a Newton-Raphson iterative method. More details about the constitutive equations, linearization and time integration can be found in [28].

Numerical implementation
Each discontinuity is able to propagate. The position and the direction of propagation are governed by two unique level sets [34]. The direction of propagation is based on the Camacho-Ortiz equivalent traction [35] defined as Here, θ is the propagation direction and t n and t s are the normal and shear traction in the direction of θ . The parameter β is a scaling factor which is typically set at 2.3 [35]. The equivalent traction t eq can be calculated for every angle θ for a given stress state σ av . The angle for which the equivalent traction exceeds the maximum allowable traction τ ult of the material is the direction of propagation. Additional details about the equivalent traction can be found in [36]. We use a non-local approach to calculate an average stress defined as [33] Here, n int is the number of integration points in the domain, σ e,i is the current effective stress state in integration point i which has a Gaussian weighting function defined as with r i being the distance between the crack tip and the integration point n i , and l a being a length scale parameter defining how fast the weight factor decays as a function of the distance between the integration points and the crack tip.
A discontinuity is assumed to propagate in a straight line through an element and always ends on the element boundary or at another discontinuity. The discontinuity can propagate through multiple elements within one single timestep. Upon the interaction of two discontinuities, there are two requirements on the numerical implementation, i.e. (i) two discontinuities must be connected once a tip is in the vicinity of another discontinuity and (ii) the connecting tip must be stopped from further propagating. We determine the event of connecting the two discontinuities by counting the number of elements between a tip and the nearest discontinuity. The number of elements is an input parameter of the simulation in which zero means that two discontinuities connect when the tip propagates into an element that already contains a discontinuity. This is often an undesirable situation since a small distance between a tip and a discontinuity may lead to an ill-conditioned system. In our simulations we therefore chose to connect two discontinuities if there is one element in between (Fig. 5).
An additional crossing checkpoint is added at the opposite side of the interacting discontinuity. An average stress (Eq. 36) is calculated in the checkpoint and is used to determine whether the discontinuity crosses the existing crack (Fig. 5). Note that the stress is only averaged at the side of the checkpoint. A new discontinuity nucleates, and thus crosses, when the average stress violates the same criterion as was used for the determination of crack propagation (Eq. 35). The new discontinuity is given an initial length of three elements to prevent instantaneous interaction with the neighbouring discontinuity.
In the case of interacting fractures, the Heaviside enrichment given in Eq. 2 is no longer valid when there are two or more enrichments present in one element [15]. The common way to solve this problem is by introducing a special junction enrichment function present in the multiple  A finite element mesh containing two discontinuities. Integration points lying in the cut off region have a Heaviside value of zero for the green discontinuity. Nodes with a double enrichment, due to the two discontinuities, are shown in green enriched elements. This would lead to additional terms in the kinematic relations in Eqs. 1, 3 and 6. To avoid this, we implement a modified Heaviside enrichment by evaluating whether an integration point belonging to multiple discontinuities is cut off by one of the discontinuities (Fig. 6). Integration points that lie in the shaded purple region do not have a contribution to the displacement field of the green discontinuity. In these integration points, the values for the Heaviside of this discontinuity are therefore set to zero. Hence, the kinematic relations do not have to be modified.
In the model, we distinguish between two types of discontinuities. Discontinuities that possess the ELP degree of freedom are hydraulic fractures. The second type are discontinuities that do not possess this degree of freedom and can therefore be considered as a closed natural fracture. When two hydraulic fractures interact, the mass balance in Eq. 34 is combined. The possibility of fluid diversion occurs when a hydraulic fracture interacts with a natural fracture. The ELP degree of freedom p d is added to the interacting element but not in the remainder of the natural fracture (Fig. 7). Fluid can divert into the natural fractures if a diversion criterion is violated. We propose that the criterion is violated if the opening displacement in all integration points in an element is positive, as shown in Fig. 7.
The nucleation criterion for an element with an added checkpoint and the diversion criterion for an element in a natural fracture are formulated as:

Results
The performance of the model is demonstrated in four examples. Pore pressure is not considered in the first three examples in order to focus on the interaction behaviour and the performance of the crossing and diversion criterion. In the first example, we illustrate the performance of the criterion, and we investigate the influence of the shear stiffness parameter in the second example. The third example consists of a comparison with an experiment and a theoretical crossing criterion. In the fourth example, we do include pore pressure and show a propagating hydraulic fracture coming in contact with a small fracture network.

Performance of the crossing and diversion criterion
To illustrate the performance of the crossing and diversion criterion, we consider a hydraulic fracture growing perpendicular into a natural fracture (Fig. 8). The hydraulic fracture is propagating from a circular hole with a radius of 7 mm. in the integration points. In the right image the criterion is violated in the right element. Therefore, the ELP degree of freedom is extended to the right We demonstrate the effect of the crossing and diversion implementation by comparing two simulations. In the first simulation, we include both criteria. In the second, we prohibit crossing and fracture growth after the hydraulic fracture merged with the natural fracture. Thus, the only path for fluid to distribute is to divert into the natural fracture. In Fig. 9, the deformed mesh after 300 s of fluid injection is shown with the pressure in the fracture as contour. It is evident that with the crossing criterion the fracture propagated through the natural fracture as expected (Fig.  9a). Without the crossing criterion the fluid diverted into the natural fracture (Fig. 9b). The pressure in the fracture with crossing included is lower than the pressure when the fluid is forced to go into the natural fracture. This is a consequence of the natural fracture being perpendicular to σ v . Injection pressure over time shows similar behaviour, as shown in Fig. 10. In the case without the crossing mechanism, this leads to improbable magnitudes of injection pressure since fracture crossing would be energetically more favourable. The first drop in pressure (after ± 90 s) in the case without crossing is caused by fluid diverting mainly in the left wing of the natural fracture. Once the left wing is filled with fluid the pressure has to slightly build up again, after which the right wing is also completely filled with fluid (after ± 180 s). The preference of which wing will fill first is in this case a numerical artefact. Even though the problem is symmetric, the FE mesh is unstructured. As a result, one of the two sides opens first.

Influence of the shear stiffness
The magnitude of the shear stiffness is one of the key parameters in transferring stress from the hydraulic fracture across the natural fracture which can lead to fracture crossing. The magnitude of stress that is transferred is finite when described by the Coulomb friction law. We use a simplified penalty stiffness law to describe the friction, as explained in Section 2.3. The penalty law leads to an infinite friction which increases linearly with the shear displacement according to stiffness parameter D s . In this example, we show the effect of this shear stiffness on the fracture propagation. We consider the same specimen and material properties as in Example 3.1. The only parameters that are varied are stiffness values of the natural fracture. Also the natural fracture is placed further away from the hydraulic fracture at distance of five centimetre in y-direction measured from the centre of the circle.
To determine the effect of the shear stress transfer, the shear stiffness is varied between 0.0 and 20.0 GPa. A penalty stiffness of 10.0 GPa is used. The length of the hydraulic fracture, including the cohesive zone, is plotted against the time in Fig. 11a. We observe almost identical fracture growth during the first stage of hydraulic fracturing. The natural fracture is too far away to influence the propagation. In the vicinity of the natural fracture, approximately 10.0 mm, there is a minor influence (Fig. 11b). There is some strength loss due to the imperfect natural fracture, leading to accelerated fracture propagation with lower shear stiffness.  There are two possibilities for the fracture to grow further after the merging with the natural fracture. The fluid can divert into the natural fracture and eventually grow the natural fracture or the hydraulic fracture must cross the natural fracture. We did not observe fluid diverting into the natural fracture. Fluid diversion is unfavourable since the maximum confining stress must be overcome to open the natural fracture. In Fig. 11c, it is shown that after merging, the pressure is building up in the hydraulic fracture. This leads to stress Fig. 12 Scheme for the interaction example. Fluid is injected in the centre of the hydraulic fracture. The centre of the natural fracture is, independent of the interaction angle, located 45 mm above the hydraulic fracture transfer across the natural fracture and eventually to fracture nucleation. The speed of stress transfer depends on the magnitude of the shear stiffness (Fig. 11d). Only with zero shear stiffness we did not observe fracture crossing.

Influence of the interaction angle and the in situ stress
It is known that the ratio between in situ stresses and the angle between a hydraulic fracture and a natural fracture are important to predict whether the hydraulic fracture crosses, diverts or arrests on the natural fracture. This relation is shown in experimental results [5,6] and is described with a crossing criterion by Gu et al. [8]. In this example, we compare our numerical results with the crossing criterion. We vary the in situ stress by varying σ h while keeping σ v equal to -20 MPa. The orientation of the natural fracture is also varied between θ i = 90 • and θ i = 15 • (see Fig. 12). The remainder of the material properties are identical to those used in Example 3.1.
To interpret the results we distinguish between four different length measurements, i.e. the length of the hydraulic fracture, the length the natural fracture, the part of the hydraulic fracture that has crossed the natural fracture, and the part of the natural fracture that is filled with fluid. These results are shown for each interaction angle measured in Fig. 14. It can be seen that high interaction angles show mainly crossing of the hydraulic fracture. In some cases, the length of the crossed hydraulic fracture (shown in the blue bar) is small. We attribute this to the loss of friction due to the crossed fracture, which leads to relaxation of strain tangential to the natural fracture. Since a penalty friction law is used, there is a lower stress transfer across the natural fracture. The bottom part of the hydraulic fracture then becomes the favourable growth direction. The results of the low interaction angles show diversion of fluid into the natural fracture. Only the right wing of the natural fracture, which is not feeling pressure exerted by the hydraulic fracture, that is filled with fluid (Fig. 13). In none of our results the natural fracture propagated. The bottom part of the hydraulic fracture propagates as soon as the right wing is filled with fluid. Interaction angles varying between θ i = 50 • and θ i = 60 • show crack arrest.
In Fig. 15, we have extracted from Fig. 14 whether crossing, arrest or diversion of the fracture occurred for the various interaction angles and stress differences simulated. Apart from one outlier, we see a trend from fracture crossing to arrest and finally diversion with decreasing interaction angle. This tendency is also observed in the experiments. The influence of the stress difference is less pronounced. This is attributed to the friction law which does not possess the same properties as a Coulomb friction, which would better represent friction in rocks. There is no explicit relation between the magnitude of the friction and the normal stress to the natural fracture in our penalty friction law.
The tendency of crossing or diversion for various friction coefficients based on the crossing criterion from Gu et al. [8] is also shown in Fig. 15. The region right to the curve represents crossing. The left region indicated diversion except close to the line, where crack arrest is expected. The criterion does not distinguish between diversion or arrest. Our results are in the proximity of μ = 2.0 and μ = 1.5. Rocks typically do not possess such a high friction coefficient [8]. It is likely that we overestimate the magnitude of the friction due to our penalty friction law. However, the trend from diversion to crack arrest and finally crossing as the interaction increases is visible in our numerical results.

Interaction with a fracture network
In this final example, we consider a hydraulic fracture propagating towards a small network consisting of three natural fractures that are in contact (Fig. 16). The two larger natural fractures have a length of 40 mm and make an angle of 30 • with the x-axis. The smaller natural fracture is perpendicular to both of the larger fractures. The specimen is now considered to be a low permeable rock with an intrinsic permeability of K int = 10 −21 m 2 . The solid grains have a compressibility of K s = 30 GPa and the fluid compressibility is taken The contour plot of the pore pressure in the bulk for four different time-steps is shown in Fig. 17. In Fig. 17a, the hydraulic fracture propagated and a low pressure in front of the cohesive zone can be seen. This is a result of tension due the pressurisation of the hydraulic fracture and leads to fluid being attracted towards the region under tension. After 8 s of fluid injection, the hydraulic fracture crossed the bottom natural fracture (Fig. 17b). The angle between the middle natural fracture and the direction of σ v is such that crossing is not likely to occur. This is also observed in Fig. 17c. Fluid diverted from the hydraulic fracture into the natural fracture network. Fluid diversion stops when the fluid front reaches the top natural fracture. Tension is being generated at the top surface due to the friction law. This can also be seen at the low pressure in Fig. 17c. Finally, a new cohesive zone nucleates and the hydraulic fracture propagates away from the natural fracture network (Fig. 17d).
A discontinuous injection pressure development over time is shown in Fig. 18. On the one hand, this is caused by Fig. 15 The numerical results whether fluid diversion, crack arrest or fracture crossing occurs due to varying interaction angles and stress differences are shown with the point markers. The crossing criterion from Gu et al. [8] for various friction coefficients is shown with the lines. The region right of a line indicates crossing

Conclusions
Based on the enhanced local pressure model, we present a numerical model to investigate the interaction of multiple cracks in hydraulic fracturing. There is no limit to the amount of fractures. Interaction and propagation can take place in arbitrary locations and directions. The criterion whether or not a fracture crosses a natural fracture is determined by an average stress criterion. Fluid is also given the possibility to divert into a natural fracture based on a criterion that examines whether the natural fracture opens. These two criteria are checked simultaneously and the competition between them determines the interaction behaviour. The fracture process is modelled by a cohesive zone model. With the ELP model the injected fluid flow goes exclusively in the fracture and steep pressure gradients near the fracture surface are reconstructed based on Terzaghi's consolidation solution.
Four examples are presented to illustrate the implementation of the criteria and to show the performance of the numerical model. The first example shows the effect of the crossing and diversion implementation. A hydraulic fracture propagated perpendicular onto the natural fracture with the in situ stress taken such that crossing of the natural fracture is preferred. As expected, the hydraulic fracture indeed crosses the natural fracture. The simulation was repeated but now the crossing criterion is not included in the model. This leads to fluid diversion into the natural fracture and has as a consequence that the required injection pressure is much higher. In the second example, the effect of the implemented friction law is shown. The penalty friction leads to an imperfect natural fracture. With a higher stiffness, the interaction took place earlier and also the fracture crosses the natural fracture at an earlier time requiring less injection pressure. The behaviour observed in these two examples is consistent with what would be expected from the chosen geometry and the boundary conditions. The influence of the interaction angle and the in situ stress conditions is studied in the third example. These results are compared with experimental data and with a theoretical crossing criterion. A trend from fracture crossing to arrest and fluid diversion is observed with a decreasing interaction angle. This is in line with the experiments and the theoretical criterion. In the last example, poro elasticity is included and a hydraulic fracture interacting with a small natural fracture network is studied. Due to tension, a low pore pressure is observed in front of the cohesive zone. Fracture crossing and fluid diversion are both shown depending which natural fracture was interacting with the hydraulic fracture.
The proposed model is suitable to simulate complex hydraulic fracture patterns in fully saturated porous media. The simplified friction law leads to some discrepancy but the global fracture behaviour satisfies expectations and experimental results. Slip behaviour is not considered but could be included by replacing the friction law with a plasticity based friction law. This and considering shear failure is left for future research.