Impact of geomechanical heterogeneity on multiple hydraulic fracture propagation

To extract more gas from shale gas reservoirs, the spacing among hydraulic fractures should be made smaller, resulting in a significant stress shadow effect. Most studies regarding the stress shadow effect are based on the assumption of homogeneity in rock properties. However, strong heterogeneity has been observed in shale reservoirs, and the results obtained with homogeneous models can be different from practical situations. A series of case studies have been conducted in this work to understand the effects of mechanical heterogeneity on multiple fracture propagation. Fracture propagation was simulated using the extended finite element method. A sequential Gaussian simulation was performed to generate a heterogeneous distribution of geomechanical properties. According to the simulation results, the difficulty of fracture propagation is negatively correlated with the Young’s modulus and Poisson’s ratio, and positively correlated with tensile strength. When each of the multiple fractures propagates in a homogeneous area with different mechanical properties, the final geometry of the fracture is similar to homogeneous conditions. When the rock parameter is a random field or heterogeneity perpendicular to the propagation direction of fracture, the fracture will no longer take the wellbore as the center of symmetry. Based on the analysis of fracture propagation in random fields, a small variance of elastic parameters can result in asymmetrical propagation of multiple fractures. Moreover, the asymmetrical propagation of hydraulic fractures is more sensitive to the heterogeneity of Poisson’s ratio than Young’s modulus. This study emphasises the importance of considering geomechanical heterogeneity and provides some meaningful suggestions regarding hydraulic fracturing designs.


