A Model for Multiple Hydraulic Fracture Propagation with Thermo-Hydro-Mechanical Coupling Effects

: In this study, a coupled thermo-hydro-mechanical model to simulate multiple hydraulic fracture propagation is presented. Fracture propagation with elastic deformation is described by using a displacement discontinuity method. The temperature distribution and induced thermal stress are calculated via a semi-analytical method in an explicit way. An iterative scheme is proposed to solve the coupling between fracture propagation with ﬂuid ﬂow and induced thermal stress. The numerical model is validated against related analytical solutions. Several numerical cases are modeled to investigate the controlling factors for uniform growth of multiple fractures. Results show that using non-uniform fracture spacings and proper increasing the spacing for fractures away from the heel of wellbore promote the uniform growth of multiple fractures by comparison with using uniform fracture spacings. Increasing the perforation diameter for the middle cluster also works. Besides, single-wing fracturing could greatly improve the uniform growth of multiple hydraulic fractures. Finally, it shows that the thermal stress has a signiﬁcant inﬂuence on fracture geometrical size but has limited effect on fracture propagation path. In addition, the thermal effect promotes the uniform growth of multiple fractures. If the variance is small, it means that the half-length of fracture is approximated to each other, and the uniform growth of multiple fractures is obtained.


Introduction
The development of hydraulic fracturing is explored in the petroleum and natural gas industry. The technology and theory of it have been investigated quite intensively during past decades. With the exploration and development of unconventional oil and gas, multi-cluster fracturing with horizontal well is becoming more and more popular. A major concern is to predict how multiple hydraulic propagates from horizontal wells. In general, the thermal effect is neglected in the modeling of hydraulic propagation since the formation temperature is not high in shallow layers. However, this is not appropriate for deep reservoirs because the temperature is much higher than that of shallow layers. The great gap of temperature between rock and fracturing fluid would give rise to extra stresses on the fracture surface, namely thermally-induced stresses. It is well known that the effect of stress shadowing plays a significant role in the simultaneous propagation of multiple fractures but the interaction between stress shadowing and thermal induced stresses on their propagation is not clearly understood so far. Therefore, it is meaningful to investigate the thermal effect on multiple fractures propagation.
Many numerical methods are presented to simulate hydraulic fracture propagation. The first kind is the finite element method. The finite element method is usually combined with the cohesive zone method to characterize fracture initiation and propagation [1][2][3], which avoids remeshing during fracture propagation, but the computational cost is high. The extended finite element method allows the crack to propagate across finite elements without remeshing, which is gradually applied to modeling hydraulic fracturing [4][5][6][7]. Zeng et al. [8] proposed a hybrid approach combining the extended finite element method and an embedded discrete fracture model to investigate the poro-elastic effect on hydraulic fracture propagation. However, when it comes to the propagation of complex fractures, the usage of the extended finite element method is still cumbersome. The second kind of method is boundary element method. The widely used one is displacement discontinuity method [9,10]. Only boundaries need to be discretized, which results in less computational cost than the finite element method. Cheng et al. [11] used the displacement discontinuity method to model simultaneous and sequential fracturing in a horizontal well. Tang et al. [12] presented a fully three-dimensional displacement discontinuity method to investigate multiple hydraulic fractures propagation in heterogeneous stress field. The third kind is discrete element method. The discrete element method is initially used to study material deformation and failure under external loading and extended to model hydraulic fracture propagation by setting fluid pressure element between particle elements [13,14]. Zou et al. [15] presented the use of discrete element modeling to investigate the effect of bedding plane on hydraulic fracture height growth. In addition, the combination of finite-discrete element method is proposed, which has the advantages of both methods and shows great promise for hydraulic fracturing modeling [16,17]. The common characteristic of these methods is that fractures are explicitly represented, which results in strong displacement discontinuity perpendicular to fractures when using these methods. Another different manner to represent fractures is proposed by using additional variables. For example, the phase field method and peri-dynamics method have been presented to resolve hydraulic fracture propagation in recent years [18][19][20]. What they have in common is that they are based on the minimization of energy, and the advantage over other methods is that fracture initiation, propagation and coalescence can be automatically captured without additional criterions. However, the discretized meshes are so refined that the computational cost is high, which makes them far from practical applications.
The effect of thermal stress on hydraulic fracturing has received lots of attention. Kim et al. [21] developed a coupled flow-thermal-geomechanical simulator hydraulic fracturing process. Li et al. [22] presented a fully coupled thermo-hydro-mechanical model for hydraulic fracture propagation. Tomac and Gutierrez [23] investigated the THM coupling effects on hydraulic fracturing using a bonded particle model. In addition, rock may exhibit plastic deformation and failure in deep reservoirs. Liu et al. [24] simulated the elasto-plastic hydraulic fracture propagation in deep reservoirs coupled with thermal effect. Zeng et al. [25] proposed an extended finite element method for hydraulic fracturing considering the thermo-hydro-elastic-plastic coupling. These studies are based on volumetric discretization, which results in high computational cost of coupling solution. Even though the thermal effect was considered, the thermal effect on multiple fractures propagation was not clearly addressed. Besides, the displacement discontinuity method were used to simulated multiple fracture propagation [26,27], but only fluid pressure is considered acting on fracture surface. Here, the additional thermal stress is imposed on fracture surface, which is induced by the temperature difference between injected fluid and reservoir rock. And a new solution scheme based on displacement discontinuity method is proposed to investigate the thermal effect on multiple hydraulic fractures propagation.
The structure of the paper is as follows: Section 2 contains the governing equations of multiple hydraulic fractures propagation in a horizontal well, including the stress distribution due to displacement discontinuity in rock, the equation of fluid flow in the fractures and wellbore, the solution of thermal induced stress as well as fracture propagation criterion. Section 3 presents the discretized form of the equations of fluid pressure and thermal Energies 2021, 14, 894 3 of 23 stress and the iterative scheme for thermo-hydro-mechanical coupling solution of fracture propagation model. Section 4 validates the numerical method via comparing with several analytical solutions, and then presents the detailed analysis for the effects of fracture spacings (uniform and non-uniform), non-uniform perforation parameters, single-wing treatment and thermal stress on hydraulic fracture propagation. Section 5 addresses some conclusions.

