Lattice Numerical Simulations of Hydraulic Fracture Propagation and Their Geometry Evolution in Transversely Isotropic Formations

Accurate prediction of the fracture geometry before the operation of a hydraulic fracture (HF) job is important for the treatment design. Simplified planar fracture models, which may be applicable to predict the fracture geometry in homogeneous and continuous formations, fail in case of fractured reservoirs and laminated formations such as shales. To gain a better understanding of the fracture propagation mechanism in laminated formations and their vertical geometry to be specific, a series of numerical models were run using XSite, a lattice-based simulator. The results were studied to understand the impact of the mechanical properties of caprock and injection parameters on HF propagation. The tensile and shear stimulated areas were used to determine the ability of HF to propagate vertically and horizontally. The results indicated that larger caprock Young’s modulus increases the stimulated area (SA) in both vertical and horizontal directions, whereas it reduces the fracture aperture. Also, larger vertical stress anisotropy and tensile strength of caprock and natural interfaces inhibit the horizontal fracture propagation with an inconsiderable effect in vertical propagation, which collectively reduces the total SA. It was also observed that an increased fluid injection rate suppresses vertical fracture propagation with an insignificant effect on horizontal propagation. The dimensionless parameters defined in this study were used to characterize the transition of HF propagation behavior between horizontal and vertical HFs.


INTRODUCTION
Hydraulic fracturing is a widely used stimulation technology in unconventional reservoirs. Due to the sedimentation process, the reservoirs and the over-and underlying formations are formed in a laminated form, which is known as the transverse isotropic medium. In general, sediments are horizontal layers, with the axis of symmetry perpendicular to the lamination; therefore, they are referred to as a vertical transversely isotropic medium (VTI or TIV). Figure 1A shows an example of fracture propagation from a stiff layer to another one, penetrating through a soft layer in between. Based on several published studies, it is well known that the fracture geometry in laminated formations is a function of stress anisotropy, the contrast between the mechanical properties of the reservoir formation and caprocks, and, to some extent, the operational aspects, such as the injection rate and fluid properties Li et al., 2017;Zeng et al., 2018;Aimene et al., 2019;Dou et al., 2019). Understanding hydraulic fracture (HF) propagation in laminated reservoirs is important to accurately predict the fracture geometry for the treatment design. Figure 1B presents the four observed interaction mechanisms (modes) when an HF intersects a natural interface (e.g., natural fracture or lamination). These are crossing, opening (HF only propagates along the interface), offsetting (HF reinitiate on the other side of the interface after a short propagation along the interface), and arresting (HF cannot cross or open the interface) mechanisms . Bakhshi et al. (2019) performed hydraulic fracture simulations based on a lattice-based method to investigate the interaction modes between the hydraulic fracture and natural interfaces considering the variable shear strength of interfaces. All simulation results are consistent with the laboratory experiment results (Sarmadivaleh and Rasouli, 2015) and proved the applicability and accuracy of the lattice-based method, which will be used in this article.
Based on the laboratory experimental and field data, some researchers concluded that the contrast in the mechanical properties of adjacent layers primarily determines the fracture height growth (Simonson et al., 1978). They outlined that HF would be contained if the stiffness of the pay zone was less than that of the adjacent caprock layers; otherwise, fracture penetration would occur. However, some scholars claimed that the fracture containment was ascribed not only to stiffness between layers but also to interface properties and interlayer stress differences (Warpinski et al., 1982;Teufel and Clark, 1984;Smith et al., 2001;Gu and Siebrits, 2005;Daneshy, 2007;Zhou et al., 2017). In addition, the interfacial shear strength and the angle of approach between the HF and natural interface may play an important role in HF containment. Essentially, the effect of the angle of the approach is also realized by affecting the shear strength of the interface. Some laboratory experiments performed to study the influence of the interfacial shear strength effect on HF interactions with natural fractures showed that a strong interface reduces the possibility of interfacial slippage but benefits HF penetration (Sarmadivaleh, 2012;Sarmadivaleh and Rasouli, 2015).
In recent decades, with the rapid development of computer science, numerical techniques have been developed to simulate HF based on different methods such as the boundary element method (BEM), finite element method (FEM), extended finite element method (XFEM), finite difference method (FDM), displacement discontinuity method (DDM), discrete element method (DEM), and some hybrid methods (Wu and Olson, 2015;Khoei et al., 2018;Vahab et al., 2019;Xie et al., 2020). Some fluid-solid models, based on the DEM, were conducted to demonstrate that the vertical heterogeneity can influence HF extension (Chen et al., 2015;Haddad et al., 2017;Wang et al., 2017;Zhang et al., 2017;Liu et al., 2019;Liu et al.,2020). The interaction modes between the HF and discontinuities were modeled to predict the possible fracture extension path (Zhang and Jeffrey, 2012;Weng et al., 2018).
In this article, XSite, a DEM-based software with high computational efficiency, was used to simulate HF propagation and its geometry in laminated formations (Damjanac et al., 2013). The tension and shear damaged zones, which are known as stimulated reservoir areas (SRAs), and the potential slippage of HF and the natural interface were determined. The impact of mechanical properties of the reservoir formation and caprocks, the stress anisotropy, natural interface properties, and injecting fluid properties was investigated. The results are presented in the following sections.

