#### Abstract

Although hydraulic fracturing has been practiced all over the world, the research on how the fracture height develops in time and space still leaves some missing gaps. The fracture height has been considered in most cases equal to the pay zone thickness, and the influence of temperature in this process has been omitted. Therefore, the aim of this paper is to study the effect of temperature, rock mechanical properties, and fluid injection rate on the development of the fracture geometry, especially on the fracture height. A multiphysics model was implemented using cohesive elements in a finite element model generated with equations in fracture mechanics. Once the model was calibrated with experimental data, it was used to conduct sensitivity studies to reveal the influence of main contributed factors such as the properties of rocks and fluids used in hydraulic fracturing, the injection rate of fracturing liquid, and especially the influence of temperature because this last aspect was omitted in literature review from previous studies. The results indicated that the fracture height depended strongly on the rock properties, not only the rock in the pay zone but also the ones in the adjacent layers. Besides, the influence of the fluid injection rate on the fracturing height is so great that it overwhelms the influence of temperature and mechanical parameters. Moreover, the impact of the leak-off coefficient is much less remarkable than that of the fluid viscosity, which demonstrates why in reality it is important to control the viscosity to achieve desirable results. This study can be applied in real life problems to predict fracture’s geometry generated in well stimulations.

#### 1. Introduction

Hydraulic fracturing is a popular well stimulation method used to increase well performance by improving the permeability. Due to the high cost of this operation, including drilling costs, proppants, and fluid injection, it is essential to be able to accurately predict the geometry of fractures generated by this process, which will lead to predicting correctly the effectiveness of the job. For this purpose, modeling the fractures’ development during hydraulic fracturing is an important tool to control the fractures’ parameters, such as the width, length, and height. Controlling the fracture’s height and length is useful not only in petroleum engineering but also in other domains such as underground storage of carbon dioxide, geophysics, and geothermal energy production stimulation ([1, 2]).

Hydraulic fracturing typically involves four key coupling processes: (1) rock deformation caused by fluid pressure on the fracture edges; (2) viscous fluid flow within the fracture; (3) fracture propagation in rock; and (4) fluid leakage from the fracture into the rock formation. The modeling of hydraulic fracturing, which has been done in the literature, has tried to cover these processes as much as possible. Some numerical modeling methods for hydraulic fracturing have been developed, including the finite element method (FEM), the discrete element method (DEM), the extended finite element method (XFEM), the numerical manifold method (NMM), the boundary element method (BEM), and the phase field method (PFM). Advani et al. [3] developed a software based on the finite element method for 3D simulation of hydraulic fracturing in multilayer reservoirs, with a focus on the propagation of hydraulic fracturing forces in elastic reservoirs’ classification. Lian-Jun et al. [4] developed a 3D hydraulic fracturing model that acknowledges the occurrence of radial flow, increasing the accuracy of the hydraulic fracturing height predictions. Lecampion [5] investigated the extended finite element method to model a hydraulic fracture in an impermeable medium. The author proposed the use of special enrichment functions at the fracture tip to capture the aperture and pressure singularities. Simoni and Secchi [6] used a remeshing algorithm and a staggered solving algorithm to cohesive fracture propagation within a nonhomogeneous porous under fluid pressure. Rho et al. [7] observed that the interface strength directly affected the fracture height’s growth as well as the fluid efficiency. These findings are important for a proper assessment of fracture height growth and better predictions of well production.