Model Formulations
The model of multiple hydraulic fractures propagation in a horizontal well is established. Four physical processes are considered, fracture propagation with rock deformation, fluid flow in the fractures and wellbore, thermal induced stress due to temperature difference between fracturing fluid and reservoirs rock, and the fracture propagation criterion. Each part is described in detail as follows.

Stress and Displacement Discontinuity in Rock
Hydraulic fracture propagation is a typical problem of displacement discontinuity. According to Crouch and Starfield [28], the analytical solutions of induced stresses due to a constant displacement discontinuity are given as: where σ xx and σ yy are the normal stresses along x, y directions, and τ xy is the shear stress. D x and D y are the displacement discontinuities along x, y directions. G is the shear modulus. f is a function of position x and y, which is expressed as: where ν represents the Poisson's ratio. Equations (1)-(3) give the expressions of stresses for a straight line of displacement discontinuity. When a curved line of displacement discontinuity is considered, the line can be discretized into finite straight segments with constant displacement discontinuity. The induced stresses are the superposition of contribution from each segment. For each segment of displacement discontinuity, the stresses expressed in Equations (1)-(3) are in the local coordinate system, which can be expressed in the global system through coordinate transformation. Therefore, the normal and shear stresses along a segment can be given as: where N is the number of segments. σ n and D n are the normal stress and displacement discontinuity. σ s and D s are the shear stresses and displacement discontinuity. C ij is the influencing matrix that accounts for contribution of segment j on segment i, which can be obtained from Equations (1)- (3). G ij is a 3D correction factor [29], which is given as: where d ij is the distance from segment i to segment j. h is fracture height. α and β are two empirical coefficients, which are usually set to be 2 and 3, respectively. As for the problem of hydraulic fracturing, the normal stress on the fracture is superposition of remote in-situ stress and fluid pressure, and the shear stress is only due to remote in-situ stress, which is written as: where θ i is the angle between segment i and the positive direction of x axis. Given the normal and shear stresses along each line segment, the displacement discontinuities D n and D s on each line segment are taken as unknows with the total number of 2N, and the same number of linear equations (Equations (5) and (6)) can be obtained. Therefore, the displacement discontinuities of each line segment can be solved, which is known as the displacement discontinuity method (DDM). The normal displacement discontinuity is the width of hydraulic fracture, which determines the fracture conductivity in the fluid flow as described below.

Fluid Pressure in Fractures and Wellbore
Fracturing fluid is injected from ground and flows into each perforation through wellbore. Due to the effect of stress shadowing, the injected fluid is not evenly distributed into each perforation, as shown in Figure 1. This is the main difference between single fracture propagation and multiple fractures propagation. In addition, the temperature distribution of each fracture is also different from each other. Two processes need to be captured, including fluid flow in fractures and wellbore.
where dij is the distance from segment i to segment j. h is fracture height. α and β are two empirical coefficients, which are usually set to be 2 and 3, respectively. As for the problem of hydraulic fracturing, the normal stress on the fracture is superposition of remote in-situ stress and fluid pressure, and the shear stress is only due to remote in-situ stress, which is written as: 2 2 sin cos where θi is the angle between segment i and the positive direction of x axis. Given the normal and shear stresses along each line segment, the displacement discontinuities Dn and Ds on each line segment are taken as unknows with the total number of 2N, and the same number of linear equations (Equations (5) and (6)) can be obtained. Therefore, the displacement discontinuities of each line segment can be solved, which is known as the displacement discontinuity method (DDM). The normal displacement discontinuity is the width of hydraulic fracture, which determines the fracture conductivity in the fluid flow as described below.