Introduction
With the development of unconventional energy sources, especially shale gas, oil and gas, companies/professionals have become increasingly interested in hydraulic fracturing, which is widely used to exploit oil and natural gas from shale (Rahm 2011;Gandossi & Von Estorff 2013;Vengosh et al. 2014;Zhang et al. 2016). Simultaneously and sequentially multi-stage hydraulic fracturing is a commonly used technique. Usually, multiple perforations are arranged in a fracturing stage and fractured simultaneously to increase the contact area between the fractures and reservoir (Miller et al. 2011;Tang et al. 2019). Understanding the dominate factors controlling the propagation of multiple hydraulic fractures is critical for the design and optimisation of hydraulic fracturing treatments.
From the analytical solutions of PKN and KGD models (Valko & Economides 1995), the influencing factors include fracturing fluid viscosity, injection rate, in situ stress, rock elastic parameters, etc. Like single fracture propagation, from the results of numerical simulation and experiments, the fracture network is also affected by the aforementioned factors and additional stress shadow (Olson 2008). The distance between perforations is so small that it induces a stress shadow effect of closed spaced fractures. The stress shadow affects parallel fractures significantly, including (i) reduction in fracture length; (ii) increase in fracture height; (iii) change in fracture width and (iv) hydraulic fractures deflect from the direction of maximum principal stress (Nagel et al. 2017;Chen et al. 2019;Kolawole et al. 2019).
In a homogeneous layer, when multiple fractures propagate, this will cause uneven fracture geometries (Li et al. 2017). To balance the propagation of fractures within one stage, Wu et al. (2017) found through numerical simulation that the fracturing fluid can be better distributed by controlling the diameter of the outer perforations and the spacing between perforations, to make the length of the fractures more uniform. In the case of the same fracturing fluid flow rate for each perforation, Dontsov & Suarez-Rivera (2020) examined the influence of different factors on the fracture geometry. In the case of a fracturing fluid with high viscosity, each fracture presents a radially symmetrical geometry. In contrast, with a low viscosity fracturing fluid or high fracture toughness of the rock, the fractures will often interact strongly, be divided and form a geometric structure similar to the petals of a flower. This influence tends to be weakened with increasing fracture spacing. With the increase in perforation amount and the influence of friction and stress shadow, the final fracture network consists of 2-3 dominated perforated hydraulic fractures and the remaining micro hydraulic fractures . In addition to the stress shadow between the fractures of a single horizontal well, Zheng et al. (2020) analysed the stress shadow between different horizontal wells. The results show that zipper fracturing significantly reduces well interference on the dynamic fracture propagation compared to parallel fracturing.
Most of the mechanical parameters in these models for the study of stress shadow and fracture network geometry are assumed to be homogeneous or of layered heterogeneity. Only a few of the mechanical parameters in the studies were completely heterogeneous. However, there is strong heterogeneity in shale reservoirs (Kumar et al. 2012;Chen et al. 2015). Mokhtari et al. (2016) analysed the heterogeneity and anisotropy of the shale of Eagle Ford using methods such as scanning-electron-microscopy, computed tomography scanning and compressional-velocity scanning. These results showed that the uneven distribution of minerals, pores and biological debris during the deposition process resulted in the heterogeneity and anisotropy of shale. Zhang et al. (2019) used X-ray fluorescence and cuttings from a horizontal borehole to study heterogeneity. Horizontal heterogeneity was confirmed and manifested in the different distribution of minerals and organic matter in the horizontal direction. The heterogeneity of the grain size leads to obvious changes in the stress shadow and a significant difference in the length of multiple fractures (Duan et al. 2020). When Young's modulus is heterogeneously distributed, the fractures are the most irregular and asymmetric. If the heterogeneity of other mechanical parameters is considered. the irregularity and asymmetry of fractures increase . The fracture propagation in homogeneous rock is faster than in heterogeneous. At the same time, the difference in fracture height in heterogeneous rock is greater than the homogeneous condition. The propagation pore pressure was higher for the later injection stage in the homogenous rock material than the heterogeneous . The high injection rate could help neutralise the effects of heterogeneity (Tang et al. 2018).
With the development of computer technology and mathematics, there are many methods for simulating multifracture propagation, such as the finite difference method, discrete element method (Fatahi et al. 2017;Nagaso et al. 2019), phase field method (Ni et al. 2020;Zhuang et al. 2020), displacement discontinuity method (DDM) (Tang et al. 2016) and finite element method (FEM) (Chen 2012). The DDM can construct the elastic equations easily since the stresses applied on the fracture elements equal the stresses induced by all displacement discontinuities plus the in situ stresses. Although the DDM has been widely accepted as an efficient and precise method for fracturing simulation, the derivation of governing equations of DDM to consider poroelasticity or elastic heterogeneity is challenging. Because the discrete element method regards the rock as a collection of bonded particles, for large-scale hydraulic fracturing design, a great number of elements will be required. At the same time, it is challenging to calibrate the parameters in the micro-constitutive models. The principle of the phase field method for simulating hydraulic fracturing is to add a continuous scalar variable to describe the smooth transition between non-fractured and completely fractured materials. Usually, the displacement discontinuities associated with the fractures on a small number of grids are approximate, that is, they are not exact. The advantage of this method is that it can simulate complex fracture geometries, stress shadow effects, heterogeneity, etc. But it is difficult to accurately simulate the displacements of fracture surfaces. The extended finite element method (XFEM) was a new method that has been used widely in various engineering simulations and hydraulic 955 fracturing in recent years (Chen 2013;Sheng et al. 2018). Compared to the cohesive method (Haddad & Sepehrnoori 2016) and the FEM where fractures must propagate along the grid boundary, XFEM does not need to remesh the grid through which the fractures may pass. Hence, the calculation takes considerably less time. However, it is challenging to simulate 3D surface fractures with complex geometries, especially fracture intersections, using XFEM due to the difficulties of constructing enrichment functions at fracture tips.
The numerical method we used for fracture propagation is XFEM, and to realise a random field of elastic properties in ABAQUS, a sequence of operations is required, including homogeneous case construction in ABAQUS, random field generation using GISLIB package (Ziegel et al. 1998), assignment of material properties to every grid in ABAQUS using an in-house developed Python code and running the heterogeneous ABAQUS model with XFEM. In this paper, we first introduce the random field generation method and the basic principle of XFEM. The fracture geometries of a single fracture and multiple fractures in a heterogeneous field are compared. Finally, the influence of random heterogeneous geomechanical parameter distribution on fracture propagation and fracture interference are investigated.