There are three typical fracture propagation mechanisms: linear elastic fracture mechanics (LEFM), elastic-plastic fracture mechanics (EPFM), and cohesive zone method (CZM). Linear elastic fracture mechanics (LEFM) was used with success to describe the mechanics of fractures in the brittle rock but should not be applied to the ductile rock. However, according to Atkinson [8], many brittle rocks exhibit ductility under high temperature and high pressure, so LEFM is not suitable for these problems. The CZM which included the cohesive elements in a finite element model was originally proposed by Barenblatt [9] and Dugdale [10] in the 1960s, but it gained interest from researchers just recently for its application in fracture stimulation. Chen et al. [11] adopted an interface element governed by a cohesive law to model the fracture propagation in an impermeable medium and compared its results with analytical solutions for a case where the fracture process is limited by rock fracture toughness. Lobao et al. [12] modeled both fluid leak-off and fracture’s propagation using the zero-thickness element method to investigate the influence of rock’s plasticity. Alfano et al. [13] show the application of CZM to predict fracture initiation and propagation in various materials, including brittle materials and also ductile materials. Yao et al. [14] used CZM to simulate hydraulic fracturing in ductile and quasibrittle rock. Elices et al. [15] pointed out the advantages, limitations, and challenges of CZM. Haddad and Sepehrnoori [16] concluded that the preferred fracture propagation model for brittle rock is LEFM, while EPFM should be used for ductile rock, and the most suitable method for quasibrittle rock is CZM. As shale is considered a quasibrittle rock, we decided to use CZM to model the propagation of hydraulic fracture in this research. Yao [17] compared the results obtained from the cohesive zone method with the ones given by pseudo-3D and Perkins and Kern (PKN) models. Shin [18] built a 3D model with XFEM for simultaneous propagation of multiple fractures in the geomechanical model, and the author found that parameters such as fluid viscosity and flow rate were key parameters affecting the hydraulic fracture height. Haddad and Sepehrnoori [19] proposed a CZM model based on XFEM in ABAQUS to simulate fracture initiation and propagation. Seyed et al. [20] built a 3D model based on the cohesive zone method to investigate the influence of parameters such as Young’s modulus, Poisson’s ratio, fracture toughness of rock, and viscosity of the fracturing fluid on the fractures initiated during hydraulic fracturing. Liu et al. [21] mentioned the limitations of LEFM in determining the hydraulic fracture height and presented a new chart to predict the hydraulic fracture height using the extended finite element method based on fluid-solid coupling equations and rock fracture mechanics. Carrier and Granet [22] considered a problem of fluid-driven fracture propagating in a permeable poroelastic medium using also a zero-thickness finite element to model the fracture, with the fracture’s propagation being governed by a cohesive zone model.

The literature review clearly showed that the linear elastic fracture mechanics (LEFM) is not suitable for ductile or quasiductile rocks. In addition, the LEFM is limited to a certain range of stress differentials between the pay zone and the adjacent formations. Furthermore, as the hardness of rock depends strongly on the temperature, the influence of temperature is therefore important in hydraulic fracturing modeling, an aspect which has been omitted largely in the literature so far. In addition, the singularity at the crack tip region, which poses considerable challenges for numerical modeling in classic fracture mechanics, can be avoided using the cohesive zone model. For the reasons mentioned above, our paper aimed at solving these existing problems by using the cohesive zone method (CZM) to model the development of fracture’s height and to study the impact of temperature during the process.

#### 2. Methodology

The cohesive zone method (CZM) can simulate the cohesive elements model for the initial loading, the initiation of damage, and the propagation of damage leading to eventual failure in the material. Unlike the normal element, the cohesive element can be as small as zero before the load is applied. The propagation of the cracks is restricted along with the layer of cohesive elements. The advantage of this method is that it does not require an existing fracture. Before the initiation damage of the cohesive elements, the constitutive relation of CZM is linear elasticity, with the elastic behavior described in the following equation [23].

are the dimensionless strains in normal, first shear, and second shear directions, given in the following.

are, respectively, the displacements in normal, first shear direction, and second shear direction. is the cohesive element thickness. For 2D simulation, the component of the second shear direction does not exist. The damage in cohesive elements follows the stress traction-separation law, which demonstrates the relation between traction and separation displacement as shown in Figure 1.

The flow of fluids in hydraulic fracturing includes two flow patterns. Tangential flow is used to describe the fluid flow in fracture, and normal flow is used to describe the leak-off of fracturing fluid into the rock matrix as shown in Figure 2.