Fluid Pressure in Fractures and Wellbore
Fracturing fluid is injected from ground and flows into each perforation through wellbore. Due to the effect of stress shadowing, the injected fluid is not evenly distributed into each perforation, as shown in Figure 1. This is the main difference between single fracture propagation and multiple fractures propagation. In addition, the temperature distribution of each fracture is also different from each other. Two processes need to be captured, including fluid flow in fractures and wellbore.  Considering the fracturing fluid as Newtonian fluid, the fluid flow in the fracture is modeled as Poiseuille flow, and the equation of motion is given by the lubrication equation as: where s and t are the space and time variables. µ is the fluid viscosity. q is the fluid flux through cross section. w and p are the fracture width and fluid pressure. Assuming the fracturing fluid as non-compressible fluid, the conservation equation is written as: In the above Equation (11), the first term accounts for the net flow out, and the second term accounts for the variation of cross section area, and the third term is the leakoff from fracture into formation. The leakoff term is given by Carter model as: where C L is the leak-off coefficient, and τ 0 is the time first exposed to fluid. As mentioned above, the fluid received by each wing of fractures is not the same. Therefore, fluid flux distribution should be modeled. The first equation is the conservation of injection rate, which is written as: where Q T is the pre-given injection rate. m is the number of fractures. Q 2k−1 and Q 2k are the injection rates into two wings of fracture k. If there is only one wing of fracture, only one injection rate is needed. According to the total injection rate, the global conservation equation is written as: Another equation is the pressure drop from wellbore to each fracture, and two kinds of pressure drop are induced by friction of wellbore and perforation. When there are two wings for each fracture, two equations of pressure drop are considered, which are given as: where p 0 is the pressure at the inlet of wellbore. p w,2k−1 and p w,2k are fluid pressures at the inlet of two wings of fracture k. ∆p cf and ∆p pf are the pressure drops induced by friction of wellbore and perforation. As for the frictional pressure drop of wellbore, it is expressed as: where x j is the distance of fracture j to the inlet of wellbore, and D is the diameter of wellbore.
As for the frictional pressure drop of perforation, it is given as: where ρ sl and ρ w are the densities of slurry and water. n p,k and d p,k are the number and diameter of perforation of fracture k.

Thermal Induced Stress
According to Ghassmi and Zhang [30], the fundamental solutions for induced stresses due to a point heat source with unit heat strength in infinite medium are given as follows [31]: where λ s and α are the thermal conductivity and volumetric thermal expansion coefficient of rock. E is the Young's modulus.
In the above Equations (20)- (22), r 2 and ξ 2 are defined as follows: where c is the thermal diffusivity of rock, which is given as: where c s and ρ s are the specific heat capacity and density of rock. In addition, Ei(u) is the exponential integral function, defined as: The temperature of the fluid is lower than that of rock, and the fluid would be heated. Whitsitt and Dysart [32] presented an analytical solution for the fluid temperature distribution in the fracture. Here, the analytical solution is adopted to calculate the fluid temperature as: where T 0 is the initial temperature of reservoir and T w is the temperature of injected fracturing fluid. Other function and variables are defined as where ρ f and c f are the density and specific heat capacity of fracturing fluid. ϕ is the porosity of reservoir. L is the length of one fracture wing, and Q is the corresponding fluid flux into the wing.

Fracture Propagation Criterion
It is assumed that fractures propagate in elastic condition, and fractures tend to extend in the direction parallel to the maximum circumferential stress direction, which is known as the criterion of maximum circumferential stress. Considering the fracture of mixing mode of I and II, the maximum circumferential stress is reached at the following angle: where θ is the polar angle at the fracture tip. K I and K II are the mode I and II stress intensity factors.
The equivalent stress intensity factor is calculated with K I and K II as It requires that the equivalent stress intensity factor exceeds fracture toughness of rock for the fracture to propagate. When multiple fractures propagate simultaneously, the growth step for each fracture wing can be given in an explicit way or retrieved in an implicit way. Here the explicit way is used and the growth step for each fracture wing is proportional to the fluid velocity into tip segment as: where ∆L is a preset growth step. v k is the fluid velocity into the tip segment of fracture wing k, and v max is the maximum value of fluid velocities into tip segments. From Equations (31) and (32), it needs to resolve the mode I and II stress intensity factors to determine fracture propagation. Based on the solution of displacement discontinuities, the stress intensity factors are calculated as [29]: where a is the half length of tip segment. D t n and D t s represent the normal and shear displacement discontinuities of tip segment.

Numerical Algorithm
To obtain solution of hydraulic fracture propagation, the governing equations presented above are solved numerically. The displacement discontinuity equations are written in discretized form as shown in Equations (5) and (6). Therefore, it remains to discretize the equations of pressure and temperature fields, and the iterative scheme is needed to obtain the final solution.

Pressure Field
The fluid pressure equation is discretized by the finite volume method. The integration form of conservation Equation (11) is given as: where q ij is the fluid flux from segment i to segment j, which is defined at the boundary. w i is the average width of segment i, and l i is the length of segment i. From the lubrication Equation (10), the fluid flux between adjacent segments is approximated as: where p i is the fluid pressure of segment i.