XSITE FORMULATION
The lattice simulation used in this study consists of a series of quasi-randomly distributed, non-linear spring-connected nodes as a representation of the nodal and rock matrices. Fluid flow throughout the node network, including pre-existing joints and FIGURE 1 | (A) Hydraulic fracture propagation in the alternation of stiff and soft layers (Afsar, 2014). (B) Interaction modes between the hydraulic fracture and a natural interface (Sarmadivaleh, 2012).
Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 787736 2 newly developed fractures as well as the rock matrix, and the resulting pressures are used to calculate the effective stresses (i.e., conditions that cause deformation and affect damage). To better simulate non-linear behaviors, such as fracture, sliding, and fracture opening/closing, all codes use display solutions. The discrete fracture network (DFN) is overlaid on a lattice spring network, assigning elasticity and strength parameters of the crack to the spring. If two nodes of a spring are placed on opposite sides of a joint plane, the spring will obey the smooth joint model (SJM) approach. During the analysis, only the direction of the discontinuity plane is considered, rather than the direction of the individual springs along the joint plane. The angle of the plane becomes the dominant direction when the joint plane cuts the spring (instead of the direction of the spring) .
The following central difference formulas of linear momentum equilibrium and displacement-velocity relation are used for each node to simulate the translational motion of each node : where, _ u (t) i and u (t) i are the velocity and position of component i (i 1, 3) at time t, respectively, and F i is the sum of all force components acting on mass within the time step of t. Likewise, the angular velocities, ω i , of the component can be calculated as follows: where, M i is the sum of all moment components acting on the node of the moment of inertia I. Flow in a crack, either pre-existing (specified as model input) or newly created (by breaking lattice springs), is solved in a network of fluid nodes connected by pipes (onedimensional flow cells). Fluid pressure is present in fluid nodes that act as microcracks located on broken springs or springs crossed by pre-existing joints. The crossed cracks are connected by pipes that allow the fluid to flow between the fluid nodes. Specifically, fluid nodes are connected to all fluid nodes within a distance equal to the lattice resolution multiplied by the fluid resolution, which is a user-specified dimensionless parameter between 0.6 and 1.2 (default value is 0.8). The flow velocity is calculated for each flow pipe. The geometry of the flow model is a function of the crack geometry in the solid model and is automatically generated and updated according to the evolution of the solid model (i.e., new cracks are generated). Initially, at the beginning of the simulation, the flow network is  generated in a specified pre-existing joint (i.e., DFN). Due to the forces in the joints, the joints' springs break, resulting in microcracks, so the code automatically creates new fluid nodes and connects them using flow pipes based on the spatial relationship with the existing flow network .