For the tangential flow described in Equation (3), is the tangential flow rate of the liquid (), is the tangential flow fluid pressure gradient (); is the fracture aperture (m), and is the tangential permeability and is calculated according to the following equation: where is the viscosity of the fracturing fluid (Pa s).

The governing equation of the normal flow is defined in the following equation: where is the flow rates into the top surface (), is the leak-off coefficient of the top surface (), is the pore pressure on the top surface (Pa), is the flow rate into the bottom surface, is leak-off coefficients of the bottom surface (), is the pore pressure on bottom surfaces, and is the fluid pressure in the fracture (Pa).

Fracture initiation criteria are applied in cases where no initial crack existed in the material. The process of fracturing begins when the stresses and/or strains satisfy certain damage initiation criteria. The fracture criteria are shown in Table 1.

For the MAXE, MAXS, QUADE, and QUADS criteria, the user can select between horizontal or vertical crack growth (in this study vertical growth, i.e., along the adhesive layers’ length). All the six aforementioned criteria are fulfilled, and damage initiates, when reaches unity.

Fracture evolution was defined as the energy required for the fracture to propagate after the material had been damaged. There are several criteria for evaluating the fracture propagation, among which two standards are often used: the power law and BK’s law (Benzeggagh-Kenane).

The power law [23] criterion states that failure under mixed-mode conditions is governed by a power law interaction of the energies required to cause failure in the individual (normal and two shears) modes. It is given in the following equation: with the mixed-mode fracture energy .

The criteria of BK’s law [24] were introduced for hydraulic fracturing when initiation fracturing occurs and is represented by the following equation: where represent the energy release rates of one normal and two shear directions. The index represents the critical energy release rate, and is the property constant of the material. Both normal stress and shear stress are affected by the hydraulic fracturing process. The components of normal stress and shear stress can be described as in the following equation: where is the stress components predicted by the elastic traction-separation behavior for the current strains without damage. Here, the indices , , and represent, respectively, the normal stress and the two shear stresses. The overall damage in the material is presented by which captures the combined effects of all active mechanisms. No damage occurs () at the start of the simulation. monotonically evolves from 0 to 1 upon further loading after the initiation of damage. The evolution of the damage variable, , has the following usual form shown in the following equation:

Here, and are the displacement components at complete damage and the effective displacement at the initiation of damage, and is the maximum effective displacement.

In brief, fully modeling the hydraulic fracturing process necessitates solving a coupled system of governing equations, which include elasticity equations that determine the relationship between the fracture opening and the fluid pressure, differential equations for fluid flow inside the fracture to the fracture opening, the fluid pressure gradient, a fracture propagation criterion that allows for quasistatic fracture growth, and an elasticity equation that determines the relationship between the fracture opening and the fluid pressure gradient. All of the workflow for the numerical modeling is summarized in Figure 3.

#### 3. Influence of the Temperature on Rock Mechanical Properties and on Liquid Viscosity in Hydraulic Fracturing

##### 3.1. Influence of Temperature on Liquid Viscosity

Viscosity varies with temperature due to two factors: cohesive force and molecular momentum transfer. In liquids, the cohesive force is the main influence on viscosity because the molecules exert an attractive force on each other through moving layers. Increasing temperature results in a decrease in viscosity because a higher temperature induces that the particles have greater thermal energy and overcome more easily the cohesive force binding them together. The results shown in Figure 4 which were obtained from the studies [25, 26] will be used in this research. In this paper, the fluid of hydraulic fracturing is assumed to be water.

##### 3.2. Influence of Temperature on Young’s Modulus

Young’s modulus depends on the temperature of the material because as the temperature of the material increases, the atomic vibrations in the crystal structure also increase. As this atomic vibration increases, the atomic spacing in the crystal also increases, and the atomic force decreases. This decrease in atomic force leads to a decrease in Young’s modulus of the material as shown in Figure 5. The results for sandstone were obtained from the study of Wu et al. [27], the ones for the granite were extracted from Chen et al. [28], the mudstone results were given by Zhang et al. [29], the limestone results by Mao et al. [30], and the marble results are given by Zhang et al. [31]. All these results are summarized and presented in Figure 5 and will constitute the basis for the input parameters and for the calibration of the simulation results in this research.