Thermal Stress Field
The fundamental solutions (Equations (20)- (22)) give the induced thermal stresses due to a constant point heat source with unit strength. As for hydraulic fracturing problem, the fractures are discretized into finite segments. The induced thermal stresses at the segment i (the center point) due to a constant line heat source segment j with unit strength can be computed by integrating the solutions of point heat source (Equations (20)-(22)) as: The integrations are evaluated numerically by further dividing segment j into finite sub-segments. The line heat source of every sub-segment is assumed to shrunk to a point source at the sub-segment center, whose strength is equal to the total strength of line heat source. These above equations give the temperature and stresses at segment i due to the fictitious line source of segment j. Besides, the stresses mentioned above are presented in the local coordinates of segment j, and they should be transformed into the global coordinate for further calculation. Therefore, the total stresses at segment i due to the fictitious line heat sources of all segments can be obtained by the accumulated effects of all line heat source as: where H j is the total heat strength of segment j. The regular subscripts n and s denote the normal and shear stresses. During hydraulic fracture propagation, the strength of heat source on each segment is not constant. A time marching scheme is needed to capture the variation of heat flow. It is assumed that the heat source can be divided into finite sub-sources, and each sub-source begins to take effect at different time. Then the total induced thermal stresses are the superposition effects of all sub-sources. For example, the induced stresses at segment i at time t n are given as: where the italic superscript n denotes time step t n . H m j is the sub-source of segment j starting to take effect at time t m .  [31,33]. As fluid flows from segment i to segment j, the temperature becomes T j from T i (shown in Figure 2). The segment j draws thermal energy from surrounding rock, and acts as a negative heat source. For duration ∆t (from t m to t m+1 ), the net volume entering segment j is written as: where q ij is the fluid flow rate from segment i to segment j, which is calculated according to Equation (37).
where the italic superscript n denotes time step tn. H m j is the sub-source of segment j starting to take effect at time tm.
It can be seen from Equations (43) and (44) [31,33]. As fluid flows from segment i to segment j, the temperature becomes Tj from Ti (shown in Figure 2). The segment j draws thermal energy from surrounding rock, and acts as a negative heat source. For duration Δt (from tm to tm+1), the net volume entering segment j is written as: where qij is the fluid flow rate from segment i to segment j, which is calculated according to Equation (37).
( ) The heat strength of segment j over the duration Δt (from tm to tm+1) is expressed as From Equation (47), the heat strength of arbitrary segment over arbitrary previous time step can be solved, and then the total induced thermal stresses for current time step are obtained (Equations (43) and (44)). Therefore, the total stresses on fracture segments (Equations (5) and (6)) are renewed accounting for the thermal stresses as: The heat strength of segment j over the duration ∆t (from t m to t m+1 ) is expressed as From Equation (47), the heat strength of arbitrary segment over arbitrary previous time step can be solved, and then the total induced thermal stresses for current time step are obtained (Equations (43) and (44)). Therefore, the total stresses on fracture segments (Equations (5) and (6)) are renewed accounting for the thermal stresses as:

Iterative Scheme
The whole process of hydraulic fracture propagation involves three physical fields, that are stress, pressure and temperature fields. Here, we adopt a strategy that stress (displacement discontinuity) and pressure are solved iteratively while the temperature and thermal induced stress are solved explicitly. The solving algorithm is illustrated as follows. For a new time step, each fracture tip is extended a length (Equation (33)) if the propagation criterion is satisfied. The temperatures on each fracture segment in previous time steps are solved (Equation (27)). The heat strengths of each fracture segment over previous time steps are calculated from equation (Equation (47)), and then the induced thermal stresses on each fracture segment can be obtained (Equations (43) and (44)). These thermal stresses are added to the total stresses on fracture segments (Equations (48) and (49)). After that, the fluid pressure and displacement discontinuities are solved iteratively. The unknowns include time step dt, fluid pressure p i , displacement discontinuities D i n D i s , injection rate Q k and fluid pressure p 0 . Considering that the equations between fluid pressure and displacement discontinuities are linear, we select the time step dt, fluid pressure p i , injection rate Q k and pressure p 0 as iteration variables. The number of iteration variables is N + 2m + 2. We need to construct the same number of equations for the use of Newton iterative solution of the problem. The Newton iteration method requires that the iteration variables are given initial values. Based on the given fluid pressure and already evaluated thermal stresses on all fracture segments, the displacement discontinuities are achieved by solving Equations (48) and (49), in which the normal displacement discontinuities give the fracture widths of all segments. N equations can be deduced by substituting Equation (37) into Equation (36). 2m equations can be deduced from Equations (15) and (16). And 2 equations are obtained from Equations (13) and (14). The number of equations set and the number of iteration variables are equal. Therefore, the selected iterative variables can be solved by using Newton iteration method. The convergence criterion requires that the relative errors of iterative variables between consecutive iterative steps are less than the tolerance ε. Therefore, the criterion for the iterative scheme is expressed as |dt where l represents the iterative step. The tolerance is set as 0.001 in the study. After the convergence, the fracture propagation criterion is checked for each fracture tip. As illustrated in Section 2.4, it requires that the equivalent stress intensity factor exceeds rock toughness for fracture extending forward. Therefore, the equivalent stress intensity factor is calculated (Equation (32)), and compared with the rock toughness. If there is a tip satisfying the criterion, the program directly proceeds to the next time step. Otherwise, if the criterion is not met, the iteration is repeated with renewed fluid pressure. The flow chart of iterative scheme for a time step is shown in Figure 3.

Discussion
In this section, the proposed numerical model is validated against the analytical solutions of two cases firstly. And then the effects of fracture spacings (uniform and nonuniform), non-uniform perforation parameters on hydraulic fractures propagation are analyzed. The result of single-wing fracturing is also compared with that of bi-wing fractur-

Discussion
In this section, the proposed numerical model is validated against the analytical solutions of two cases firstly. And then the effects of fracture spacings (uniform and non-uniform), non-uniform perforation parameters on hydraulic fractures propagation are analyzed. The result of single-wing fracturing is also compared with that of bi-wing fracturing. Finally, the effect of thermal stress is investigated.