Generation of heterogeneous fields
The first step to realise multi-fracture propagation in heterogeneous material is random field generation. Initially, according to the literature (Rybacki et al. 2015;Pakzad et al. 2018), the parameters affecting the propagation of multiple hydraulic fractures most include: Young's modulus (E), Poisson's ratio ( ), tensile strength (T 0 ), porosity ( ) and permeability (k). Among them, the mechanical parameters (i.e. Young's modulus and Poisson's ratio) are chosen in our work to investigate the impact of heterogeneity on fracture propagation. The reason we only consider the random field of elastic properties is that in practice, the field of Young's modulus and Poisson's ratio can be obtained more directly from the logging data such as density, V p and V s (P and S wave velocities). Other mechanical parameters, such as tensile stress, are often calculated with empirical functions with less accuracy. The choice of elastic properties is based on the higher chance to obtain their spatial distribution in practice. The porosity, permeability and other parameters are constant and homogeneous in the following cases.
The SGSIM software in the GISLIB software package can generate geomechanical random fields using the sequential Gaussian simulation method (SGSIM) (Ziegel et al. 1998). After inputting the values of the required parameters, such as domain size and variogram type, SGSIM can generate random field data with the specified correlation length, variance and other parameters. Then, the random field of Young's modulus and Poisson's ratio can be obtained.

Principle of XFEM
The XFEM is a new method that overcomes the shortcomings of the traditional FEM. It can model any discontinuities and display the fracture propagation path without remeshing. Belytschko & Black (1999) proposed a new method that required the re-division of a smaller grid to solve the fracture propagation problem in FEM. In the following year, Daux et al. (2000) developed an XFEM by adding specific extended functions and additional nodal degrees of freedom to standard FEM.
Level set function is an essential part of the XFEM. Osher & Sethian (1988) proposed the level set method (LSM) to capture propagating fractures. The LSM uses a level set function to express the low-dimensional changes of the interface as high-dimensional changes to separate the discontinuous description from the finite element grid.
In the extended finite element simulation, the grid is divided into three types. The first type of grid has no fractures and fracture tips, which can be considered a standard finite element grid. The second type is a grid through which a fracture passes. When processing this type of grid, it is necessary to introduce step functions as enrichment functions based on the finite element.
where x o is the grid node closest to the point x in the calculation grid area and dist() is the vector distance from the node to the fracture.
The last type of grid has fracture tips. The following four equations are necessary for enrichment functions, as shown next: (2) in which (r, ) is the local polar coordinate system of the fracture tips. The specific form can be expressed as: where tip is the angle of the fracture tip in the orthogonal coordinate system. 956

Model definition
Fluid flow in ABAQUS can be divided into the following three types: (i) flow in fractures, (ii) flow in the matrix, (iii) leak-off from fractures to matrix, where w is the fracture width, m; q is flow rate, m/s; is the viscosity, Pa⋅s; k is permeability, m 2 ; c l is the dimensionless coefficient of leak-off and ∇p is the pressure difference, Pa.
According to the principle of virtual work, if the object is in equilibrium under the action of an external force and the object has a slight virtual displacement, the virtual work of an external force, i.e. volume force (f ) and surface force (t), must be equal to the virtual work of the stress (∫ V : dV).
In the equation (7), = v x , = ′ + p. v is the virtual velocity field; p is the pore pressure; is the virtual rate of deformation; ′ is the effective stress and is the total stress.
The continuity equation for fluid phase is: where n is the free liquid; n t is the trapped liquid; is the fluid density.
The maximum principal stress criterion is used as the failure criterion (ABAQUS 2017): o max is the maximum allowable stress: When the value of equation (9) evolves to 1, the material will be damaged.
Viscous regularisation is a default operation in ABAQUS. Softening behavior and stiffness degradation in material will lead to convergence difficulties. After using viscous regularisation, the tangent stiffness matrix of the softening material will be positive in a small time step to overcome the convergence difficulty.