##### 3.3. Influence of Temperature on Poisson’s Ratio

Poisson’s ratio is the ratio of the relative lateral strain to the relative axial strain (in the direction of force application). Poisson’s ratio is affected by temperature as shown in Figure 6 which clearly indicates that Poisson’s ratio decreases with increasing temperature. The results presented in Figure 6 were extracted from [32].

##### 3.4. Influence of Temperature on Rock Fracture Toughness

Fracture toughness of materials which are exposed to high temperatures will be reduced, which will weaken the rock’s ability to resist fracturing. Figure 7 shows that fracture toughness is not significantly affected at the temperature less than 200°C and begins to decrease sharply at conditions more than 200°C. The results of Westerly granite were extracted from Atkinson et al. [33], and the sandstone’s results are given by Feng et al. [34], black gabbro and Westerly granite are given by Meredith and Atkinson [35], stripa granite, FS gabbro, and FS marble are given by Zhang et al. [36], and BS granite is given by Wang [37].

##### 3.5. Influence of Temperature on Leak-off Coefficient

Equation (10) [38] shows that the leak-off coefficient is affected by the viscosity, which in turn is affected by the temperature. Hence, this equation will be used in this paper to model the development of fracture height, based on the results of the relation between viscosity and temperature stated in Section 3.1.

is the leak-off coefficient (), is the fluid compression coefficient (1/Pa), is the viscosity of fracturing fluid (Pa s), is the reservoir permeability (), is the pressure difference between a hydraulic fracture and rock matrix (Pa), and is the porosity. According to Salimzadeh and Khalili [39], if the leak-off coefficient is less than a certain value (), it will have little effect on the fracture height.

#### 4. Simulation Procedure

In this study, the hydraulic fracturing process will be modeled using a three-dimensional numerical model of the upper barrier–reservoir–lower barrier as shown in Figure 8. The size of this model is , with the thickness of the reservoir being 10 m, and the height of each upper and lower barrier is 20 m. The modeling took into account temperature change, rock dynamics, rock porosity, permeability, pore pressure, and fracture surface filtration.

To consider how hydraulic fracturing height is affected under the influence of temperature, in this study, a data set for calculation and simulation is used, presented in Table 2. This data set was extracted from a real-life hydraulic fracturing job.

The cases studied in this paper are summarized in Table 3. These case studies were proposed to effectuate a sensitivity study to understand the influence of temperature on the fracture’s height, as well as the influence of uncontrolled parameters such as rock’s properties and the influence of controlled parameters such as injection rate of fracturing liquid.

#### 5. Results and Discussions

##### 5.1. Case 0

In this basic case (Figure 9), we used input data mentioned in the previous section while ignoring the influence of the temperature in the modeling. This case will be used as a basic case for comparison with results from other cases. Furthermore, Figure 10 shows that this model was validated using a data set employed in earlier studies, and high accuracy was observed while comparing with the results extracted from [40]. For the units used in the results’ presentation in this paper, one time-step is equal to 100 seconds in real time, and the pressure unit is in Pa.

**(a)**

**(b)**

**(c)**

**(d)**

In the initial stage of hydraulic fracturing when the maximum tensile stress at the tip of the hydraulic fracturing is greater than the tensile strength of the rock, the fracture will be initiated. Firstly, the propagation of the fracture’s height is limited in the reservoir layer because the upper and lower barriers constrain the development of the fracture’s height between the interfaces (Figure 9(a)). After the hydraulic fracturing has reached the interface, the stress at the fracture tip continues to develop. Before the maximum principal stress reaches the tensile strength of the barrier rocks, the longitudinal propagation of the hydraulic fracture is restricted, which is manifested as width increase and unchanged hydraulic fracturing height (Figure 9(b)). The propagation speeds of hydraulic fracturing in the barriers zone are lower than that in the reservoir (Figure 9(c)), which can be explained by the difference in in situ stress. Teufel and Clark [41] found out that an increase of minimum horizontal stress could inhibit the vertical growth of the fracture into the barriers zone. This observation in Figure 9(c) can be confirmed again in Figure 9(d), where the existence of a restriction barrier presents a clear blocking effect on the vertical propagation of the hydraulic fracture, leading to a slow increase in the height of the fracture once it passes through the interfaces.