Model Validation
The first case is a vertical fracture case as shown in Figure 4, in which the fracture is imposed a constant fluid pressure inside. The fluid pressure would induce extra stresses around the fracture, and the induced stresses at arbitrary point is given as [34]: The variables in above equations are defined as: We discretize the fracture into finite segments, and obtain the induced stresses using the displacement discontinuities method. The obtained results of dimensionless stresses are compared with that of the analytical solutions as presented in Figure 5. It shows great agreement between the numerical solution and analytical solution, which indicates that the stress and displacement discontinuity are correctly captured with the DDM.  The fracture in the first case is stationary, and the second case is to model fracture propagation. To validate the calculation of the hydraulic fracture propagation, the numer- We discretize the fracture into finite segments, and obtain the induced stresses using the displacement discontinuities method. The obtained results of dimensionless stresses are compared with that of the analytical solutions as presented in Figure 5. It shows great agreement between the numerical solution and analytical solution, which indicates that the stress and displacement discontinuity are correctly captured with the DDM.  The fracture in the first case is stationary, and the second case is to model fracture propagation. To validate the calculation of the hydraulic fracture propagation, the numer- The KGD model is analytically solved and expressed as [35]: where L is the fracture half-length. G is the shear modulus. q 0 is half of the injection rate. The hydraulic fracture half-length and width are solved numerically and analytically, which are shown in Figure 6. The results show that good consistencies are obtained. The KGD model is analytically solved and expressed as [35]: where L is the fracture half-length. G is the shear modulus. q0 is half of the injection rate. The hydraulic fracture half-length and width are solved numerically and analytically, which are shown in Figure 6. The results show that good consistencies are obtained. It should be mentioned that there are no thermal effects for the verification cases. The verification cases are intended to verify that the presented numerical method to the basic problem of fracture propagation is correct. The calculation method of thermal stress is referred to Shen's study [31], and it is proved to be effective in their study. Therefore, the thermal effect on hydraulic fracture propagation can be correctly captured.

