Investigation of Hydraulic Fracturing Behavior in Heterogeneous Laminated Rock Using a Micromechanics-Based Numerical Approach

: Understanding hydraulic fracturing mechanisms in heterogeneous laminated rocks is important for designing and optimizing well production, as well as for predicting shale gas production. In this study, a micromechanics-based numerical approach was used to understand the physical processes and underlying mechanisms of fracking for di ﬀ erent strata orientations, in-situ stresses, rock strengths, and injection parameters. The numerical experiments revealed a very strong inﬂuence of the pre-existing weakness planes on fracking. Geological models for rock without weakness planes and laminated rock behave very di ﬀ erently. Most simulated fractures in the rock without weakness planes were caused by tensile failure of the rock matrix. In an intact rock model, although a radial damage zone was generated around the injection hole, most of the small cracks were isolated, resulting in poor connectivity of the fracture network. For rock models with pre-existing weakness planes, tension and shear failure of these structural planes formed an oval-shaped network. The network was symmetrically developed around the injection well because the strength of the pre-existing weakness planes is generally lower than the rock matrix. The research shows that the angular relations between the orientation of the structural planes and the maximum horizontal stress, as well as the in-situ stress ratios, have signiﬁcant e ﬀ ects on the morphology and extent of the networks. The strength of the pre-existing weakness planes, their spacing, and the injection rate can dramatically inﬂuence the e ﬀ ectiveness of hydraulic fracturing treatments.


Introduction
Hydraulic fracturing is an important technology for shale gas developement [1]. This technology depends on multiple parameters, and fracture formation is dominated by several factors, e.g., the in-situ stress, the rock mass properties, and the injection rates [1,2]. The seepage pathway of fluid flow in a fractured rock mass is mainly controlled by the geometry, pattern, and heterogeneity of the hydraulic fracture network [3][4][5][6]. Therefore, the quality of the artificially modified fracture network is a critical factor in practical applications [7,8].
The anisotropic mechanical properties of rocks play an important role in the damage and stability of rock structures. Zhang et al. [9] revealed that the deformation of Callovo-Oxfordian claystone depends on the orientation of the major principal stress to the bedding and the claystone strength depends on the loading path and direction with respect to the bedding. Valente et al. [10] analyzed the mechanisms of parallel-to-bedding cracking through a sensitivity analysis and using a numerical paper because of its advantages and ability to simulate the interactions between induced hydraulic fractures and pre-existing fractures in heterogeneous rock.
In this study, the evolution of the hydraulic fracturing process is simulated in a model of anisotropic shale at the laboratory scale in which the shale microstructure was explicitly represented by inserting one or two sets of parallel discontinuities into the model. Through a series of modeling experiments, the influence of the discontinuities on the hydraulic fracturing was extensively investigated. In addition, hydrofracturing simulations were conducted to assess the influences of the in-situ stress ratio, the mechanical properties of the discontinuities, the discontinuity spacing, and the injection rate on the evolution of hydraulic fractures.