Hydraulic fracturing is a process that lasts from several hours to several days depending on the specification of each well. In this study, after trying simulation at 300 steps time (each step is equal to 100 seconds in real time), we observed that the stability of the fracturing fluid’s pressure that rendered the hydraulic fracturing height gradually stabilized, while the length and width of the fracture continue to propagate, as shown in Figure 11. In this study, we focused on the simulation of the fracture’s height; hence, a simulation time of the 50 steps time was chosen for this study because a clear change of height can be observed during that duration, while the stability of the fracture pressure can be obtained after that time as shown in Figure 12.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.2. Case Study 1

In this case, we studied the influence of mechanical parameters of the upper and lower barriers on the development of hydraulic fracturing height. The simulation results are presented in Figure 13.

**(a)**

**(b)**

**(c)**

**(d)**

From the simulation results, it is easy to see that the development of hydraulic fracture’s height is similar to the case zero when the fracture propagates to the interfaces (Figures 13(a) and 13(b)). The influence of Young’s modulus on the fracture’s height can be observed in Figures 13(c) and 13(d), where the results clearly showed that the fracture develops more rapidly in the zone with a higher Young modulus than in the zone with a lower Young modulus. This was the same observation in case 0. The results can partially be explained using the observations given in [42], where the fracture toughness of the materials is proportional to Young’s modulus. However, as the temperature increased in-depth, the ductility of the rock was also changed by temperature; hence, we needed another case to identify the impact of the temperature on the results, which will be treated in case 6 below.

##### 5.3. Case Study 2

The influence of viscosity was studied in this case to understand how the viscosity of the fracturing liquid affects the development of the fracture height, in order to complement the existing studies about the impacts of fracturing liquid’s viscosity on the fracture’s propagation [43].

The results presented in Figure 14(a) showed that when viscosity less than 0.1 Pa s is used for hydraulic fracturing, the longitudinal propagation of hydraulic fracture in the barrier zone is weak, which is beneficial to control the fracture height. However, as the fracture width is small, this is unfavorable to the movement of the proppant. Fracture width increases noticeably when the viscosity of fracturing fluid exceeds 5 Pa s (Figure 14(b)), which indicates that increasing the viscosity of fracturing fluid can improve fracture conductivity, and at the same time, fracture height also increases with the increase of fracturing fluid viscosity (Figure 14(c)). When the viscosity of fracturing fluid exceeds 20 Pa s, fracture height increases continuously, which indicates that the further increase of fracturing fluid viscosity will increase fracture propagation capacity in the barrier (Figure 14(d)). These results are interesting for application in real life because they showed that we can monitor a fracture’s height by controlling the fracturing liquid’s viscosity, and we do not want to increase viscosity unlimitedly because the fracture’s height can be out of control and reach into water or gas zones; we do not either use a very low fluid’s viscosity because of the unfavorable insert of the proppant.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.4. Case Study 3

In this case, we studied the impact of leak-off coefficient on fracture height. The results in Figure 15 indicate that the fracture height is influenced by the leak-off coefficient with permeability smaller than , although this influence is not significant (Figure 15(a)). When the permeability exceeds , the fracture height begins to slow down dramatically (Figure 15(b)). The fracture’s height is significantly affected when the permeability is greater than (Figure 15(c)), because of the fracturing permeability into the rock matrix upper and lower barriers, causing the fracturing fluid’s pressure to be insufficient to break the rock’s fracture toughness, resulting in no formation of crack height or width. When the permeability is greater than , no creation of fracture’s height can be observed (Figure 15(d)). These results showed that we need to be able to keep the permeability as low as possible in real-life applications because the results indicated clearly that the larger the permeability is, the higher its effect on the fracture dimension is. However, when the permeability reaches a certain limit, no more fracture height can occur. The conclusion here is that leak-off coefficient has little effect on the fracture height when the permeability of the rock bed is less than , which confirmed the finding of Salimzadeh and Khalili [39].

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.5. Case Study 4