Uniform Fracture Spacings
To  It should be mentioned that there are no thermal effects for the verification cases. The verification cases are intended to verify that the presented numerical method to the basic problem of fracture propagation is correct. The calculation method of thermal stress is referred to Shen's study [31], and it is proved to be effective in their study. Therefore, the thermal effect on hydraulic fracture propagation can be correctly captured. The fluid flux from wellbore into each fracture is not the same, and it varies with the injection time as shown in Figure 8. The left fracture closer to the heel of wellbore receives more fluid than the right fracture. As the fracture spacing increases, the difference of flux distribution between them decreases. The fluid flux from wellbore into each fracture is not the same, and it varies with the injection time as shown in Figure 8. The left fracture closer to the heel of wellbore receives more fluid than the right fracture. As the fracture spacing increases, the difference of flux distribution between them decreases.
Since the left fracture receives more fracturing fluid, the left fracture length is generally larger than right fracture length, which means that the length of each fracture is not the same. The variance of fracture half-length is calculated to investigate how fracture spacings affect the uniform growth of multiple hydraulic fractures as follows: where L i is the half-length of fracture i, and L is the average of fracture half-lengths. m is the number of fracture. The fluid flux from wellbore into each fracture is not the same, and it varies injection time as shown in Figure 8. The left fracture closer to the heel of wellbore more fluid than the right fracture. As the fracture spacing increases, the differenc distribution between them decreases. Since the left fracture receives more fracturing fluid, the left fracture length i ally larger than right fracture length, which means that the length of each fractu the same. The variance of fracture half-length is calculated to investigate how spacings affect the uniform growth of multiple hydraulic fractures as follows: The average half-lengths and corresponding variances for different fracture spacings with the same injection time is shown in Figure 9. As the fracture spacing increases, the average half-length increases slightly, but the variance decreases greatly. The variance measures the spread between fracture half-lengths. If the variance is small, it means that the half-length of fracture is approximated to each other, and the uniform growth of multiple fractures is obtained. The average half-lengths and corresponding variances for different fracture sp with the same injection time is shown in Figure 9. As the fracture spacing increas average half-length increases slightly, but the variance decreases greatly. The va measures the spread between fracture half-lengths. If the variance is small, it mea the half-length of fracture is approximated to each other, and the uniform growth o tiple fractures is obtained. The case of three-cluster fracturing is also conducted. The fracture propagation and widths distribution are shown in Figure 10. It can be seen that the length of m fracture is much smaller than those of the left and right fractures when the fractur ing is 20 m. As the fracture spacing increases to 40 m, the length of each fracture is equal. The fluid flux distribution for different spacings is presented in Figure 11. T fractures receive most of fracturing fluid. As the fracture spacing increases, the diff of flux distribution between each fracture decreases. The case of three-cluster fracturing is also conducted. The fracture propagation paths and widths distribution are shown in Figure 10. It can be seen that the length of middle fracture is much smaller than those of the left and right fractures when the fracture spacing is 20 m. As the fracture spacing increases to 40 m, the length of each fracture is nearly equal. The fluid flux distribution for different spacings is presented in Figure 11. The side fractures receive most of fracturing fluid. As the fracture spacing increases, the difference of flux distribution between each fracture decreases. and widths distribution are shown in Figure 10. It can be seen that the length of middle fracture is much smaller than those of the left and right fractures when the fracture spacing is 20 m. As the fracture spacing increases to 40 m, the length of each fracture is nearly equal. The fluid flux distribution for different spacings is presented in Figure 11. The side fractures receive most of fracturing fluid. As the fracture spacing increases, the difference of flux distribution between each fracture decreases.  Moreover, the average half-length and variance for different spacings is shown in Figure 12. The average half-length with fracture spacing 30 m is slightly larger than that with fracture spacing 40 m. The reason is that larger frictional pressure drop is produced in the wellbore when using larger fracture spacing. However, as the fracture spacing increases, the variance of fracture half-length decreases.  Moreover, the average half-length and variance for different spacings is sh Figure 12. The average half-length with fracture spacing 30 m is slightly larger th with fracture spacing 40 m. The reason is that larger frictional pressure drop is pr in the wellbore when using larger fracture spacing. However, as the fracture spa creases, the variance of fracture half-length decreases. Moreover, the average half-length and variance for different spacings is shown in Figure 12. The average half-length with fracture spacing 30 m is slightly larger than that with fracture spacing 40 m. The reason is that larger frictional pressure drop is produced in the wellbore when using larger fracture spacing. However, as the fracture spacing increases, the variance of fracture half-length decreases.
Moreover, the average half-length and variance for different spacings is sh Figure 12. The average half-length with fracture spacing 30 m is slightly larger th with fracture spacing 40 m. The reason is that larger frictional pressure drop is pro in the wellbore when using larger fracture spacing. However, as the fracture spac creases, the variance of fracture half-length decreases.   Figure 13. It can be seen that when using the non-uniform spacings of 30 m and 20 m, the left fracture has the length much larger than the middle one and the right one, and the middle fracture is nearly as long as the right fracture. However, when using the non-uniform spacings of 20 m and 30 m, the right fracture is the longest among the three fractures, and the length of it is slightly larger than that of the left fracture.  Figure 13. It can be seen that when using the non-uniform spacings of 30 m and 20 m, the left fracture has the length much larger than the middle one and the right one, and the middle fracture is nearly as long as the right fracture. However, when using the non-uniform spacings of 20 m and 30 m, the right fracture is the longest among the three fractures, and the length of it is slightly larger than that of the left fracture. To quantify the effect of non-uniform fracture spacings on the uniform growth of multiple fractures, the resulting variances for each case are compared as shown in Table 1. Little difference of the average half-length between each case is obtained. It can be shown that the variance of case 1 using the non-uniform spacings of 30 m and 20 m is slightly smaller than that of case 3 using the uniform spacing of 20 m, but is greatly larger than that of case 4 using the uniform spacing of 30 m. However, when using the non-uniform spacings of 20 m and 30 m, the variance drops greatly in comparison with case 3, but is still larger than that of case 4. It indicates that proper increase of spacing for fractures away from the heel of wellbore contributes to decreasing the variance of fracture half-lengths, and enhancing the uniform growth of multiple hydraulic fractures.  To quantify the effect of non-uniform fracture spacings on the uniform growth of multiple fractures, the resulting variances for each case are compared as shown in Table 1. Little difference of the average half-length between each case is obtained. It can be shown that the variance of case 1 using the non-uniform spacings of 30 m and 20 m is slightly smaller than that of case 3 using the uniform spacing of 20 m, but is greatly larger than that of case 4 using the uniform spacing of 30 m. However, when using the non-uniform spacings of 20 m and 30 m, the variance drops greatly in comparison with case 3, but is still larger than that of case 4. It indicates that proper increase of spacing for fractures away from the heel of wellbore contributes to decreasing the variance of fracture half-lengths, and enhancing the uniform growth of multiple hydraulic fractures.

Non-Uniform Perforation Parameters
As presented in the model formulation, frictional pressure drop of perforation is caused when fracturing fluid flows into each cluster. Non-uniform perforation parameters for multiple hydraulic fractures propagation are simulated. Two main perforation parameters are the number and diameter. Increasing the number and diameter of perforation decreases the frictional pressure drop. To promote the uniform growth of multiple fractures, the number and diameter of perforation for the middle cluster is increased. Three cases are modeled, and the perforation parameters are presented in Table 2 as well as the above presented case with the uniform perforation parameters. The fracture propagation paths and widths distribution for these three cases are shown in Figure 14. It can be shown that the fractures are nearly equal in case 5 when using non-uniform perforation diameters and uniform diameter numbers. In case 6, only the diameter number is increased for the middle cluster, the fractures configuration is similar to that of case 8 with uniform perforation parameters. Further, we increase both of perforation number and diameter for the middle cluster, an interesting phenomenon is observed that the middle fracture is much longer than the left and right fractures, which is remarkably different from case 8 with uniform perforation parameters. The average half-length and variance for each case are compared as shown in Table 2. The variance of case 5 using non-uniform perforation diameters drops greatly in comparison with the case using uniform perforation parameters. Therefore, it is recommended that increasing the perforation diameter for the middle cluster to promote the uniform growth of multiple hydraulic fractures.