Model validation
We compare the results of our model with the analytical solution of the KGD model (Valko & Economides 1995). The parameters used are listed in Table 1. The results calculated with our model and analytical solutions are compared in figure 1. As shown in figure 1, our model results agree well with the KGD model.
Besides comparing with the theoretical model, we further compare our model with the laboratory data ). This experiment is a true triaxial hydraulic fracturing experiment conducted on a sandstone sample of 300 mm in length, width and height. First, a hole with diameter 25 mm is drilled in the sample, and then 10 mm fractures are generated on the opposite side of the drilled hole. The parameters used in this experiment are shown in Table 2. The fracture path predicted by our model is compared with the experimental results of sample S-3 in the literature.
The laboratory results and our simulation results are shown in figure 2. The dotted red line in the figure is our simulation result and the light black line is the experimental result. A good agreement can be observed.

Single fracture propagation
The influence of formation heterogeneity on fracture propagation and inter-fracture interference is investigated with the following cases. Table 3 lists the parameters of these cases.
Cases 1.0-1.3 are the propagation of a single fracture under regional heterogeneity. The 50 × 50 m 2 area is divided into different areas, and different rock properties were assigned to reveal the influence of rock-mechanical properties on fracture propagation qualitatively. In these cases, the number of grids is 10 000 and the length of the grid is 0.5 m. Cases 1.0-1.3 consider the heterogeneity of three different rock-mechanical properties. Cases 2.0-2.4 consider the propagation of multiple hydraulic fractures. In these 957 Figure 1. Comparison of (a) fracture width and (b) half-length between our results and the KGD model.  cases, the influence of heterogeneity perpendicular and parallel to the direction of fracture propagation is studied. Case 3.1 is a simulation of multi-fracture propagation in a random heterogeneous field with a spatial correlation. In this case, geostatistics software (SGSIM) is used to generate the cor-responding random field. Hence, the spatially correlated random field is used for research instead of a completely random field. Case 1 was designed to explore the influence of different rock-mechanical parameters. Figure 3a presents a schematic of model settings and boundary conditions. The principal stresses, 11, 22 and 33, are 6, 8 and 15 MPa, respectively. For comparison, the simulations were performed in a homogeneous and heterogeneous matrix. Table 4 lists the parameters used for simulation. The fluid total injection time is 70 s. The viscosity regularisation coefficient (0.001) is used in this section as numerical damping to increase the convergence of models.
The models are divided in to three regions with height ratio of 2:1:2, and each region is assigned with different materials. The fracture will pass through at least two of the regions before reaching the boundaries. The relationship between the fracture geometry and the rockmechanical parameters can be obtained by measuring the fracture lengths of different wings initiating from the perforations.
A higher Young's modulus assists in fracture propagation (Table 5 and figure 4). Similarly, a larger Poisson's ratio is also helpful to the initiation of fractures. In Case 1.3, the area with a lower tensile strength is easier to damage because the rock with lower tensile strength is more likely to yield. The development of the fracture requires the generation of sufficient tensile stress in the unbroken element connected to the fracture tip to achieve tensile strength. Since the material with a larger Young's modulus and Poisson's ratio has a stronger ability to resist deformation and can transfer enough stress to the next element with a smaller deformation, these areas are more likely to be damaged. Figure 5 shows the results of the fracture length and width of Cases 1.0-1.3. In heterogeneous rock, the trend of the fracture width and length of Cases 1.0-1.3 are similar, and the geometry of the fractures is basically the same. Compared with 958    the homogeneous rock, the fracture width at the perforation is no longer the largest.