MODEL SETUP
In this study, a two-caprock layer model with laminations in a horizontal direction was established for simulation purposes. The model, as shown in Figure 2A, has a length (along the X-axis), width (along the Z-axis), and height (along the Y-axis) of 10, 8, and 8 m, respectively. The reservoir layer of 2 m height is located in the middle of the model. The caprocks are homogenous materials with 1 m thickness placed above and below the reservoir formation. The six zero-thickness interfaces 1 m apart from each other, as labeled from 1 to 6 in Figure 2B, were used to characterize the cemented natural interface between layers. The principal stresses were considered as σ h 5 MPa, σ v 10 MPa, and σ H 8 MPa in X, Y, and Z directions, respectively. It is important to note that in order to speed up the simulations, the stress differences were applied to the model, rather than the total stresses. This will have no impact on the fracture geometry but  Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 787736 only fracture pressures. The fracture is initiated from a cluster in the middle of a horizontal wellbore, which is placed along the X-axis (or σ h ) direction. The spherical cluster is the point of fracturing fluid injection, and a small starter crack (notch) is placed perpendicular to σ h in order to facilitate fracture initiation. The mechanical properties of this model, which are typical values of a tight unconventional formation, are listed in Table 1. The injection fluid is slick water with a viscosity of 0.002 Pa.s. The simulation was run in a mechanical step for 0.1 s to achieve an initial model equilibrium and continued in the fluid-solid coupling mode after starting the fluid step.

RESULTS AND ANALYSIS
Complex fracture geometries consist of tensile and shear fractures (Settgast et al., 2017). To quantify the simulation results of HF propagation in a laminated reservoir, the concept of the stimulated area (SA) is proposed to evaluate the fracture area generated by HF. The shear stimulated area (SSA) is defined as the region containing natural discontinuities that are subjected to shear slippage. By contrast, the tensile stimulated area (TSA) refers to the area where tension fractures form. These two parameters are used to describe the HF extensional patterns in laminated formations. The SSA may represent the area along the natural interface that is stimulated by HF. A higher SSA along an interface means larger fracture connectivity. On the other hand, the TSA may be more representative of the vertical extension of HF (i.e., fracture height). The greater the TSA in the vertical direction, the larger will be the fracture penetration. In the following sections, the effect of caprock Young's modulus, stress anisotropy, interface properties, and injecting fluid rate on SSA and STA, and hence, fracture geometry will be modeled and the results discussed.