Single-Wing Fracturing
In above cases, two-wing fracturing is modeled. The single-wing fracturing is simulated in this section. Two-cluster fracture with 10 m spacing and three-cluster fracture with 20 m spacing are modeled under the condition of single-wing fracturing. The fracture propagation paths and widths distribution are shown in Figure 15. Compared with two-wing fracturing, the variance of half-length decreases from 33.5 to 0.06 for the twocluster fractures and decreases from 121.9 to 7.4 for the three-cluster fractures. It indicates that single-wing fracturing could greatly improve the uniform growth of multiple hydraulic fractures.
ters and uniform diameter numbers. In case 6, only the diameter number is increased for the middle cluster, the fractures configuration is similar to that of case 8 with uniform perforation parameters. Further, we increase both of perforation number and diameter for the middle cluster, an interesting phenomenon is observed that the middle fracture is much longer than the left and right fractures, which is remarkably different from case 8 with uniform perforation parameters. The average half-length and variance for each case are compared as shown in Table 2. The variance of case 5 using non-uniform perforation diameters drops greatly in comparison with the case using uniform perforation parameters. Therefore, it is recommended that increasing the perforation diameter for the middle cluster to promote the uniform growth of multiple hydraulic fractures.

Single-Wing Fracturing
In above cases, two-wing fracturing is modeled. The single-wing fracturing is simulated in this section. Two-cluster fracture with 10 m spacing and three-cluster fracture with 20 m spacing are modeled under the condition of single-wing fracturing. The fracture propagation paths and widths distribution are shown in Figure 15. Compared with twowing fracturing, the variance of half-length decreases from 33.5 to 0.06 for the two-cluster fractures and decreases from 121.9 to 7.4 for the three-cluster fractures. It indicates that single-wing fracturing could greatly improve the uniform growth of multiple hydraulic fractures.

Effect of Thermal Stress
To get an intuitive understanding, the thermal effect on single fracture propagation is analyzed firstly. And then two-cluster and three-cluster fracturing are simulated to discover the thermal effect on the simultaneous propagation of multiple fractures. The thermal parameters are set as follows: rock thermal conductivity λs = 20 W/(m·°C), rock volumetric thermal expansion coefficient α = 2.4 × 10 −5 1/°C, rock specific heat capacity cs = 800 J/(Kg·°C), rock density ρs = 2500 Kg/m 3 , fluid specific heat capacity cf = 4200 J/(Kg·°C), fluid density ρf = 1000 Kg/m 3 . The temperature of injected fracturing fluid is set as 30 °C.
For the single fracture propagation case, the evolution of fracture half-length with different initial reservoir temperatures is shown in Figure 16. Compared with the case without thermal effect, considering thermal effect results in shorter fracture half-length at the same injection time. Furthermore, fracture half-length decreases as the initial reservoir temperature becomes higher. With higher initial reservoir temperature, it would give rise to greater thermal stress on fracture surface, which relieves the situ-stress on fracture segments. Given the same injection rate, larger fracture width would be obtained, and then shorter fracture length is got. It indicates that the thermal effect decreases fracture propa-