Multiple fractures propagation
Based on the previous cases, a multi-fracture propagation model with three perforations is designed, as shown in figure 3b. Table 6 lists the case settings. The spacing between perforations is 12.5 m, and the injection rate of each perforation is the same. In practice, due to the difference of fracture pressure at the perforation, the flow rate among different perforations within one stage can be different   (Kumar & Ghassemi 2016;Tang et al. 2019). However, the flow within the wellbore has not been simulated in our current model, thus the unbalanced distribution of injected fluid is not considered. Such even distribution of injection fluid might be realised when extremely limited entry technology is used. We are still working on a new UDF in ABAQUS to consider the fluid flow within the borehole. The fracture geometries obtained by our model under homogeneity assumptions are similar to Haddad & Sepehrnoori's (2016) work: they also assumed the injection fluid is evenly distributed in all perforations. The model partitions are the same as Cases 1.1-1.3, and each partition was assigned a different Young's modulus to investigate multi-fracture propagation in a heterogeneous area. The results of other mechanical parameters can be inferred from the Young's modulus and the conclusion of Case 1. Figure 6 presents the result of Case 2.0. In ABAQUS, positive stress represents tensile stress and negative stress represents compressive stress.
The fracture width at middle perforation is the smallest, but it has greater length. Because the fluid in the fractures on both sides can release energy through matrix deformation, i.e. squeeze the middle fracture, the length of the fractures on both sides is not long, but the width of the fracture at the perforation is the largest. In practice, it is more difficult for internal fractures to propagate. Moreover, it is difficult for the fluid in the wellbore to enter the internal fractures, which hinders the propagation of internal fractures (Tang et al. 2016;Rao & Wang 2020).  When approaching the terminal of fracturing, the length of the middle fracture suddenly increases, which is caused by the squeezing of the fractures on both sides.
With continuous injection of fracturing fluid into perforations, length of the middle fracture eventually exceeds the other fractures.
Cases 2.1 and 2.2 are designed to further understand the impact of Young's modulus heterogeneity on multiple fractures propagation. In Case 2.1, the Young's modulus from top to bottom is 32, 35 and 37 GPa, respectively. In Case 2.2, the counterpart is 35, 32 and 37 GPa. Except for the different Young's moduli, other parameters are the same as Case 2.0. Figures 7 and 8 demonstrate the results of Cases 2.1 and 2.2.
The conclusions of Case 1 suggest that fractures prefer to propagate toward the area with a higher Young's modulus. However, the middle fracture in Cases 2.1 and 2.2 chose the opposite direction. In contrast to the middle fracture, fractures on the two sides propagate mainly toward the region with a high Young's modulus.
Due to the squeezing of the outer fractures on the inner fractures, the side fractures grow faster and reach the boundary of the middle area earlier. Since the fractures are prone to propagate to the area with high Young's modulus, the side fractures grow to the high Young's modulus area instead of moving upward. The possibility of the inner fracture to propagate to the high Young's modulus area is suppressed, and the inner fracture can only grow in the opposite direction.  In a heterogeneous matrix, the geometry of multifractures will become more complicated. According to the simulation results in Cases 2.1 and 2.2, if there is heterogeneity perpendicular to the direction of fracture propagation, the fractures on both sides will propagate preferentially to the area with the larger Young's moduli, while opposite is observed with middle fracture. In the two fracture tips of a perforation, there is always one fracture tip that propagates farther than the other fracture tip in a perforation caused by the stress shadow and heterogeneity. The stress shadow shortens the propagation distance of the fracture tip, and the heterogeneity determines the main propagation direction of the fracture tip.
Cases 2.3 and 2.4 are designed to understand the influence of the heterogeneity parallel to the propagation direction on the fracture geometry. Except for the different Young's mod-uli, the parameters of the fracture propagation model in Cases 2.3 and 2.4 are the same as those in Case 2.0. Figure 9 presents the simulation results of Cases 2.3 and 2.4. The final geometries of the three fractures are similar to 0.0. Judging from the length of the fractures on both sides, the conclusion from Case 1 can be verified, i.e. the fracture is easier to propagate in the area with larger Young's modulus. Furthermore, the longer fracture has a smaller width because the flow rate injected into each fracture is constant.
From Cases 2.1-2.4, the influence of heterogeneity in two directions on fracture geometries is demonstrated. If there is heterogeneity perpendicular to the direction of fracture propagation, it will cause all fractures to be no longer symmetrical along the wellbore. There is a main propagation direction for each fracture, and the midpoint of the fracture and the area with maximum width will no longer be the 962 Figure 11. Propagation process of each fracture in Case 3.1: (a) left, (b) middle and (c) right. perforation. The other case is that the heterogeneity is parallel to the direction of fracture propagation. This is equivalent to the three fractures propagating in a different matrix. The difference is that the boundaries of the three rectangles are interconnected, and there are stress, strain and fluid exchanges. Therefore, the homogeneity conclusion will cause a larger error when heterogeneity exists in the rock perpendicular to the fracture propagation direction.