Young's Modulus
The effect of Young's modulus on HF geometry is presented here. The caprock Young's modulus varied from 20 MPa to 27.7, 40, 50, and 60 MPa, while reservoir formation Young's modulus remained unchanged at 27.7 MPa. The front view (X-Y plane) and side view (Y-Z plane) of the simulation results are shown in Figures 3A,B. As shown in Figure 3A, the yellow wellbore is November 2021 | Volume 9 | Article 787736 5 oriented horizontally, two dark gray layers serve as the caprock of the model, HF is shown in blue propagating vertically, and some fluid flows into the natural fracture along the interface. The results illustrate that the caprock inhibits HF vertical propagation but promotes its extension along the interfaces between adjacent layers, particularly the inner interfaces near the injection point. In order to examine the effect of the angle of the approach, the models were run at inclinations of 0°, 10°, and 20°. The results show the greater the inclination, the larger will be the slippage area. This is in agreement with previously published results (Goldstein and Osipenko, 2015). Figure 4 presents the TSA and SSA corresponding to the results of Figure 3. The total SA is the sum of the TSA and SSA. The results of this figure show that TSA and SSA values are increasing in lockstep as the caprock Young's modulus increases. However, the change in inclination for TSA and SSA shows the opposite trend. Larger inclinations correspond to higher SSAs but lower TSAs. This is due to the fact that when the HF intersects the natural interface orthogonally (i.e., the layers' inclination is 0°), as opposed to the lower angle of approaches, the crossing interaction mode is more favorable to occur. This finding agrees with previous research works ( Gu and Weng, 2010).
When the caprock's stiffness is increased, an increase in the TSA and SSA is observed. To elucidate this conclusion further, a detailed analysis of the HF aperture was performed, and the corresponding results are shown in Figures 5A,B. The aperture profiles in the Y-Z plane for cases involving caprocks with a Young's modulus of 20 and 60 MPa are shown in this figure. Figure 5A shows a smooth decline of the aperture along the Y axis, whereas an obvious drop occurs at the interface between the reservoir and caprock layer. This is because Young's modulus is a mechanical property that quantifies the relationship between tensile stress and axial strain; thus, when a fracturing fluid with uniform fluid pressure transits from a high-to low-Young's modulus layer, the high-modulus rock will exhibit a relatively small strain. Additionally, to illustrate an increase in the  When natural interfaces in the same position are compared, the lower stiffness case exhibits a larger aperture value. Additionally, for a single case, one can observe that the interfaces between the reservoir and caprock layer exhibit a greater aperture than those farther from the injection point, particularly in the area contacted by HF interfaces.

Vertical Stress Anisotropy
To investigate the effect of vertical stress anisotropy (or differential stresses), the vertical stress (σ v ) was changed from 8 MPa to 10 MPa, 12 MPa, and 14 MPa, while the minimum horizontal stress remained constant at σ h 5 MPa. The front (X-Y plane) and side view (Y-Z plane) of the simulation results are shown in Figure 6. Figure 7 presents the TSA and SSA corresponding to Figure 6. From these figures, it is seen that as σ v increases, the total-SA decreases. Specifically, the SSA has the largest proportion of the contribution to the total SA decline compared to the TSA, implying that high-stress anisotropy results in less stimulation of the interface during HF operation. This finding is also reported by previous research, which claimed a higher differential stress causing a shorter offset in the interfaces (Chuprakov et al., 2010).

Tensile Strength
In this section, the impact of the tensile strength (T) of both caprocks and interfaces on HF propagation is studied by changes in the T values from 0.5 MPa to 3.5 MPA and 5.0 MPa. The front and side views of the models are depicted in Figures 8A,B. Figure 9 presents the TSA and SSA corresponding to Figure 8. A noticeable reduction in the total SA is observed when the tensile strength increases from 0.5 MPa to 5.0 MPa. This is contributed largely by reduction of the SSA from approximately 25 to 5 m 2 , while a little change is observed in the TSA. The induced stress of the fractures in the reservoir and the weak cementation of the natural fractures make the natural fracture faces susceptible to shear slip, thus creating the so-called hydraulic aperture to facilitate fluid leakage into the natural interfaces. Therefore, the magnitude of natural interface tensile strength is lower than that of the main formation; tensile normal stress on the natural fracture may cause its opening and activation in preference to local growth along its original path (Daneshy, 2019).  Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 787736

Injection Rate
The injection rate (Q) is one of the controllable operational parameters during the HF operation; therefore, it is important to understand its impact on fracture propagation. Five different models with variable injection rates of 0.02 m 3 /s, 0.04 m 3 /s, 0.06 m 3 /s, 0.08 m 3 /s, and 0.1 m 3 /s were considered while keeping other parameters the same for simulation purposes. Figures 10A,B show the results of the aperture as a function of the injection rate. It is seen that as the injection rate increases, the HF aperture in the reservoir increases, but the TSA decreases. The results in Figure 11 indicate that the total SA decreases as the fluid injection rate increases, and it is seen that this reduction is more contributed by the TSA and a little impact by the SSA. However, the general consensus regarding the effect of the injection rate is that an increase in the injection rate contributes to the tensile failure and results in more tensile fracture, whereas a low injection rate favors the formation of natural interface shear failure and act as a lubricant to reduce  Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 787736 8 friction and accommodate shearing (Beugelsdijk et al., 2000;Zou et al., 2016;Yu et al., 2019). This can be explained by Figures  12A,B, which illustrate that a high-injection rate case forms a high-pressure zone, with a breakdown pressure of 87.6 MPa, compared to 27.5 MPa, in a short time around the injection point and has a larger aperture in the reservoir layer, resulting in less fluid energy available to pressurize and extend the fracture in the vertical direction. Also, from the results of Figures 12C,D, one can observe that a high injection rate results in a greater interfacial extension and wider fractures than the case of lower injection rates. This is because the high-flow rate HF has higher hydraulic pressure at the intersection point (x −5 m in Figures  12C,D) when it contacts with natural interfaces than the low-rate case; hence, the HF prefers to propagate along the interface.

DISCUSSION
To assess that the previous results hold true universally, the considered parameters are transformed into a dimensionless space. The elastic modulus was defined as the ratio of the modulus in the caprock layer to the modulus of the reservoir, denoted as effective Young's modulus (E E ); the ratio of the tensile strength of the natural interface and caprock layer (which was assumed the same in this study) to the tensile strength of the reservoir formation was denoted as effective tensile strength (T E ); the effective stress anisotropy (S E ) was defined as the ratio of the vertical to the minimum horizontal principal stresses; and the effective injection rate (Q E ) was defined as the ratio of the injection rate to the minimum value used in each case. All of these dimensionless parameters were calculated to see if they present a meaningful trend for the quantitative analysis of HF propagation and its geometry.
The stimulated area ratio (SAR), defined as the ratio of the TSA to the SSA, was also used here. These larger SARs indicate more penetrability of the fracture versus connectivity or larger height (vertical extension) as opposed to horizontal extension. The relationship between effective Young's modulus and the SAR is presented in Figure 13A, which shows minor changes to the SAR as the stiffness of the caprock increases. This may be explained that fluid pressure restricted by the stiff caprock boosted both the TSA and SSA. Figure 13B shows the change in the SAR as a function of the effective stress anisotropy. It is seen that the SAR increases as stress anisotropy increases, a result which is consistent with the previous findings of this study which showed that high maximum horizontal stress reduces the SSA along the natural interfaces. The relationship between the SAR and the effective tensile strength is presented in Figure 13C. The results of this figure show a significant increase in the SAR when tensile strength increases. This is due to the fact that high tensile strength prevents further opening of the fractures formed by the interfacial slip. Figure 13D illustrates the effect of the injection rate on the SAR, indicating that there is a general tendency for the SAR to decrease as the injection rate increases, with slight variations due to changes in inclinations. In a nutshell, the previous results suggest that the larger Young's modulus and injection rate favor horizontal propagation of HF ( Figure 14B), whereas larger stress anisotropy and tensile strength favor vertical propagation of HF ( Figure 14A).
The evolution of HF from vertical to horizontal propagation modes can be visualized schematically in Figure 14D, where the bottom left and top right vertexes represent the vertical hydraulic fracture (VHF) and horizontal hydraulic fracture (HHF) propagation, respectively. In this figure, the bottom right and top left vertices represent the tension and shear stimulated areas, respectively. The transition between VHF and HHF propagation regimes is determined by the SAR. This means that the tensile strength and stress anisotropy enable HF to evolve from HHF to VHF via the path close to the SSA, whereas the injection rate Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 787736 evolves from the VHF to HHF via the path along the TSA. Furthermore, Young's modulus has an effect on the propagation mode through the TSA and SSA, causing it to follow a nearly diagonal path along the rectangle. Further investigation is needed to better understand the relationships between the SAR, the total stimulated area, and the HF network (HFN).

CONCLUSION
Based on the lattice-based numerical simulation of HF propagation in a laminated model presented in this work, the following conclusions were drawn: 1) Higher caprock Young's modulus promotes HF vertical and horizontal extensions. Higher vertical stress anisotropy inhibits HF horizontal extension. The higher tensile strength of the interface and caprock prevents HF horizontal propagation. A higher injection rate promotes HF width growth in the reservoir. 2) The tension stimulated area (TSA) and shear stimulated area (SSA) were used to determine the vertical and horizontal extendibility of HF. The TSA shows HF's ability to penetrate the interface and propagate vertically, whereas the SSA represents the ability of HF to form interface slippage or horizontal propagation of the fracture. These two parameters can quantitatively describe the behavior of HF propagation when it interacts with natural interfaces.
3) The ratio of the TSA to SSA, known as the stimulated area ratio (SAR), was used as a parameter to show how HF propagation transition occurs between vertical and horizontal directions. The greater the SAR, the greater will be the HF penetrability. A dimensionless space was used to demonstrate how mechanical and operational factors affect the propagation of HFs along different pathways in a laminated formation.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
DQ contributed to writing-original draft-and Methodology. JZ assisted with project administration and writing-review and editing. YL contributed to writing-review and editing. JL contributed to writing-review and editing. MR contributed to writing-review and editing. VR helped with supervision and writing-review and editing. BD contributed to resources, Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 787736 software, and writing-review and editing. RH helped with writing-review and editing.