Effect of Thermal Stress
To get an intuitive understanding, the thermal effect on single fracture propagation is analyzed firstly. And then two-cluster and three-cluster fracturing are simulated to discover the thermal effect on the simultaneous propagation of multiple fractures. The thermal parameters are set as follows: rock thermal conductivity λ s = 20 W/(m· • C), rock volumetric thermal expansion coefficient α = 2.4 × 10 −5 1/ • C, rock specific heat capacity c s = 800 J/(Kg· • C), rock density ρ s = 2500 Kg/m 3 , fluid specific heat capacity c f = 4200 J/(Kg· • C), fluid density ρ f = 1000 Kg/m 3 . The temperature of injected fracturing fluid is set as 30 • C.
For the single fracture propagation case, the evolution of fracture half-length with different initial reservoir temperatures is shown in Figure 16. Compared with the case without thermal effect, considering thermal effect results in shorter fracture half-length at the same injection time. Furthermore, fracture half-length decreases as the initial reservoir temperature becomes higher. With higher initial reservoir temperature, it would give rise to greater thermal stress on fracture surface, which relieves the situ-stress on fracture segments. Given the same injection rate, larger fracture width would be obtained, and then shorter fracture length is got. It indicates that the thermal effect decreases fracture propagating velocity. Regarding to two-cluster fracturing, similar phenomenon is observed in the tion of fracture half-length with different initial temperatures as shown in Figu when considering the thermal effect no matter what fracture spacing is chosen (1 20 m). Correspondingly, the fracture width increases as shown in Figure 17b whe sidering the thermal effect. It indicates that the thermal stress has a significant inf on fracture geometrical size. On the other hand, the fracture propagation paths are in Figure 18. Little difference can be seen between them, which indicates that the th stress has trivial effect on the fracture propagation path. The reason for this may b the thermal induced stress is mainly concentrated around fracture surface, and g decays away fracture surface, which would not interfere the intrinsic stress shad between fractures.  Regarding to two-cluster fracturing, similar phenomenon is observed in the evolution of fracture half-length with different initial temperatures as shown in Figure 17a when considering the thermal effect no matter what fracture spacing is chosen (10 m or 20 m). Correspondingly, the fracture width increases as shown in Figure 17b when considering the thermal effect. It indicates that the thermal stress has a significant influence on fracture geometrical size. On the other hand, the fracture propagation paths are drawn in Figure  18. Little difference can be seen between them, which indicates that the thermal stress has trivial effect on the fracture propagation path. The reason for this may be that the thermal induced stress is mainly concentrated around fracture surface, and greatly decays away fracture surface, which would not interfere the intrinsic stress shadowing between fractures. Regarding to two-cluster fracturing, similar phenomenon is observed in the evolution of fracture half-length with different initial temperatures as shown in Figure 17a when considering the thermal effect no matter what fracture spacing is chosen (10 m or 20 m). Correspondingly, the fracture width increases as shown in Figure 17b when considering the thermal effect. It indicates that the thermal stress has a significant influence on fracture geometrical size. On the other hand, the fracture propagation paths are drawn in Figure 18. Little difference can be seen between them, which indicates that the thermal stress has trivial effect on the fracture propagation path. The reason for this may be that the thermal induced stress is mainly concentrated around fracture surface, and greatly decays away fracture surface, which would not interfere the intrinsic stress shadowing between fractures.  As for three-cluster fracturing, the thermal effect on fracture half-length is als lar as shown in Figure 19a. In addition, the differences of half-length between the fi the second fractures are plotted in Figure 19b, which shows that the difference o length decreases greatly if the thermal effect is considered, especially when using fracture spacing. The average half-length and variance of each case are calcula shown in Table 3. Compared with the former case that neglects the thermal effe variance decreases from 121.9 to 22.7 and 8.8 with using the spacing of 20 m and de from 17.2 to 1.9 and 0.9 with using the spacing of 30 m if the thermal effect is cons It indicates that the thermal effect promotes the uniform growth of multiple fractu (a) (b) Figure 19. Evolution of fracture half-length and its difference with different initial reservoir atures. (a) Half-length of the first fracture, (b) Difference of half-length between the first second fractures.  As for three-cluster fracturing, the thermal effect on fracture half-length is also similar as shown in Figure 19a. In addition, the differences of half-length between the first and the second fractures are plotted in Figure 19b, which shows that the difference of half-length decreases greatly if the thermal effect is considered, especially when using small fracture spacing. The average half-length and variance of each case are calculated as shown in Table 3. Compared with the former case that neglects the thermal effect, the variance decreases from 121.9 to 22.7 and 8.8 with using the spacing of 20 m and decreases from 17.2 to 1.9 and 0.9 with using the spacing of 30 m if the thermal effect is considered. It indicates that the thermal effect promotes the uniform growth of multiple fractures. As for three-cluster fracturing, the thermal effect on fracture half-length is also similar as shown in Figure 19a. In addition, the differences of half-length between the first and the second fractures are plotted in Figure 19b, which shows that the difference of halflength decreases greatly if the thermal effect is considered, especially when using small fracture spacing. The average half-length and variance of each case are calculated as shown in Table 3. Compared with the former case that neglects the thermal effect, the variance decreases from 121.9 to 22.7 and 8.8 with using the spacing of 20 m and decreases from 17.2 to 1.9 and 0.9 with using the spacing of 30 m if the thermal effect is considered. It indicates that the thermal effect promotes the uniform growth of multiple fractures.  The comparison result of fracture propagation paths between the case that considers thermal effect and the one without thermal effect is shown in Figure 20. Although there is some deviation between fractures paths with spacing 20 m, the differences of fracture  The comparison result of fracture propagation paths between the case that considers thermal effect and the one without thermal effect is shown in Figure 20. Although there is some deviation between fractures paths with spacing 20 m, the differences of fracture propagation directions are rather small. It verifies that the thermal stress has little effect on fracture propagation paths again.

Conclusions
A coupled thermo-hydro-mechanical model is presented to simulate multiple hydraulic fractures propagation in a horizontal well. An iterative scheme is proposed to solve the coupling between fracture deformation, fluid flow and thermal stress. The proposed method is validated through comparison with analytical solution of two problems. Numerical results show that stress shadowing effect results in the non-uniform growth of multiple hydraulic fractures, especially when using small fracture spacings. However, using non-uniform fracture spacings and perforation parameters could improve this situation. Increasing the spacings for fractures away from the heel of wellbore and perforation diameter for the middle cluster promotes the uniform growth of fractures. Further, the thermal effect is of great significance to fracture geometrical size, but of trivial significance to fracture propagation paths. The thermal effect decreases fracture propagating velocity, but it is conducive to the uniform growth of multiple hydraulic fractures.

Conclusions
A coupled thermo-hydro-mechanical model is presented to simulate multiple hydraulic fractures propagation in a horizontal well. An iterative scheme is proposed to solve the coupling between fracture deformation, fluid flow and thermal stress. The proposed method is validated through comparison with analytical solution of two problems. Numerical results show that stress shadowing effect results in the non-uniform growth of multiple hydraulic fractures, especially when using small fracture spacings. However, using non-uniform fracture spacings and perforation parameters could improve this situation. Increasing the spacings for fractures away from the heel of wellbore and perforation diameter for the middle cluster promotes the uniform growth of fractures. Further, the thermal effect is of great significance to fracture geometrical size, but of trivial significance to fracture propagation paths. The thermal effect decreases fracture propagating velocity, but it is conducive to the uniform growth of multiple hydraulic fractures.