The previous cases treated the influence of uncontrollable parameters (rock’s mechanical properties), so in this case, we studied the impact of the controllable parameter which is the injection rate of the fracturing liquid. It is noted that in this case, the influence of temperature is ignored.

The results in Figure 16 show that when the injection rate is increased, the pressure of hydraulic fracturing increases rapidly, forming high pressure in a short time, leading to easy penetration of the fracture into the upper and lower barriers. This result is confirmed in Figure 17 where we observed a higher fracturing pressure overtime at a higher injection rate.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.6. Case Study 5

As shown in the case study 2, increasing the viscosity causes the fracture height to increase. In case study 3, on the other hand, after the permeability reaches a certain limit, it will have no effect on the growth of fracture height. Hence, in this case, we will study the impact on the fracture height if the reservoir has both a high permeability and a high viscosity, in order to better understand the link between these two parameters as shown in Equation (10). The results of this case study can be used in practice for reservoirs with high permeability.

Figure 18(a) shows that the growth of fracture height can be noticed in the early stages of the simulation. This result indicated that when both the viscosity and leak-off coefficient are increased at the same time, the influence of viscosity on the fracturing height is so great that it overwhelms the influence of the leak-off coefficient, resulting in a development of hydraulic fracturing height, which can be confirmed again in Figures 18(b) and 18(c). Furthermore, viscosity has such a strong influence that the height rises above the upper restrictions (Figure 18(d)). These results are useful in practice because it showed that we should focus more on the fracturing fluid’s viscosity, a parameter that we can control easily in real life, than the leak-off coefficient, a parameter that is more difficult to control as it depends on many uncontrollable factors.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.7. Case Study 6

In this case 6, we wanted to clarify the influence of temperature on the development of hydraulic fracturing height by varying the temperature from 100°C to 250°C along the depth, while the mechanical parameters of layers were the same as in case 0.

Figure 19 indicates clearly that for a range of temperature from 100°C to 250°C, the influence of temperature on fracture’s height is not significant; that is why we observed the same results as in the original case (case 0), as shown in Figure 20. These results confirmed that the influence of mechanical parameters on the fracture’s height is much more considerable than that of the temperature. This can be explained by the fact that the mechanical parameters are not affected by temperature much in the temperature range from 100°C to 250°C, as analyzed in Section 3. Although this range of temperature is normally encountered in reality [44], some reservoirs can manifest higher temperatures, especially in the case of the high-pressure high-temperature reservoirs (HPHT); therefore, case 8 will study the influence of high temperature on the results.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.8. Case Study 7

This case study continues the research from previous cases because now we include the effect of temperature in the investigation. Now, the viscosity was initially taken to be equal to 20 Pa s, and its value will change in function of the temperature. A comparison between results drawn from Figures 21 and 14 presented in Figure 22 shows that the viscosity is indeed influenced by the temperature, which consequently leads to slower growth of hydraulic fracture’s height. A more detailed analysis of Figures 21(d) and 14(d) indicates a substantial difference in fluid pressure between the two cases. This indicates that as the temperature rises, the viscosity decreases, thus lowering the hydraulic fracturing’s height.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.9. Case Study 8

In this case, the temperature is increased above 500°C to verify the impact of high temperature on the development of hydraulic fracturing heights. Previous studies have shown that at temperatures above 500°C, the rocks will be affected greatly by temperature’s variation ([45, 46]).