Fracture propagation in a random field
In this case, based on existing geological information (Mokhtari et al. 2016), the SGSIM (Ziegel et al. 1998) is used to generate a Young's modulus random field with a mean value of 37.5 GPa. The other parameters remain the same as Case 2.0.
The data generated by this algorithm is constrained by a distribution function, variogram, correlation coefficient and search radius. The stochastic parameters for the random field in Case 3.1 are as follows: the search radius in all three directions is 1000 m; the correlation lengths in the x and y directions are 100 m. Figure 10a shows the generated random field, and figure 10b shows the statistical distribution of the Young's modulus. In-house developed Python code is used to assign the generated random data to different grid in the XFEM model in ABAQUS.
The simulation results are presented in figures 11 and 12. The data of the three fractures in figure 12a shows that the widths are similar, but the propagation directions  differ significantly among the three fractures. Figure 11 parts a, b and c present the evolution of fracture geometries over time. First, the fracture is symmetrical along the center of the perforation (y = 0). The fracture tips start to propagate in different directions with the injection of the fracturing fluid. At the same time, position of maximum width moves gradually to the main propagation direction and deviates from the central axis. When the material properties of each grid were no longer the same, it required more time steps to balance the creep strain increment and the total elastic strain in ABAQUS. Therefore, the efficiency of frac-ture propagation under random field conditions is lower than Cases 2.0-2.4.
In hydraulic fracturing, the width of the fracture is also very important. The fracture must have a specific width to ensure that the proppant can enter the fracture to increase the conductivity. The width of the perforation varied with time, as shown in figure 12c. The fracture width at the perforation will decrease slowly after reaching the peak because the fractures squeeze the fracture in the middle from both sides. Figure 12d shows that all the fractures are similar in length. The length of the fracture changes almost linearly with time. 964 Figure 15. Young's modulus (with unit GPa) distribution with different variances (vr, unit: GPa 2 ) and correlation lengths (CL, unit: m).
With the same stochastic parameters, multiple realisations can be generated. To decrease the uncertainties of results, more results based on other realisations are also conducted. The other two generations of Young's modulus random fields and the corresponding results are shown in figure 13.
Although the inner fracture propagates in opposite direction, asymmetric distribution of fracture length can also be observed in different realisations with similar rock displacements.
The impact of the heterogeneity of Poisson's ratio on multiply fracture propagation is also investigated. The correlation length of the random field of Poisson's ratio is the same as the Young's modulus, but the variance is 0.0018. The distribution of Poisson's ratio and the statistical data are shown in figure 14a and b.
The heterogeneity Poisson's ratio will cause asymmetrical fracture propagation, which is similar to the Young's modulus. Also, the fractures prefer to going into the area with higher Poisson's ratio, which also agrees with our previous analysis.
In these cases, the stochastic parameters of random field are fixed. Here, we further analyse the role of random field parameters on fracture geometries. The random fields with different parameters are shown in figure 15.
The calculated fracture geometries are shown in figure 16. Similar to Case 3.1, the inner fracture propagates toward the opposite direction compared to the side fractures. Increasing correlation length (figure 15g-l) can lead to greater compressional stress in the central area. Besides, the two outer side fractures present more similar geometries due to 965 Figure 16. Fracture geometry with different variances and correlation lengths. the increase of spatial continuity. With the increase of variance, the difference between two side fractures becomes more significant (figure 15a-f). Meanwhile, taking the ratio of the longest fracture length to the shortest fracture length at the end of fracturing as the y axis, the influence of heterogeneous field parameters generated by different initial random number seeds on fracture properties is analysed. As shown in the figure 17, different colours of points represent different realisations of random fields. In figure 17a, the correlation length is 100 m for all cases. The input variance is fixed to be 10 GPa 2 in figure 17b. With the increase of variance or the decrease of correlation length, the ratio increases gradually. The reason is that with these changes, the difference of Young's modulus in different fracture areas increases.
On the basis of Case 2.1, we further conduct several case studies with the elastic properties in the upper and middle areas to be constant, and the lower part changes. The Young's modulus in these two parts is 35 GPa and the Poisson's ratio is 0.25. The difference of Young's modulus between these two parts and the lower part decreases 966 Figure 17. The ratio of the longest fracture length to the shortest fracture length with (a) different variances (unit: GPa 2 ) and (b) correlation lengths (unit: m) with multiple realisations. exponentially from large to small to find the threshold of symmetrical propagation.
It can be seen from figure 18 that the threshold of Young's modulus and Poisson's ratio are 0.03 GPa and 8e-5 in the partitions. In random field, the threshold standard deviation of Young's modulus and Poisson's ratio are 0.05 and 1e-5, respectively. The results obtained from the partition cases are similar to the random field cases. For the cases discussed here, a small variance of elastic parameters can result in asymmetrical propagation of multiple fractures, which means the assumption of homogeneity can only be correct for very specific formations. In addition, the asymmetrical propagation is more sensitive to the heterogeneity of Poisson's ratio (the coefficients of variation: standard deviation/mean value) for Young's modulus and Poisson's ratio are around 0.001 and 0.0002, respectively).