Modeling Methodology
The FSD coupled approach was derived from damage mechanics and statistical theory by considering the deformation of an elastic material containing a random initial distribution of micro-factures [35,[38][39][40][41][42]. This modeling approach can simulate nonlinear rock behaviors and macroscopic rock fractures without knowledge of where and how the fractures will occur. Compared with other numerical methods, there are two distinct features: (1) By introducing the heterogeneity of rock parameters, nonlinear behavior in rock can be simulated, and (2) by introducing elastic modulus reduction after element failure, discontinuum mechanics can be simulated in a continuum mechanics model. For a more detailed description of the method, readers are referred to the references.
To simulate the progressive failure of the heterogeneous, permeable rock mass, the FSD coupled model assumes [38][39][40][41][42] that (1) the flow of fluid in the model follows the Biot consolidation theory; (2) the model material is assumed to be elastic-brittle with residual strength and its failure can be described by elastic damage theory; (3) the Mohr-Coulomb criterion and the maximum tensile strength criterion can be used to define the tensile failure and shear failure of the rock; and (4) the permeability coefficient of the rock material is varied with the stress in an elastic state and increases dramatically according to a deformation-dependent law when the element progresses towards the failure stage. The continuum damage principle is applied and the fracture is presented as the damage zone. The basic equations used in the analysis are as follows: Equilibrium equation: Strain-displacement equation: Constitutive equation: Seepage equation: where σ ij is the total stress in the ij-plane, σ ij is the effective stress in the ij-plane, ρ is the density, X j is the body force in the j-th direction, u is the element displacement, ε ij is the strain in the ij-plane, ε v is the volumetric strain, p is the pore water pressure, λ is the Lame coefficient, δ ij is the Kronecker constant, G is the shear modulus, k is the coefficient of permeability, α is the pore-fluid pressure coefficient (dimensionless), H and R(Q, α) represent the Biot's constant, β is a coupling parameter that reflects the influence of stress on the permeability coefficient, ξ is a damage factor of permeability (dimensionless), 1/R represents the water capacity change due to changes in water pressure, 1/H represents the change in the overall volume of the media due to changes in water pressure, and 1/Q represents the amount of water squeezed into the porous media under the action of water pressure without changing the volume of the porous media.
The relationship between the Q, R, H, and α is where E is the elastic modulus, K is the bulk modulus, and µ is the Poisson's ratio. In the process of steady flow, the flow tends to be stable as the pore water dissipates. Therefore, the value of Q is very large; assuming Q = ∞, Equation (7) can be simplified as A pore volume change of ∆n causes a flow rate change. The permeability coefficient k is a function of the pore volume change ∆n. Therefore, the coupling equation can be described as Equation (9) is a negative exponential function that is applied to describe the relationship between stress and damage. The most important hypothesis reflected in the FSD model is that a heterogeneous distribution of rock strength causes rock failure behavior. The rock medium is assumed to be locally heterogeneous by randomly assigning the Young's modulus and compressive strength to each element in a Weibull distribution [39]: where s is the mechanical parameter of the elements, such as elastic modulus or strength. The scale parameter s 0 is related to the average value of the element parameter, and the parameter m defines the homogeneity index; a higher value of m represents a more homogeneous material. The strength distribution with different m values is illustrated in Figure 1. In the FSD coupled model, the coupled effects of flow, stress, and damage on the fractures and the permeability of rock mass are considered [44]. Both the tensile and shear failure modes are considered in the model calculation. The maximum tensile strain criterion and the Mohr-Coulomb criterion are used to define the types of deformation and breakdown. According to the elastic damage mechanics, the elastic modulus of the damaged material can be defined as [39] where D represents the damage variable, and E and E 0 are the elastic modulus of the damaged and the undamaged elements, respectively.
When the tensile stress reaches its tensile strength, where f i is the tensile strength of the rock element. Then, the damage variable can be described as [39] where f tr is the residual tensile strength of the element, ε is the equivalent principal strain, ε to is the threshold strain, and ε tu is the maximum tensile strain of the element. Then, the permeability can be described as k(σ, p) = ξk 0 e β n = ξk 0 e β( p where ξ is a damage factor of permeability, k 0 is the initial coefficient, α is the pore-fluid pressure coefficient (dimensionless), β is a coupling parameter, and p is the pore water pressure.
In the FSD coupled model, tensile failure of the element occurs if the element's strength is smaller than the minor principal stress, and shear failure occurs when the compressive or shear stress satisfies the Mohr-Coulomb failure criterion: where σ 1 and σ 3 are the maximum and minimum principal effective stresses, φ is the effective angle of friction, and f c is the compressive failure strength. In this case, the damage factor under uniaxial compression can be described as where f cr is the residual compressive strength and ε cu is the maximum compressive strain. In this case, the permeability can be defined by

Model Setup for Hydraulic Fracturing Simulation
For this study, bedding joints and crossing joints were the structure planes considered in the modeling simulation. Figure 2 shows the basic geometry and model setup, which represents a 2D horizontal section of a reservoir. The model is composed of 250,000 (500 × 500) identical 2 mm square elements. The element size is small enough for the accuracy requirement of the model, and no further re-meshing was needed. A vertical wellbore with a radius of 20 mm is in the center of the model and serves as the injection hole during the simulation. The material mechanical parameters used in the FSD coupled model were calibrated through a series of trial-and-error experiments [44]. The microparameters were calibrated to match the macro properties of the Longmaxi shale rock, including the elastic modulus, peak strength, and Poisson's ratio. The mesoscopic physical and mechanical parameters of the model used in this study are shown in Table 1.

Fracture Distribution for the Base Model
The failure process and hydraulic fracture characteristics in heterogeneous rocks without the influence of bedding joints (base model) were studied first; the results serve as a base-case for comparison with other results presented in this paper. Figure 3 shows the progressive fracturing process of the model (in-situ stress ratio σ H / σ h = 1.5, injection rate = 0.15 m 3 ·d −1 ·m −1 ). This model is the base model for the following comparative study. As illustrated in Figure 3, the fracturing can be divided into five stages, depending on the fracture formation process: (1) Stress accumulation. A few micro-units clustered around the injection hole with a lower strength and stiffness are damaged ( Figure 3b). Since the tensile strength is far less than the compressive strength, and there is no preferential location around the injection hole for fractures to initiate, the location and orientation of the fracture initiation are unpredictable; (2) Small fracture formation. The continuous injection increases the pore pressure around the injection hole, and many small fractures occur in all directions around the injection hole (Figure 3c). Although a radial damage zone was generated, most of the small fractures were isolated with poor connectivity; (3) Hydraulic fracture initiation. With further injection, an embryonic network of connected hydraulic fractures gradually forms around the injection hole due to propagation and coalescence of the fractures ( Figure 3d); (4) Hydraulic fracture propagation. The dominant hydraulic fractures grow and extend parallel to the direction of σ H (Figure 3e). Meanwhile, many small hydraulic fractures occur around the main fractures as the injection continues; (5) Decay in fracture propagation. The rate of fracture propagation and the growth in the area with new fractures gradually decline because the constant injection rate cannot sustain the fracture growth far away from the injection hole. A large and growing proportion of the injection fluid is scattered into the fractures and the rock strata. The characteristics of fluid leak-off are related to fracture connectivity and the injection rate. In this example, both the size of the simulated area and the number of fractures achieved a peak and an oval-shaped fracture network eventually developed (Figure 3f).

Morphology of Fracture Networks in Laminated Rocks
The presence of structural weakness planes such as bedding in reservoir rock affects the response and effectiveness of hydraulic fracturing [7,8]. The pre-existing structural planes, such as bedding planes and tectonic joint sets, influence the volume and shape of the fractured zone in laminated shale formations [3,17]. In general, the strength of natural structural planes is lower than the strength of the rock matrix, and these structures are particularly vulnerable to damage and failure during hydraulic fracturing [17]. Figure 4 shows the modeled fracture networks under the influence of different pre-existing structural planes (bedding joints). In contrast with the base case model, the most remarkable feature is that the orientation of the structural planes exerts an important influence on the propagation direction of the hydraulic fractures. Under a constant injection rate, the structural planes were activated and damaged to form the framework of the hydraulic fracture network. In these models, the pore water pressure propagated throughout a large area of the model through these fractures. The results are in good agreement with most of the laboratory experiments and numerical results. For instance, Guo. et al. [7] studied the propagating rules of fractures in shale through a true triaxial test system and found that the hydraulic fractures easily propagate along with the natural fractures ( Figure 5). Liu et al. [45] studied the influence of natural fractures on the propagation geometry of hydraulic fractures through a tri-axial fracturing system and found that the horizontal differential stress and the angle between the maximum horizontal principal in situ stress and the pre-fracture are the dominating factors for the initiation and propagation of hydraulic fractures. He et al. [46] studied the hydraulic fracturing through the displacement discontinuities method and found that the pre-existing fractures were favorable to the formation of the hydraulic fracture network. Zhou et al. [47] studied the near-borehole fracture propagation in a laminated reservoir rock using PFC2D and found that the hydraulic fracture propagation in a laminated reservoir is controlled by both in situ stress and the bedding surface. Apart from the above common features, there are some intriguing differences between the models with different structural plane orientations. For models with low angles (0 • -15 • ), the structural planes control the initiation and propagation of the hydraulic fracture networks. However, as the angle increases, more small-scale transverse fractures occur in conjunction with the main fracture networks. When the orientations of the structural planes are approximately parallel to the direction of σ H (0 • to 15 • , Figure 4a,b), the initiation and propagation of hydraulic fractures tend to develop along the structural planes and form straight fractures. When the angle increases slightly (30 • to 45 • , Figure 4c,d), hydraulic fractures still occur along the structural planes, but many secondary fractures also occur. These secondary fractures are aligned in the direction of σ H and are interconnected with the main hydraulic fractures. The increase in the pore water pressure causes a high tensile stress area at the fracture tips along the structural planes. As the tensile strength is exceeded, the fractures can grow in a direction parallel to σ H .
For models with steeply inclined structural planes (60 • -75 • , Figure 4e,f), some obvious changes occur in the morphology of the fracture networks. The maximum horizontal stress plays an increasingly important role in the initiation and propagation of the hydraulic fractures. As the angle increases, more hydraulic fractures occur through the rock between the structural planes in conjunction with the fracturing along the structural planes. Furthermore, the roughness of hydraulic fractures increased significantly with the increased angle of the structural planes. In this case, both the structural planes and the maximum horizontal stress determine the evolution of the simulated hydraulic fractures.

Evolution of Hydraulic Fractures in Laminated Rocks
Hydraulic injection and fracturing is a dynamic failure process that involves a complicated hydro-mechanical coupled interaction [1]. When high-pressure fracture fluids are injected into a target formation, a sharp increase in pore pressure and adjustment of the stress field will be induced in the affected area [22]. In this case, pre-existing open fractures and new fractures are the two main ways to develop a fracture network for the desorption and migration of shale gas. The dramatic changes in pore water pressure resulting from the fluid injection can generate significant changes in the effective stresses [20]. In addition, the total stress changes induced by the high-pressure fluid are influenced by the geometry of the fractures in the reservoir and the variations in the physical-mechanical properties of the fractured reservoir and surrounding strata [44]. The critical pressure necessary for the initiation of a new hydraulic fracture or fracturing along a natural weakness plane depends on the stress state in the rock and the structural planes. In shale, tensile failure plays an important role because the tensile strength of the structural planes can be very low. In laboratory experiments, tensile failure of the rock matrix usually shows up as the initiation of a fracture and the accumulation and coalescence of many small tensile failures, leading to macroscopic shear failure of the rock mass. This is why combined tension-shear failure was considered as a widespread phenomenon during hydraulic fracturing. When the elevated pore pressure exceeds the sum of the tensile strength of the rock and the minimum principal stress, tensile failure can occur and propagate. Figure 6 shows the minimum principal stress field caused by the initiation and propagation of hydraulic fractures in the model with different structural plane orientations. As expected, when small fractures occur, they are associated with tensile stresses in the adjacent elements near the fracture tip [23]. The tensile stress and the heterogeneity in the mechanical properties initiate the fractures and drive the propagation of these fractures, resulting in the coalescence of small fractures.
It is generally known that acoustic emission occurs due to strain energy release caused by rock failure and slip along rock mass discontinuities. Acoustic emission events can be induced by fracture initiation and propagation and by failure along natural bedding planes [20]. The model can track the location and type of acoustic emission events. The element failures tend to occur in the laminated layers around the injection hole due to their lower strength and Young's modulus. Detailed information on how acoustic emission events are modeled can be found in reference [38]. Figure 7 shows a case where acoustic emissions occur during the simulation of hydraulic fracturing. Although a fracture network dominated by shear failure is expected to create the most efficient seepage pathways for shale gas recovery [43,48,49], the simulation results show that nearly all the failed elements experienced tensile failure and the hydraulic fracturing of rock is essentially a tensile failure phenomenon at the micro-level. In practice, tensile fractures close easier than shear fractures after fluid injection. This is the primary reason that the fracture fluids contain propping agents, such as sand or aluminum oxide pellets, which are suspended in the fracture fluids and used to hold the fractures open after the pumping stops.  In the early stages of fracturing, nearly all the small fractures are initiated at the weak structural planes. The locations of weak structural planes and reservoir heterogeneity determined the initiation points for the small fractures. The subsequent growth and propagation of fractures are affected by the angular relation between the orientations of the structural planes and σ H . The hydraulic fractures tend to grow parallel to the orientation of σ H . When the structural planes are inclined by a small angle to σ H , the main hydraulic fractures propagate along the weak structural planes. In this case, the weak structural planes dominate the propagation direction. In contrast, when the structural planes are inclined by a large angle to σ H , both the maximum horizontal stress and the weak structural planes dominate the propagation of hydraulic fractures. As the fractures propagate, more fractures cut across the weak structural planes and propagate parallel to the direction of σ H .

Fracturing Response to Geological and Operational Variables
In addition to a qualitative evaluation of the simulation results, the fracturing response and the efficiency of hydraulic fracturing under different geological conditions were quantitatively investigated. Key parameters, such as the stress ratio, the structural plane spacing, the strength of the structural plane, the number of sets of structural planes, and the injection rate were studied to reveal their effects on hydraulic fracturing. Here, two parameters were defined to quantitatively analyze the effects of the different factors on hydraulic fracturing. These indices are (1) the hydraulically fractured area, defined as the area of the pre-existing structural planes that experience failure, and (2) the length of hydraulic fractures, defined as the total length of the created fractures.
The hydraulically fractured area corresponds to an area with a significant pore pressure increase due to injection. The length of hydraulic fractures is the total length of the failed structural planes and newly created fractures. These parameters were quantified through computer image recognition and statistical techniques to ensure their reliability and validity.

Effect of In-Situ Stress Ratio
The in-situ stresses in a reservoir significantly affect the performance of hydraulic fracturing [26,50]. The stress ratio, which is defined as the ratio of σ H /σ h , was studied to evaluate its effect on the response to fluid injection during hydraulic fracturing. Figure 8 shows the simulated fracture networks with stress ratios equal to 1.5, 1.25, and 1.0. The fracture network zone decreases in size as the stress ratio approaches unity. In addition, the simulation results indicate that the orientation of the long axis of the fractured zone also changes with different stress ratios. Specifically, when the stress ratio equals 1.0, the direction of the fracture zone is about 70 • to σ H ; when the stress ratio equals 1.25, the direction of the fracture zone is about 45 • to the maximum horizontal stress; and when the stress ratio equals 1.5, the direction of the fracture zone is about 20 • to σ H . In addition, the simulation results show that a large stress ratio is conducive to the development of relatively more branching fractures in the direction of σ H , which can enhance the connectivity of the hydraulic fractures and add to the complexity of the fracture network morphology. Figure 9 shows the impact of the stress ratio on the hydraulically fractured area and the length of fractures. The results indicate that the fracture pattern is significantly affected by the stress ratio. For the same injection rate and duration, a higher stress ratio achieves a larger number of simulated fractures and a larger fractured area.

Effect of Structural Plane Mechanical Properties
It is generally known that the properties of joints or bedding planes have a significant impact on the deformation and failure behavior of a rock mass. In shale fracturing, the mechanical properties of the structural planes represent an important issue, which cannot be ignored in the design of hydraulic fracturing. The structural planes in the numerical models are elements with lower mechanical properties than the rock matrix. Different elastic moduli of structural planes were assigned to the models to evaluate the effect of the stiffness and strength of the structural planes on the hydraulic fracturing. Figure 10 shows the quantitative results of the simulated fracture network when the structural planes have different elasticity moduli. The fractured area and length of fractures changed slightly with an increase of the elastic modulus of structural planes, which were 30%, 50%, and 70% of the rock matrix, respectively. There was no significant change in the morphology of the fractured zone with an increase of elasticity moduli. The boundary stress and injection rate were the same as for the base model. Figure 11 shows the fractured area and length of fractures when the structural planes have different uniaxial compressive strengths. The default tensile-compression strength ratio of material in the FSD is 10. Therefore, changes in the uniaxial compressive strength of the structural planes are essentially paired with changes in the tensile strength in this study. The results suggest that the fractured area increases with a reduction in the uniaxial compressive strength of the structural planes, but the length of fractures in the fracture network decreases. Figure 12 shows the qualitative results of the simulated fracture network when the structural planes have different compressive strengths. The strength of the structural planes has an important influence on hydraulic fracturing. At a strength of 100 MPa, the hydraulic fractures initiated and propagated rapidly along the weak structural planes, and only a few small branching fractures were induced between the main hydraulic fractures, resulting in poor connectivity of the fracture network. At a strength of 150 MPa, although the zone of failed structural planes is smaller than the former case, there are more branching fractures between the structural planes. In this case, the structural planes still dominated the evolution of hydraulic fractures. At a strength of 200 MPa, which is very close to the strength of the rock matrix, the morphology of the fracture network is distinctly different. Here, the impact of the structural planes on hydraulic fracturing is not apparent, and the morphology of the fractured zone is similar to the base model.

Effect of Structural Plane Spacing
As stated before, natural structural planes typically form a key part of a hydraulic fracture network. Therefore, the structural plane spacing and the number of sets of planes can inevitably affect the morphology of the fracture network and the effectiveness of the hydraulic fracturing. A group of structural planes with different spacings was studied. For example, with an orientation of the structural planes equal to 15 • , Figure 13 shows the distribution of the simulated fracture network when the structural planes have different spacing. The structural plane spacing influences the morphology and connectivity of the fracture network. The failed structural planes dominate the development of the hydraulic fractures and control the pattern of the fracture network. Simulation results with a smaller spacing generate more complex fracture networks and better fracture connectivity. Figure 14 shows the fractured area and the total length of fractures for structural plane spacings of 1, 2, and 3 times the radius of the injection hole. The fractured area did not vary significantly for different spacings, but a reduction in spacing resulted in a decrease in the total length of fractures that were created. There may be multiple sets of structural planes in a reservoir. To study this more complex geological condition, the effect of the structural plane spacing on the hydraulic fracturing for two sets of weakness planes was also simulated. Figure 15 shows the qualitative characteristics of the simulated fracture networks for two sets of orthogonal structural planes with different spacings. In this case, the structural planes still dominated the morphology and connectivity of the fracture network. Although the model with wider spacing for the structural planes achieved a larger fractured area, it had fewer hydraulic fractures and poor connectivity compared to the model with closer spacing.

Effect of Injection Rates
The injection rate is one of the most important operational parameters when conducting hydraulic fracturing in shales [3,5]. Although it was believed that large injection rates can enhance the propagation length of fractures [8], the quantitative relationship between the injection rate and the reservoir fracturing is still largely unknown [4]. In this study, different injection rates were tested to study their effects on hydraulic fracturing. Figure 16 shows the qualitative results of the simulated fracture network with injection rates equal to 0.12, 0.15, and 0.18 m 3 ·s −1 ·m −1 , with the structural planes parallel to σ H . The injection rate affects the growth and development of the fracture network during the hydraulic fracturing. High injection rates create a larger fractured area and a longer total fracture length. The fractured area and number of fractures increased with the injection rate. For a high injection rate, the injection fluids not only fracture the structural planes, but also initiate small fractures in the rock matrix. In contrast, for a lower injection rate, the time required to reach the breakdown pressure of the shale may be longer, thus resulting in less fractured structural planes and high rates of fluid leak-off into the formation.  Figure 17 shows the quantitative results of the fractured area and length for different constant injection rates. The results show that higher injection rates result in better hydraulic fracturing. At a low injection rate, the fracturing fluid tends to only flow along the structural planes in the direction of σ H . However, at a high injection rate, multiple random branches of hydraulic fractures form, creating a more complex and connected fracture network.

Conclusions
The evolution of hydraulic fracturing in bedding or laminated rocks was investigated through a micromechanics-based 2D numerical modeling approach. The numerical models incorporated heterogeneity in the mechanical properties of the rock and all models other than the base model included planes of weaker elements to simulate the presence of structural features, such as bedding planes in the rock. The mechanical-hydraulic coupled process between the rock and the injected fluid was also incorporated in the models. A series of simulations were conducted to explore the effect of the rock's mechanical properties and fluid injection rates on the initiation and propagation of fractures and the resulting morphology of the fracture network created.
The process of hydraulic fracturing in laminated rock can be divided into five stages: stress accumulation, small fracture formation, hydraulic fracture initiation, hydraulic fracture propagation, and decay in fracture propagation. Heterogeneity in rock properties influences the fracture initiation and propagation, and the location and orientation of fracture initiation are unpredictable.
The angle between the orientation of the structural planes and the maximum horizontal stress has a significant effect on the propagation and evolution of the fracture network. When the structural planes are nearly parallel with the direction of σ H , the weak structural planes dominate the propagation direction and the morphology of the fracture network. When the orientation of the structural planes is significantly different from the direction of σ H , the extent and morphology of the fracture network are controlled by both the in-situ stresses and the orientation of the structural planes. In this case, many small transverse secondary fractures can occur in conjunction with fracturing along the structural planes. The acoustic emission events induced by the hydraulic fracturing demonstrated that most of the failed elements are caused by tensile failure at the microscale.
The in-situ stress ratio has a great impact on the morphology and extent of the fracture network. For a fixed injection rate and duration, a high-stress ratio creates a larger area of fractured rock with a better connectivity. The stiffness of the structural planes has minimal influence on the extent and morphology of the fracture network, but the mechanical strength of pre-existing weakness planes, the plane spacing, and the injection rate have a significant effect on the hydraulic fracturing. Small structural plane spacing, a high strength contrast, a high-stress ratio, and a high injection rate result in a complex interconnected network of hydraulic fractures.