A Numerical Study on the Effect of Anisotropy on Hydraulic Fractures

In this paper, we present a two-dimensional numerical model for modelling of hydraulic fracturing in anisotropic media. The numerical model is based on extended finite element method. The saturated porous medium is modelled using Biot’s theory of poroelasticity. An enhanced local pressure model is used for modelling the pressure within the fracture, taking into account the external fluid injection and the leak-off. Directional dependence of all the rock properties, both elastic and flow related, is taken into account. A combination of the Tsai–Hill failure criterion and Camacho–Ortiz propagation criterion is proposed to determine the fracture propagation. We study the impact on fracture propagation (in both magnitude and direction) due to anisotropies induced by various parameters, namely ultimate tensile strength, Young’s modulus, permeability and overburden pressure. The influence of several combinations of all these anisotropies along with different grain orientations and initial fracture directions on the fracture propagation direction is studied. Different regimes are identified where the fracture propagation direction is controlled by the degree of material anisotropy instead of the stress anisotropy.


Introduction
Hydraulic fracturing is a process of inducing fractures in rock structures by injecting fluid at high pressures. Interest in hydraulic fracturing has been increasing in recent years due to its applications in the oil and gas industry for enhanced oil recovery in conventional reservoirs, for heat recovery from geothermal reservoirs and for the profitable extraction of oil and gas from unconventional reservoirs. A better understanding of the fracture growth phenomena will enhance productivity and also reduce the environmental footprint as less fractures can be created in a much more efficient way.
Several models have been developed for this purpose, the earliest of which was by Khristianovic and Zheltov (1955) 1 3 and later extended by Geertsma and Klerk (1969) to form the KGD model, which proposed an analytical solution to the problem by assuming plane strain conditions. Although this model has been accepted as a standard case for hydraulic fracturing due to its simplicity, it suffers from the assumption that there is no leak-off from fracture into formation and also the formation is assumed to be solid. Another analytical solution with a different geometrical assumption (fracture length ≫ fracture height) was given by Perkins and Kern (1961) and extended to include fluid loss by Nordgren (1972). Analytical asymptotic solutions were formulated (Adachi and Detournay 2002;Garagash and Detournay 2005;Detournay 2016), based on a parametric space to identify the significant parameters in a given regime. Later, these asymptotic solutions have been modified to take into account the leak-off, fracture toughness and fluid viscosity (Adachi and Detournay 2008;Kovalyshen 2010;Dontsov 2017).
The first numerical models for modelling hydraulic fracturing were developed by Boone and Ingraffea (1990) using finite elements for modelling the formation and using finite volume for modelling flow with cohesive zones along element edges describing the fractures. In recent times there have been several models developed based on different numerical techniques. A linear elastic fracture mechanics (LEFM)-based finite element model (FEM) was proposed by Hossain and Rahman (2008). To avoid the singularity problems at the crack tip in LEFM, FEM models with zerothickness elements (describing the fractures using cohesive zones) were developed (Carrier and Granet 2012;Chen 2012). The FEM approach requires re-meshing to capture the fracture propagation accurately, whereas extended finite element (XFEM)-based models allows for fracture propagation in arbitrary directions without the need for re-meshing (Mohammadnejad and Khoei 2013;Remij et al. 2015;Meschke and Leonhart 2015). A novel approach in which the asymptotic behaviour near the fracture tip was resolved with extended finite element method Peirce 2013a, b, 2015). An alternate approach based on phase field modelling which combines FEM with continuum damage mechanics has been developed (Mikelic et al. 2015;Miehe and Mauthe 2016) which provides a convenient way for modelling complex fracture interactions. But all the above hydraulic fracture models assume the rock formation to be isotropic in nature.
Most rocks (especially shales, which are the most common rock type to be hydraulically fractured) are highly anisotropic in nature (Jaeger et al. 2009;Barton 2007). Kaarsberg (1959) and Sayers (1994) observed that shales have a bedding plane along which the grains are oriented causing the properties along the grain direction to be vastly different from the properties perpendicular to the grain direction. This causes a special type of anisotropy, called transverse isotropy, where the material properties in any direction in the plane can be obtained by using the material properties along any two mutually perpendicular set of directions in that plane. Although there are several studies (Abousleiman et al. 2008;Zhubayev et al. 2015) experimentally obtaining the anisotropic parameters, there are few papers by Cheng (1997) analytically deriving the anisotropic poroelastic coefficients. Abousleiman et al. (1996) modelled the deformation and pressure in a transversely isotropic porous medium without any fracture. Porous material with a stress-driven fracture in an orthotropic medium was modelled by Remij et al. (2015). More recently, the influence of rock anisotropy on tensile fractures was studied experimentally by Mighani et al. (2016).
In this paper, we enhance the aforementioned XFEM model by Remij et al. (2015) to include the effects of anisotropy on hydraulic fracturing. Using the model we analyse the effect of anisotropic rock properties (Young's modulus, ultimate tensile strength and permeability) on the fluid-driven fracture propagation and also the impact when combined with loading.

Mathematical Formulation
In order to solve a poroelastic problem, we need to solve for the solid deformation ( ) and fluid pressure (p) at all points within the porous media. Biot's theory of poroelasticity (Biot 1941) is used to describe the porous media, and fluid flow is described using Darcy's law. The porous medium is assumed to be saturated.
In addition for the hydraulic fracture problem, we make use of an enhanced local pressure (ELP) as proposed by Remij et al. (2015) to model the pressure inside the fracture at all points along the fracture length. This additional degree of freedom ( p d ) enables us to model the complicated phenomenon happening within the fracture, namely the injection of external fluid, moving boundaries of fracture surface and the leak-off. Leak-off from fracture to formation is described using the 1-D Terzaghi's consolidation equation.
For solving the unknowns, a set of governing equations along with auxiliary equations are used. The governing equations used to describe the poroelastic problem are of two types: solid deformation-based momentum balance and fluid flow-based mass balance. We consider an additional equation to ensure the mass balance inside the fracture. The auxiliary equations are used for relating these governing equations with the unknowns and also for coupling them. A schematic flow chart of the mathematical formulation is represented in Fig. 1.

3
A schematic of a body with a discontinuity d , which splits the body into two domains + and − , along with prescribed boundary conditions is represented in Fig. 2.

Discretisation
We need a numerical method in order to solve the set of coupled differential equations and also incorporate the discontinuous jump in various parameters due to the fracture. Hence, we make use of XFEM which models the discontinuous jump due to fractures by using additional degrees of freedom. XFEM allows for the fracture to propagate through the elements, thus ensuring accuracy in capturing the fracture even with a coarse mesh making it computationally very efficient. Similar to traditional FEM, the unknowns are obtained at certain control points (nodes) by solving the differential equation and the unknowns in the intermediate regions are obtained by interpolation using shape functions. The discretised form of the unknown displacement and pressure in a porous medium intersected by the fracture is expressed as: where H d is the Heaviside step function given as: u and p represent the regular nodal degrees of freedom for displacement and pressure, respectively, while ũ and p represent the enhanced nodal degrees of freedom representing the discontinuous jump in displacement and pressure across the fracture. p d represents the pressure inside the fracture at points where the fracture intersects the element edges (Fig. 3). N, L are two-dimensional interpolation or shape functions for the displacement and pressure fields, whereas V is a one-dimensional shape function for interpolation of the pressure along the fracture length. (1) Schematic of a domain with discontinuity d 1 3

Solution
The governing equations are combined with the auxiliary equations as shown in Fig. 1. Weak form of this set of equations is obtained by integrating them along with a test function. By substituting the discretised unknowns given by Eqs.
(1), (2) and (3) into the weak form, we convert the set of differential equations into a set of algebraic equations. In order to solve this set of equations simultaneously, we make use of the Newton-Raphson iterative solver in combination with Euler's forward scheme for obtaining the time derivative, and Euler's implicit scheme for time-independent parameters. A detailed description of the solution procedure is given in Remij et al. (2015).
The unknown X = ûũppp d T degrees of freedom are solved at each grid point (nodes) for every time step.

Propagation
A propagation criterion is needed in order to determine the propagation initiation and also the magnitude and direction of propagation. Hence, we make use of the Camacho-Ortiz criterion (Camacho and Ortiz 1996) which can be used for mixed mode fractures as well. This criterion states the condition for fracture propagation as: where t eq is the equivalent traction ratio in a specific direction, ult and S ult represent the ultimate tensile and shear strength of the porous media, respectively, t n and t s are the normal and shear tractions along that orientation which are obtained from average stresses as: where is the coefficient of friction, and n and s are the unit normal and tangent vector to the direction. av is the average stress, used to obtain a better approximation of the stress state in the vicinity of the fracture tip. The average stresses are obtained by assigning weight functions to the Gaussian integration points within a certain distance (generally three times the characteristic element length) from the fracture tip, as derived by Jirasek (1998). As an averaged stress is used, this may lead to a slight delay in the onset of fracture propagation. The direction of propagation is taken to be the direction in which the equivalent traction is maximum. The fracture is assumed to propagate through the entire element length in a single time step in a straight line. Further details on the implementation of the solution are described by Remmers (2006) and Remmers et al. (2003).

Anisotropic Parameters
In this section, we highlight the parameters which need to be modified to incorporate the effect of anisotropy.

Constitutive Relation
The effective elastic stress ( e ) in the solid grain is related to the elastic strain by means of the generalised Hooke's law.
where C is the constitutive relationship matrix and is the elastic strain in the solid grains. The coefficients of the constitutive matrix depend on the material type and geometrical assumptions. The constitutive relation matrix can be obtained as the inverse of the compliance matrix, C = S −1 . We assume a plane strain case with a transverse isotropic material for which the compliance matrix in the grain direction is given as: where E ∥ and E ⟂ are the Young's moduli parallel and perpendicular to the grain direction and in is the in-plane Poisson's ratio, representing compressive strain perpendicular to the grain direction due to a tensile stress parallel to the grain direction and out represents the out-of-plane Poisson's ratio.
The constitutive matrix at any arbitrary direction is obtained from the following expression: where T is the transformation matrix given as a function of the angle ( ) between the global direction and the grain orientation direction.

Ultimate Tensile Strength
The ultimate tensile strength ( ult ) is an important parameter which determines the fracture propagation. The ultimate tensile strength is maximum in the grain orientation direction and minimum perpendicular to it. Some previous studies (Remij et al. 2015;Yu et al. 2002;Lee and Pietruszczak 2015) have assumed cosine functions to interpolate the ultimate tensile strength at arbitrary directions from the values along the grain direction and the perpendicular direction.
Here we make use of the Tsai-Hill failure criterion (Jones 1998) to model the directional dependence of the ultimate tensile strength. The Tsai-Hill failure criterion for a 2-D transverse isotropic material is given as: sin 2 2 sin cos sin 2 cos 2 −2 sin cos − sin cos sin cos cos 2 − sin 2 ⎤ ⎥ ⎥ ⎦ where || , ult || and ⟂ , ult ⟂ are the stresses, ultimate tensile strengths parallel and perpendicular to the grain direction, respectively, and s and S ult indicate the shear stress and ultimate shear strength in the plane.
Assuming a fracture propagating at an angle with respect to the x-axis, the ultimate tensile strength perpendicular to the fracture propagation direction frac is given as the normal stress perpendicular to fracture propagation direction which can satisfy the Tsai-Hill failure criterion. To obtain frac , this stress state ( frac = 0 frac 0 T ) is rotated to the grain orientation direction by using the stress transformation relations and then substituted into Equation (12).
where = − , is the grain orientation angle with respect to the x-axis, and is the fracture propagation angle with respect to the x-axis (Fig. 4). Hence, for a randomly oriented fracture, we obtain the ultimate tensile strength in the direction perpendicular to the fracture propagation direction ( + 90 • ) as: In order to validate the proposed Tsai-Hill-based theory, we make a comparison with experimental results. Mighani et al. (2016) conducted tensile fracture experiments on Lyon's sandstone and pyrophyllite rocks and observed the ultimate tensile strength variation with which is the angle between the grain orientation direction and the fracture propagation direction. We make use of the experimental values for maximum and minimum ultimate tensile strengths and interpolate for various angles. As one can observe in Fig. 5, Tsai-Hill failure criterion provides a much better fit for the experimental values when compared to the previously assumed cosine functions (Remij et al. 2015;Lee and Pietruszczak 2015). Cheng (1997) derived the analytical expressions for the transverse isotropic poroelastic coefficients based on the constitutive relationship matrix. In a fully saturated porous medium, the external stresses on the porous media are partly taken by the fluid pressure in the pores and partly by deformation of the solid grains. This can be represented mathematically as:

Poroelastic Coefficients
where is the total stress, e is Terzaghi's effective stress, and is Biot's coefficient matrix.
The anisotropic Biot's coefficient matrix reduced to the 2-D form is given by the expression: The other poroelastic constants such as bulk modulus and compressibility modulus are given as: where n f is the porosity of the porous media, and K s and K f are the bulk modulus of the solid and fluid, respectively.

Permeability
Permeability enters into the formulation through Darcy's law which describes the fluid flow in the porous medium as: where is the intrinsic permeability tensor, is the dynamic viscosity of the pore fluid, q is the flux, and p refers to the pressure in the porous media. Permeability tensor for an anisotropic rock in the principal grain directions is given as where permeability in the grain direction ∥ is larger than the permeability perpendicular to the grain direction ⟂ . For permeability in any arbitrary orientation we make use of the transformation matrix:

Validation
Since there are no studies which exactly deal with hydraulic fracturing in anisotropic media, we divide the validation into two parts: (1) Mandel's problem which compares the numerical results with an analytical solution for a transverse isotropic porous medium without fractures and (2) the standard KGD problem which compares the numerical results with the ELP model (Remij et al. 2015) for a hydraulic fracture problem in an isotropic medium. Abousleiman et al. (1996) provided an analytical solution for Mandel's problem in a transversely isotropic porous medium. Mandel's problem (Fig. 6) consists of an infinitely long rectangular block with the left and right ends free from stresses and the fluid is free to flow, whereas an external force is applied on the top and bottom boundaries. The external force is taken as 10.5 MPa. The Young's moduli along horizontal and vertical directions are 20 and 10 GPa,

Mandel's Problem
cos . respectively. In-plane Poisson's ratio is assumed to be 0.30, whereas the out-of-plane one is taken as 0.20. Similarly, permeability values are taken to be 10 −19 and 10 −17 m 2 along horizontal and vertical directions. The bulk moduli of the solid grains and the pore fluid are assumed to be 36 and 3 GPa. A time step of 500 s is used.
We compare the pore pressure solution from the numerical model with the analytical solution at different time periods in Fig. 7. The numerical pore pressure decay from the centre of the specimen to the free edges is found to be consistent with the analytical solution with relative errors (< 5%). The displacement in the x-direction along the centre line of the specimen is plotted and compared with the analytical solution.

KGD
In this validation case we consider a KGD problem (Fig. 8) which is a standard test case for hydraulic fracture problems. When the rock is assumed to be isotropic, there exists a theoretical solution given by Geertsma and Klerk (1969).
The Young's modulus and Poisson's ratio are taken as 20 MPa and 0.2. The ultimate tensile strength is assumed to be 2 MPa, while the toughness is 120 N/m . The permeability and viscosity are given as 10 −19 m 2 and 0.1 Pa s . The initial fracture at the boundary is injected at the rate of 25 mm 2 ∕s for a time period of 100 s with a time step of 0.1 s. The KGD problem considered here lies in the viscosity-storage propagation regime.
The fracture propagates on a non-predefined path. In Fig. 9, we compare the fracture profiles at various time steps from the current numerical model and the ELP model (Remij et al. 2015). As observed the current model accurately reduces to the ELP solution for isotropic values of the parameters. The fracture mouth opening pressure variation with time is also plotted and compared. A mesh of 50 × 50 mm 2 is used for the purpose.

Vertical Hydraulic Fracture Problem
The test case (Fig. 10) that is used here is very similar to the previous KGD problem, but an initial crack is placed at the bottom of the model, to model the vertical fracture growth representative of the hydraulic fractures. The model is simulated for 10.5 s with a time step of 0.1 s. Also the model is subjected to external stresses of 40 MPa, due to the overburden pressures or in situ stresses existing at depths of around 1.5 km. The fluid is free to flow from the top and the right boundaries, whereas no-flow conditions exist at the left and bottom boundaries. The isotropic values of the parameters used in the model are specified in Table 1. All angles are with respect to the horizontal axis. For cases of anisotropy, these isotropic values are perturbed depending on the degree of anisotropy. The degree of anisotropy is defined as: By using the parameters in Table 1 to obtain the nondimensional parameter (  k ) described in Bunger et al. (2005), we observe that the hydraulic fracture problem described here lies in the viscosity-dominated regime and closer to the storage edge. In the following Sects. 5.3 and 5.6.1 we try to understand the influence of anisotropy in each individual parameter by keeping all other parameters isotropic. We also look at the possible combination of anisotropy in these parameters in Sects. 5.4 and 5.6.2.

Parametric Anisotropy
In this subsection, we vary one parameter at a time to find out the fracture propagation variation with anisotropy in each individual parameter. In all the considered test cases we assume that the grains are oriented along the horizontal direction ( 0 • ) and the initial fracture is oriented in the vertical direction ( 90 • ).

Anisotropy Due to Young's Modulus
We consider three possible scenarios for varying Young's modulus: (1) E-Parallel: anisotropy caused by increasing the From the fracture length plot in Fig. 11, we can observe that E ∥ has a much greater effect on fracture propagation than E ⟂ . This is due to the fact that the propagation of the initial vertical fracture is dependent on the stresses which are perpendicular to it. Hence, an increase in E ∥ results in higher stresses perpendicular to the initial fracture, which promotes fracture growth significantly, whereas a decrease in E ⟂ only has a smaller Poisson's effect on the stress. Since E ∥ > E ⟂ always, the fracture prefers to orient itself perpendicular to the grain orientation which is observed in all the three scenarios in the fracture orientation plot. Also we plot the pressure at the mouth of the fracture, which is a much more easily measurable quantity in the field. Since all the three scenarios favour fracture propagation in the same initial fracture direction, there is little variation ( < 5% ) in the pressure required to open the fracture.

Anisotropy Due to Ultimate Tensile Strength
Similar to the Young's modulus variation, we consider the same three scenarios for understanding ultimate tensile strength-induced anisotropy. Fracture propagation is resisted by the ultimate tensile strength perpendicular to the fracture orientation. Hence, the fracture tends to propagate along the direction perpendicular to the minimum ultimate tensile strength. Since ult ⟂ is always lower than ult ‖ , the fracture tends to propagate parallel to the grain orientation. But since the initial fracture is oriented in an unfavourable direction (perpendicular to the grain direction), the fracture continues to propagate in its initial direction until a threshold level where the effect of anisotropy becomes significant to rotate the fracture as observed from the fracture orientation plot in Fig. 12.
Also looking at the fracture length variation we observe that there is a significant increase in the fracture length when the fracture re-orients itself from its initial direction to the favourable direction. Since the fluid inside the fracture has to go through steep rotation ( ∼ 80 • ), much higher pressures are required in order to drive the fracture.

Anisotropy Due to Permeability
Anisotropy in the permeability of the rocks was considered in the formulation. As indicated in Table 1, the isotropic permeability of shales was assumed to be 10 −19 m 2 (100 nd). It was varied within two orders of magnitude, i.e. 10 −18 to 10 −20 m 2 (1000-10 nd). However, its impact on the fracture growth was found to be very negligible since shales already have very low permeability values (almost impermeable).

Degree of Material Anisotropy
In the following subsections, we focus only on the fracture propagation direction. This is because variation in fracture length or width due to various anisotropies and combinations can be overcome by varying the fluid injection time, but the fracture orientation direction cannot be modified by means of any external influence as it is solely dependent on the field conditions. Hereafter, all the anisotropies considered are by varying both the values parallel and perpendicular to the grain direction as given by Eq. (22).
In this set of cases, we study the interplay between anisotropy due to Young's modulus and the anisotropy due to ultimate tensile strength by varying them from 0 to 75% individually. As seen from the previous Sect. 5.3, both these anisotropies have contrasting effect on the fracture propagation direction. Therefore, it is important to identify the regions where one parameter has a higher degree of influence than the other. Figure 13 represents a schematic of the material anisotropy contour plots shown in Fig. 14. As can be seen there are two distinct regimes: (1) Regime A, influenced by anisotropy in ultimate tensile strength causing the fractures to propagate parallel to the grain direction, (2) Regime B, influenced by anisotropy in Young's modulus resulting in fractures finally getting oriented perpendicular to the grain direction. There is a sudden transition from one regime to another once a threshold value is crossed.
In Fig. 14, we observe that contour plots are represented for four different grain orientation ( ) directions. refers to a grain orientation with respect to the global horizontal direction (x-axis). The non-smooth variations in the threshold values for transition between regimes in the contour plots are due to the variation of anisotropies in step sizes of 5%. Looking at the different contour plots we can see that the size and shape of the different regimes vary with varying grain orientation angles.
From the = 0 • plot, we observe that beyond a threshold value of 50% DOA in ult the fracture moves from Regime B to Regime A (red regions representing < 20 • ). Looking at the plot for = 30 • , we observe that the fracture tries to align itself perpendicular to the grain direction (blue regions representing > 100 • ) even for relatively low values of Looking at all the contour plots we can observe that the area of Regime A increases with grain orientation angle. The influence of ult anisotropy is much higher as the grain orientation angle increases due to the fact that the fracture needs to be rotated by a smaller angle from its initial vertical orientation to align itself with the grain orientation direction. The converse is true for the influence of Young's modulus anisotropy.

Angle of Orientation
In this set of cases, we vary the degree of anisotropy of the material parameters (namely both Young's modulus and ultimate tensile strength simultaneously) between 0 and 75% as well as the grain orientation direction from 0 • to 90 • . The results of these cross-sets of cases are represented in the form of a contour plot as shown in Fig. 15.

Grain Orientation Direction
We observe that along the vertical axis (for a particular grain orientation direction), after a certain degree of combined 1 3 material anisotropy the crack always tends to align with the grain direction. This is due to the fact that ult anisotropy has a greater impact than the E anisotropy at larger values of DOA ( > 30%), whereas E anisotropy has a much larger impact when the DOA is lower ( < 30%). Along the horizontal axis (for a particular combined material DOA), we observe that the fracture tends to re-orient itself with relative ease as the angle ( = − ) between the initial fracture orientation and the grain orientation reduces.
In Fig. 16, the fracture propagation at a constant combined material degree of anisotropy of 55% for various grain orientation angles from 0 • to 90 • is plotted.

Initial Fracture Orientation
In all the previous and later cases, the initial fracture is assumed to be along the vertical direction ( 90 • ). In this subsection, three initial fractures at varying initial orientation angles of 90 • , 60 • and 0 • are included simultaneously. The grain orientation angle is assumed to be 60 • with a combined material DOA of 15%.
From Fig. 15 for a material DOA of 15%, we observe that the fracture aligns along the grain orientation direction when < 50 • and aligns perpendicular to the grain direction when > 50 • . For the initial fractures 90 • , 60 • and 0 • corresponding theta values are 30 • , 0 • and 60 • , respectively. Hence, the first two initial fractures ( 90 • and 60 • ) are expected to be aligned along the grain orientation direction ( 60 • ), while the 0 • initial fracture is expected to be aligned perpendicular to

Anisotropy Due to Lithostatic Stress
Similar to other parametric variations, we consider the three possibilities of stress-induced anisotropy due to (1) an increase in vertical stresses, (2) a decrease in horizontal stress, (3) a combination of both. Note that in real hydraulic fracture scenarios in the field, the vertical overburden pressure is assumed to be slightly larger than the horizontal pressure. The fracture orients itself preferably in the direction parallel to the maximum compressive stress, which is same as the initial fracture direction ( 90 • ) in all the three scenarios. The decrease in horizontal stresses results in much lower compressive stresses perpendicular to the fracture, thereby promoting fracture growth and also requiring less pressure to open the fracture as seen from Fig. 18. Similarly, the increase in vertical stresses only causes a mild increase in the tensile stresses acting perpendicular to the fracture due to Poisson's effect.

Combined Anisotropy
In this final scenario we look at combining all the abovediscussed anisotropy scenarios. The material anisotropy refers to varying both Young's modulus and ultimate tensile strength simultaneously from 0 to 75% by using Eq. (22). The stress anisotropy refers to variation of the overburden pressures from 0 to 50% given by Eq. (22). The contour plots are obtained for four different grain orientation angles ( ), namely 0 • , 30 • , 60 • and 90 • . Looking at the schematic (Fig. 19) for the combined anisotropy contour plots shown in Fig. 20, we observe that apart from Regimes A and B, there is an additional Regime C in which the in situ stresses play a major role in determining the final fracture orientation. The fractures get oriented parallel to the maximum compressive stress direction in this regime. Looking at the plot for = 0 • , as already seen in Fig. 14 the fracture tends to re-orient itself parallel to grain direction after a combined material DOA of 50%. But this re-orientation is possible only when the stress anisotropy is less than 20%, beyond which the fracture transitions into Regime C. In this case, Regime B and Regime C coincide with each other as perpendicular to the grain direction and maximum compressive stress are both in the vertical direction.
The plot for = 30 • is the most complete plot with all the regimes which are influenced by the various anisotropy parameters. The dark blue regions representing fracture orientations ( ) larger than 110 are indicative of Regime B at low material DOA ( < 25% ) and low stress DOA ( < 25%). The high material DOA ( > 25% ) and low stress DOA ( < 30% ) regions represented by varying intensities of red colour (Regime A) are indicative of the ultimate tensile strength influence. For high stress DOA ( > 30% ) (Regime C), the fractures are more influenced by the stress-induced anisotropy. But the fracture orientation angles ( ) are not exactly 90 • which it is supposed to be as the vertical overburden pressures are maximum. This is because of the fact that although external stresses are maximum in one direction the local stress state has maximum values in a different direction as a result of the Young's modulus anisotropy. Therefore, the fracture ends up oriented at angles ( ∼ 100 • ) in between perpendicular to the grain orientation ( 120 • ) and the vertical direction ( 90 • ). When = 60 • , the ultimate tensile strength has great influence over most of the regions except the regions with high stress DOA ( > 30% ) and high material DOA(> 40% ) where both the stress and Young's modulus anisotropy combine. When = 90 • , both the combined material DOA and the stress anisotropy prefer the fracture to propagate along its initial vertical direction causing all the regimes to coincide.
Looking at all the four contour plots together one can observe the reduction in the influence of stress-induced anisotropy as the grain orientation direction increases, or in other words when the angle ( = − ) between the initial fracture orientation and the grain direction decreases, while the converse is true for material-based anisotropy.

Conclusion
From all the above simulations, we observe that the hydraulic fractures are greatly influenced by the anisotropy arising due to various parameters. Some of the important observations are as follows: (a) Young's modulus anisotropy promotes fracture growth perpendicular to the grain direction. (b) Ultimate tensile strength anisotropy promotes fracture growth parallel to the grain direction. (c) Stress-induced anisotropy promotes fracture growth parallel to the maximum overburden pressure. (d) At high degrees of material anisotropy, ultimate tensile strength has a greater influence than the Young's modulus, while the converse is true for low degrees of anisotropy. (e) Most important angle influencing fracture orientation is the angle between grain orientation and the initial orientation. When this angle decreases, the influence of ultimate tensile strength anisotropy increases, while the influence of Young's modulus anisotropy and stressbased anisotropy decreases.
Combination of stress anisotropy, material anisotropy and the initial fracture orientation with respect to the grain orientation is observed to determine the final fracture propagation direction.