Conclusions
In this paper, the fracture geometries and fracture interference in a heterogeneously mechanical field were studied 967 using the XFEM. The main conclusions of this paper are summarised as follows.
(i) Using the XFEM and assigning different material properties to different grids in Python, the simulation of simultaneous multiple fractures propagation under heterogeneous conditions can be realised. (ii) In the case of single fracture propagation, the fracture is more likely to propagate in the region with larger Young's modulus, larger Poisson's ratio and smaller tensile strength. (iii) When the rock heterogeneity is parallel to the fracture propagation direction, the final geometries of fracture are similar to the results under homogeneous conditions. Therefore, when heterogeneity only exists in this direction, the simulation results under homogeneous conditions are still applicable. (iv) When the three fractures propagate simultaneously, and there is heterogeneity in the direction perpendicular to the fracture propagation, the two wings of the fracture lose symmetry. Hence, the homogeneity assumption may lead to an incorrect estimation of the fracture geometry. (v) When heterogeneity of the formation is strong, compared to the conclusions drawn from the homogeneous models, the length of different fractures is more uniform. Therefore, in some areas, more fractures might be able to propagate simultaneously than predicted by the homogeneous models. (vi) The larger variance and smaller correlation length of geomechanical properties could result in greater asymmetrical propagation of fractures in different clusters. Only a small variance of elastic properties can result in asymmetrical fracture geometries. Compared to the Young's modulus, the asymmetry growth of fractures is more sensitive to the variance of Poisson's ratio.