At a high temperature, Figure 23 shows that the hydraulic fracturing height starts to slow down, accompanied with a decrease in the pressure of hydraulic fracturing in comparison with the above cases. It is then deduced that the high temperature influences greatly the development of fracture’s height, and the effect can be observed at different times. This observation can be explained as follows. As the temperature increases, the thermal motion of mineral particles, crystals, and atoms in the sandstone gradually increases. Consequently, the number of weak mineral crystal faces inside the rock increases, and the potential for local plastic deformation increases. These factors result in the formation of river-like and tear-like fracture patterns. As the temperature increases highly enough, the internal structure of the rock gradually deteriorates, and the number of fractures increases. This is an important cause for the deterioration of the mechanical properties of the rock. However, it is unlikely to encounter such high temperatures; in reality, the contribution in this case study 8 is, therefore, more academic than practice.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.10. Case Study 9

This case 9 is a complement study for case 4 because the influence of temperature was taken into account in this case 9. Therefore, the temperature varied from 100°C to 250°C along the depth of the layers, and the injection rate was increased.

The results in Figure 24 showed quite similar results of the hydraulic fracturing height as in the case 4. This indicated that the influence of temperature on the fracture’s height is dominated by the influence of the injection rate. Therefore, while running hydraulic fracturing in high temperature reservoirs, it is possible to consider using an increase in the injection rate to limit the influence of temperature on the height of the fracture.

**(a)**

**(b)**

**(c)**

**(d)**

#### 6. Conclusions

The cohesive-zone-finite-element-based model was proposed to simulate hydraulic fracture. Hence, it is possible to predict hydraulic fracturing height. Once the model was calibrated with experimental data, it was used to conduct sensitivity studies to reveal that the influence of the main contributed factors, therefore, brought academic contribution as well as application aspect because the modeling workflow and the analysis of results were given thoroughly and can be applied in reality to enhance the effectiveness of hydraulic fracturing.

Firstly, this study showed that the hydraulic fracturing height is clearly affected by temperature, an aspect usually omitted in this domain of research and also in reality. The higher the temperature is, the higher its impact on the fracture’s height is.

Secondly, mechanical parameters also have a clear impact on the propagation of the hydraulic fracture. The influence of mechanical properties on the fracture’s height is more considerable than that of the temperature.

Thirdly, the study also clearly demonstrated the relationship between the fracturing fluid’s viscosity and the leak-off coefficient, as well as their influence on the development of fracture height. Some interesting conclusions can be deduced for practice in real life. Firstly, we need to keep the leak-off coefficient as low as possible to achieve higher performance in hydraulic fracturing, and when the leak-off coefficient reaches a certain limit, no more fracture’s height can occur. Moreover, the influence of fracturing fluid’s viscosity on the fracturing height is so great that it overwhelms the influence of the leak-off coefficient; hence, in practice, we should pay more attention on the fracturing fluid’s viscosity, a parameter that can be controlled easily, than the leak-off coefficient, a parameter which is harder to control.

Last but not least, this study showed that increasing the injection rate will have a remarkable effect on the propagation rate and amplitude of the hydraulic fracture’s height. The influence of the injection rate on the fracturing height is so great that it overwhelms the influence of temperature and mechanical parameters. This shows that, in reality, in addition to the mastering of information about the mechanical parameters of the rock, controlling the properties of the fracturing fluid and the injection rate constitute the most important factor to control the desired height of hydraulic fracture.

#### Abbreviations

BEM: | Boundary element method |

BK: | Benzeggagh-Kenane law |

CZM: | Cohesive zone method |

DEM: | Discrete element method |

EPFM: | Elastic-plastic fracture mechanics |

FEM: | Finite element method |

HPHT: | High-pressure high-temperature |

KGD: | Khristianovich-Geertsma-de Klerk fracture model |

LEFM: | Elastic fracture mechanics |

MAXE: | Maximum nominal strain damage |

MAXPE: | Maximum principal strain damage |

MAXPS: | Maximum principal stress damage |

MAXS: | Maximum nominal stress damage |

NMM: | Numerical manifold method |

PFM: | Phase field method |

PKN: | Perkins-Kern-Nordgren fracture model |

QUADE: | Quadratic nominal strain damage |

QUADS: | Quadratic nominal stress damage |

XFEM: | Extended finite element method |

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.