Analysis of Thermal Stimulation to Enhance Shale Gas Recovery through a Novel Conceptual Model

State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China Department of Oil-Gas Field Development Engineering, College of Petroleum Engineering, China University of Petroleum, Beijing 102249, China State Key Laboratory of Shale Oil and Gas Enrichment Mechanisms and Effective Development, SINOPEC Group, Beijing 050021, China ERE &BIC-ESAT, College of Engineering, Peking University, Beijing 100083, China School of Energy Resources, China University of Geosciences, Beijing 100083, China


Introduction
Shale gas is widely distributed with huge reserves and is an important component of the unconventional energy resources mix [1].With the growth of energy demands, the enhancement of shale gas recovery, i.e., the ratio of cumulative gas production and original gas in place, attracts increasing attentions worldwide [2].However, it is a challenging task to enhance the shale gas recovery due to the ultralow permeability, porosity, and high adsorption capacity of shale gas reservoirs [3].The determination of economic recoverable reserves of shale gas depends largely on how well the horizontal well drilling and hydraulic fracturing technologies have been applied.It has been proved that the multistage fracturing of horizontal wells is an effective technique to develop the shale gas reservoirs [4,5].However, the production rate of shale gas may still decline rapidly in a short period.Therefore, it is necessary to apply well stimulation to maintain stable production and economic recovery.Gas in shale reservoir exists in two forms: the adsorbed gas on the clay and kerogen and the free gas in the matrix pores and natural fractures.The free gas porosity and the porosity consumed by adsorbed gas around clay and organics are not independent but interconnected [6].Most of the free shale gas is located in large pores and can be relatively easily produced.It is the main component in the initial high production rate stage.In the long-term production, the adsorbed shale gas contributes as the primary source to the gas production.A large part of gas within shale (approximately 20% to 85%) exists as the adsorption phase [7,8]; thus, the production of adsorbed gas will have a profound impact on ultimate recovery.Many methods to enhance unconventional gas recovery have been investigated with the adsorbed gas as the target, such as the CO 2 injection enhanced gas recovery method [9] and thermal stimulation method [10].
In the gas injection methods, because of the stronger adsorption to CO 2 or N 2 compared to CH 4 , the injection of CO 2 or N 2 has been studied intensively to enhance coalbed methane and shale gas recovery, and the injection of CO 2 is an attractive strategy since it can also help reduce the level of CO 2 on the atmosphere through geological storage [11][12][13][14][15][16][17].
In the thermal stimulation methods, an important factor that may need special attention is the nonisothermal adsorption characteristics in the adsorption curve [7,18].The amount of adsorbed gas in shale can decrease with the increase of temperature when pressure is fixed, which provides the basis to use thermal stimulation technology in shale to enhance the shale gas recovery.Thermal recovery methods have been successfully applied to improve the heavy oil recovery but with the target to reduce the viscosity of heavy oil when heating formation [19].Many methods to raise the formation temperature have been developed in oil reservoir, such as commonly used steam injection [20] and steamassisted gravity drainage (SAGD) [21].An alternative method to heat the reservoir is the electromagnetic heating, which is a popular method to enhance oil recovery in recent years [22].As early as 1969, the alternating current heating method has been used to enhance oil recovery, and the heating mode relies on the frequency of electrical current [23].Sahni et al. [24] discussed the advantages and disadvantages of the low-frequency resistance heating and high-frequency electromagnetic heating.Their simulation results showed that utilizing a 60 kW microwave source, the formation temperature can be increased to a temperature of 422 K near the source within one year.With respect to the heating method in gas reservoir, Salmachi and Haghighi [18] confirmed that heat injection can increase the gas recovery in coal seam by more than 58% within 12 years, and the peak gas production rate was almost 6.8 times higher than conventional recovery methods.Recovery can also be greatly improved by injecting hot water, but the additional resistance brought by the blockage of gas flow channels when injecting hot water needs to be evaluated.Wu et al. [25] presented a model to simulate the injection of high temperature CO 2 into a coal seam to promote desorption of coalbed methane while simultaneously sequestering CO 2 in the coal seam.However, the coexistence of CH 4 and CO 2 gas flow makes the production forecast more difficult.Several literatures have considered the thermal impact on the cumulative production of shale gas by heating hydraulic fractures [26][27][28].Thermal stimulation technology has been applied to increase shale formation temperature around the hydraulic fractures, and shale gas recovery can be greatly improved due to the increased desorption rate at high temperature.However, two important factors that may cause estimation error of shale gas production have not been taken into account in these previous thermal stimulation studies.Firstly, only a single level of hydraulic fractures in the stimulated reservoir volume (SRV) region was considered into their models.Some researchers have confirmed that the hydraulic fractures in shale need to be characterized by fractal geometry [29][30][31].To accurately describe the fractal geometry of the hydraulic fracture distribution in shale reservoir, the simulation model should have the capability to deal with multiple levels of hydraulic fractures.Secondly, instantaneous local thermal equilibrium was assumed in the previously proposed models, which ignores the heat loss between fluids and rock.
A novel conceptual model is developed in this study to investigate the effect of thermal stimulation on the shale gas recovery in a complex hydraulic fracture network system.The flow mechanisms in shale gas reservoirs are unique due to its ultralow permeability.Knudsen diffusion, slippage effect, and non-Darcy flow have been incorporated into the shale gas flow model to accurately characterize the flow process.The nonisothermal temperature-dependent adsorption process is considered in the thermal stimulation process, and the local thermal nonequilibrium can explicitly describe the convective heat exchange between rock and fluids.Such adsorption model has not been coupled in shale gas flow model with complex fracture networks in the previous studies.The fluid flow equation is composed of gas flow equations through hydraulic fractures and matrix.The complexity of the hydraulic fractures has been accounted for by considering multiple levels of hydraulic fractures.Multiple levels of fractures are characterized by different fracture parameters, such as fracture permeability, fracture length, and fracture compressibility in the fracture flow equation.The thermal effect is brought into the fluid flow equation through adsorption term and gas viscosity term.The local thermal nonequilibrium between rock and fluid on the thermal-pressure coupling process has been also considered in heat transfer model.
The remainder of this paper is organized as follows: Section 2 provides the non-isothermal adsorption model, governing equations of matrix and hydraulic fractures, heat transfer equations of fluids and rock and the final form of the coupled solution.The model validation is described in Section 3, and it also presents a series of numerical experiments to investigate influence of thermal properties, fluids properties and reservoir parameters on thermal recovery from shale gas reservoirs.The conclusions are drawn in Section 4. The obtained results are helpful to understand the effects of thermal stimulation on ultimate shale gas recovery.

Mathematical Model
The schematic diagram of the reservoir simulation model is depicted in Figure 1.Due to the influence of hydraulic fracturing, the fracture network in the SRV region can be characterized by a fractal geometry with the tree shape [29,30].For the demonstrative purpose, the hydraulic fracture network is characterized by using two levels in this work, i.e., the primary fractures and the secondary fractures, respectively.More levels of hydraulic fractures can be easily accounted for by introducing additional fracture flow equations.The fractured network is composed of the primary fractures in the near-well zone and the secondary fractures extended 2 Geofluids from the tips of the primary fractures.The combination of hydraulic fractures and natural fractures creates a zone with high permeability around a horizontal well, called SRV region.
To differentiate the flow characteristics caused by the formation properties in and out of this region, the entire domain is divided into the SRV region and non-SRV (NSRV) regions as shown in Figure 1.For example, a high matrix permeability value in SRV region represents comprehensive effects of natural and hydraulic fractures.The permeability in the NSRV region maintains its original value without being affected by hydraulic fracturing.The reservoir temperature can affect gas adsorption capacity in shale.In general, the rate of gas desorption will increase when shale rock is heated, which provides a possible means to improve the shale gas recovery.A novel model to describe the influence of reservoir temperature on the gas production rate is proposed here.The proposed model is associated with a few assumptions: (1) it is assumed that the produced composition of shale gas is single-phase methane, and the multiple compositional effects on flow and sorption behaviors are not considered; (2) the homogeneous fluid and reservoir properties and uniform SRV region are used for the demonstrative purpose.The heterogeneity can be considered in the model, but it has been disregarded to improve the computational efficiency of the conceptual model; (3) the thermal properties, such as heat transfer coefficient and heat capacity, are characterized by constant and does not change during the heating process.Dynamic changing properties can be considered by using iterative process, but this nonlinearity may cause some numerical stability issue to the numerical solution.
2.1.Nonisothermal Adsorption Model.The Langmuir adsorption isotherm model describes the single molecular layer adsorption process under the state equilibrium condition [32].The absorbed gas content could be described by using Langmuir isotherm theory: where V L is Langmuir volume, m 3 /ton; b is adsorption coefficient, Pa −1 ; and p is pore pressure, Pa.However, during the process of thermal stimulation in shale gas reservoir, the adsorption characteristics cannot be accurately described by Langmuir isothermal adsorption model since the temperature may change dramatically in a short period of time.In this case, the nonisothermal adsorption model [26,33] must be used by modifying the adsorption content as a temperature-dependent variable: And the temperature-dependent adsorption coefficient can be expressed as where K 0 is a constant which is independent of temperature, K 1/2 /Pa; T is the temperature, K; E is the characteristic adsorption energy, J/mol; and R is the universal gas constant, 8.3145 J/mol/K. Figure 2 shows the influence of temperature on the adsorption capacity.The parameter values of K 0 , E, and V L are set to be 3.2288E − 9, −20,936 J/mol, and 200 scf/ton, respectively.It can be found that the adsorption capacity of shale gas decreases sharply with the increase of temperature.It implies that thermal stimulation can have a profound impact on the ultimate recovery, and the effect of temperature on the adsorption capacity becomes more significant when the formation pressure decreases.It can be also inferred that under the pressure of 10 MPa, the adsorption volume decreased dramatically from 140 scf/ton to 16 scf/ ton when the temperature increases from 350 K to 600 K.It is encouraging since it implies that the thermal stimulation method may have a remarkable potential to enhance the shale gas recovery due to the large portion of the adsorbed shale gas.

Continuity Equation for Fluid Flow in Shale Matrix.
Different from conventional gas reservoirs, gas flow in shale reservoir can be affected by Knudsen diffusion, slippage effect, and non-Darcy flow.The Darcy's law can be modified as where μ g is viscosity, Pa ⋅ s; ∇ = ∂x, ∂y is the gradient operator.
The Knudsen diffusion, slippage effect in the matrix can be taken into account by the apparent permeability: where k ga is the apparent permeability, m 2 ; k g is the intrinsic gas permeability, m 2 ; D g is the diffusion constant, m/s; c g is

Horizontal well
Primary fractures Secondary fractures SRV region the compressibility, Pa −1 ; K n is Knudsen number; and F D is the permeability coefficient.
The temperature-dependent diffusivity equation for shale gas transport in the matrix of SRV region can be expressed as The temperature-dependent diffusivity equation for shale gas transport in the matrix of NSRV region can be expressed as where k gai is apparent permeability in the SRV region, m 2 ; k gao is apparent permeability in the NSRV region, m 2 ; ρ g is the gas density; ρ sc is the density of gas under the standard conditions, kg/m 3 ; ρ s is the density of shale rock, kg/m 3 ; and t is the time variable, s.According to (2) and ( 3), the derivatives of the adsorption capacity V with respect to pressure can be written as And the derivatives of the adsorption capacity V with respect to temperature can be written as where T is the temperature, K and p is the matrix pressure, Pa.
The gas density ρ g can be calculated by using equation of state for real gas: where M is the average molecular weight of mixed gas, kg/mol; Z factor [34] can be estimated through the corresponding pseudoreduced pressure p pr and temperature T pr : The compressibility of gas can be determined through critical pressure p c : And the compressibility of gas at the critical pressure c gp can be determined by [35] The pressure-temperature-dependent gas viscosity can be obtained by [36] μ g = ae bρ c g 14 The coefficients a, b, and c [36] can be compute as

Continuity Equation for Flow in Fracture.
A large set of fractures in shale reservoir can be produced by hydraulic fracturing, and the created fractures can be divided into primary fractures and secondary fractures as shown in Figure 1.The generation of hydraulic fractures can improve shale permeability and thus increase the production of shale gas.By using Darcy's law, velocity of gas flow in the primary and secondary fractures can be, respectively, described as 4 Geofluids where k f and k fs are the permeability of primary and secondary fractures, respectively, m 2 and p f and p fs are formation pressure in primary fractures and secondary fractures, respectively, Pa.
According to (18) and ( 19), the continuity equations of primary fractures and secondary fractures can be, respectively, written as where c f is primary fracture compressibility, Pa −1 and c fs is the secondary fracture compressibility, Pa −1 .

Heat Transfer Equations.
To simulate the influence of thermal recovery on the shale gas production, it is required to understand the temperature distribution in the shale gas reservoir.Thus, the heat transfer equation needs to be established.There are three main methods in the literatures to build the heat conduction equation: equivalent temperature method which treats the temperature in the porous medium as a single value, matrix-fracture temperature method which assumes that thermodynamic equilibrium between rock and fluids at all times, and rock-fluid temperature method which treats temperature individually between rock and fluids [37,38].The first two methods are local isothermal methods.The heat conduction between rock and fluids under local nonisothermal condition is considered in the proposed model in this paper.By using Fourier's law, the heat transfer equations for fluid and rock can be described as where T is fluid temperature, K; T r is rock temperature, K; λ f is thermal conductivity of fluids, W/m/K; λ r is thermal conductivity of rock, W/m/K; h a is heat transfer coefficient between fluids and rock, W/m 3 /K; ϕ is porosity; c pf is heat capacity of fluids, J/kg/K; c pr is heat capacity of rock, J/kg/ K; and v g is fluid velocity vector in (4), m/s.It is noted that the fluid flow and temperature are coupled together through (22) and (23), where the gas flow velocity is incorporated into the heat transfer equation.

Boundary Conditions.
To obtain the final solution, proper boundary conditions need to be set.Here, the gas flow rate and formation pressure at the boundary of rock matrix and hydraulic fractures should satisfy the following conditions: Flow rate and pressure conditions at the junction of the SRV and NSRV zones can be given as follows: where the vectors v g , v m , and v s represent velocity in the matrix, the primary fracture, and the secondary fracture, respectively.The vector n represents outward directional vector normal to the boundary.The superscripts m and s represent the primary and secondary fractures, respectively.Γ m and Γ s represents the boundaries of matrix and primary fractures and these of matrix and secondary fractures.Γ SRV and Γ NSRV represent the outer boundaries of SRV region and inner boundaries of NSRV region.A finite rectangular region is used in the simulation of this paper.All outer boundaries are closed and no external heat source affects the temperature distribution in the simulation domain.Therefore, the outer boundary conditions for flow rate and temperature should meet the following conditions: Inner boundary conditions for hydraulic fracture wall are assumed to be constant production pressure and constant temperature heating at fracture walls, which could be measured and known in some specific cases.
The matrix flows (6) and ( 7), fracture flow (20) and ( 21), and heat transfer (22) and ( 23) together with boundary conditions (24), ( 25), ( 26), ( 27), ( 28), ( 29), (30), (31), and (32) constitute the system of coupled equations.Due to the complexity of the above equations, it is difficult to obtain an explicit analytical solution.In this work, the system of coupled equations is solved by using the finite element method (FEM).Triangular meshes are used in the twodimensional model.The fractures are assumed to be fully penetrated; thus, this two-dimensional model can be representative of the real situation.The coarsening of the grid is applied throughout the matrix reservoir.However, mesh refinement must be used around the hydraulic fractures in 5 Geofluids order to obtain accurate and reliable results as shown in Figure 3.A flow chart of the overall numerical procedure is presented in Figure 4.

Results and Discussion
The accuracy of the coupled model is validated by comparing the numerical solution with an available analytical solution in the literature [39].After the validation of the solution accuracy, a series of numerical experiments are performed to investigate the influences of thermal properties, fluid properties, and reservoir parameters on the performance of thermal-stimulated gas recovery in shale gas reservoirs.

Model Validation.
To validate the proposed model, a published analytical solution in [39] to deal with the similar problem has been adopted as a benchmarking model to compare with the proposed model.In the benchmarking model, an analytical solution is derived to simulate the pressure transient, as shown in equation (46) in the reference, and production behaviors of the fractured horizontal wells in unconventional shale gas reservoir.The temperature of thermal stimulation is set as initial temperature, and the reservoir pressure before the implementation of thermal stimulation is set as initial pressure.The production rate and cumulative production are used as model output so that these two models can be compared in a consistent manner.However, no previous model can be as theoretically comprehensive as the proposed one here.The benchmarking model uses isothermal adsorption and simple fracture geometry without secondary fractures; therefore, the proposed model is simplified to be consistent with the benchmarking model for comparative purpose.The settings of the reservoir and fracture properties are as shown in Table 1. Figure 5 shows the validation results by comparing the numerical results based on the proposed model with the analytical solution.It can be found that the production rate and cumulative gas production obtained from the proposed model agree with these obtained from the analytical solution very well.The slightly higher values obtained from the analytical solution relative to the proposed model, in terms of both cumulative production and production rate, are caused by trilinear flow hypothesis considered in the analytical solution.2. The distribution of temperature during the process of thermal stimulation is depicted in Figure 6.The temperature on the walls of both primary and secondary fractures is set to 520 K, which is higher than the initial temperature of shale gas reservoir so that it can represent the effect of the injected thermal fluid.The pressure of inner boundary, i.e., the walls of fractures, is set to 3 MPa, which is lower than initial pressure of shale gas reservoir in order to represent the pressure drop induced by the gas recovery process.The numerical simulation results in Figure 6 show that as the time increases, the distance of the temperature diffusion increases gradually.At a given time, the temperature diffusion caused by the thermal effect weakens in the area far from the primary and secondary fractures.
Figure 7 shows the temperature profile averaged in the transverse direction around the horizontal well on a crosssectional area from the minimum x coordinate 0 m to the maximum x coordinate 1200 m in the heated reservoirs (T in = 520 K).As shown in Figure 7, there exist three peaks in the temperature profile, which represent the temperature values on the intersections of the horizontal well and hydraulic fractures due to the constant temperature heating.The temperature between two adjacent stages of hydraulic fractures tends to decrease as it moves away from the hydraulic fractures.As shown in Figure 6, at a given location within SRV region, the temperature increases with time.6 Geofluids After 20 years, the temperature change caused by thermal stimulation spreads to the boundary of the SRV region.Figure 8 shows pressure distribution in the shale formation during the thermal stimulation process with the temperature of the thermal stimulation as T in = 520 K.The pressure of both primary fractures and secondary fractures is set to be 3 MPa, which is lower than the initial pressure of shale gas reservoir.The temperature of fracture wall is set to be the same as the thermal stimulation.The results in Figure 8 show that as the time increases, the rate of the pressure propagation increases rapidly.The fastest pressure drop rate occurs in the vicinity area of both primary and secondary fractures.The formation pressure drop spreads to the boundary of SRV region after just 1 year and tends to propagate further at a slower rate into the NSRV region.This is attributed to the fact that the hydraulic fracturing greatly improves the permeability of the rock matrix in the SRV region due to the reopening of the natural fractures, in addition to create large scale hydraulic fractures.

Horizontal well Fractures
Figure 9 shows the pressure profile with the temperature of thermal stimulation as T in = 520 K, and the pressure profile is obtained along the same cross section as the temperature distribution in Figure 7.The time of thermal stimulation is set to be 1 year, 5 years, 10 years, and 20 years, respectively.As shown in Figure 9, there exist three local minima in the graph, which represent the pressure distribution on the intersections of horizontal well and fractures due to constant pressure on the fracture walls.The pressure remains high in the SRV region between two adjacent stages of the hydraulic fractures.With the proceeding of the gas production, the peak pressure values between hydraulic fractures decrease.The sharp pressure changes near the boundaries of the SRV region (i.e., x = 150 m and x = 1050 m) can be attributed to the difference of the formation properties in the SRV and NSRV regions.Compared with temperature spreading, the rate of pressure propagation is much faster.After 20 years, the pressure on the outer boundary is reduced from initial pressure 20 MPa to 14 MPa.
The pressure and temperature distributions during the shale gas production process have been investigated independently in the previous section.To analyze the effect of the  7 Geofluids thermal stimulation on the formation pressure distribution, a series of test cases with different thermal stimulation temperature have been studied.The data used here is the same as those listed in Table 2 except that the thermal stimulation temperatures are set to be 393.15K (the same as the initial reservoir temperature without thermal stimulation), 520 K, and 620 K. Figure 10 shows the effect of thermal stimulation temperature on the pressure profile at different production time, e.g., 1 year, 5 years, 10 years, and 20 years, respectively.The pressure profile is obtained along the direction of horizontal well as before.It can be observed that the pressure drop becomes less severe in general when the temperature of the thermal stimulation is elevated.This is caused by the fact that the thermal stimulation can enhance the desorption effect of shale gas, which can increase the formation pressure as the desorbed shale gas enter the pore of shale matrix.The pressure increased caused by shale gas desorption can compromise the pressure drop during the production process in the hydraulic fractured reservoir.It indicates that the shale gas production can greatly benefit from the thermal stimulation method.One can also find that the thermal stimulation mainly affects the pressure distribution in the SRV region in the early stage of the thermal-stimulated shale gas production, and the thermal effect on pressure change in the NSRV region becomes observable after 10 of shale gas production as shown in Figure 10(c).This is attributed to the temperature diffusion process similar to what has been shown in Figure 6, where the thermal stimulation mainly affects the area in the vicinity of hydraulic fracture in the early stage and it takes a long time for the temperature to diffuse to the boundary of SRV region before it can have an effective influence on the pressure distribution in the NSRV region.Therefore, it is necessary to account for the coupling effect of both temperature and pressure.

Effect of Thermal Stimulation on Shale Gas Recovery.
Shale gas recovery largely depends on thermal properties, such as heat transfer coefficient, gas thermal conductivity, rock thermal conductivity, and heat capacity of rock.8 Geofluids walls) are set to be 3 MPa and 520 K in this analysis.The other parameter values are the same as those listed in Table 2.
Figure 11 shows the effect of volumetric heat transfer coefficient h a (the product of convective heat transfer coefficient and surface area) on shale gas recovery, which is used to investigate the effect of local thermal nonequilibrium on gas recovery and pressure field.Volumetric heat transfer coefficient h a represents rock-gas exchange in the heating reservoir and a large h a value can reflect strong heat exchange between rock and gas.As shown in Figure 11, h a values are set to 0 W/m 3 /K, 0.0001 W/m 3 /K, 0.01 W/m 3 /K, and 1 W/ m 3 /K, respectively.When h a is equal to 0 W/m 3 /K, there is no heat exchange between the rock and gas.It can be found from Figure 11 that the gas recovery increases as the production progresses and it reduces with the increase of h a value at a given time, which implies that the gas recovery will be overestimated if the heat exchange between rock and gas is ignored.The recovery rate of the shale gas decreases with the increase of the volumetric heat transfer coefficient.The reason is that a larger volumetric heat transfer coefficient leads to more heat transferred into the gas phase during the heating process and the lowered rock temperature results in less shale gas desorption.There exists a limit of the effect of volumetric heat transfer coefficient on gas recovery.In these selected values of volumetric heat transfer coefficient, gas recovery will stay the same when h a value is higher than 0.01 W/m 3 /K.This threshold value may change with the thermal property parameters, such as thermal conductivity, but the recovery tendency does not change under the different thermal properties.
Figures 12 and 13 show the effects of thermal conductivity of gas, λ f , and thermal conductivity of rock, λ r , on the shale gas recovery.The thermal conductivity of gas, λ f , is set to be 0.01 W/m/K, 0.05 W/m/K, and 0.10 W/m/K, respectively.The thermal conductivity of rock, λ r , is set to be 1 W/ m/K, 10 W/m/K, and 20 W/m/K, respectively.The thermal conductivity value represents the capability of heat transmission.The thermal conductivity is related to the structure, density, humidity, temperature, pressure, and other factors.Generally, the thermal conductivity of rock increases with the increase of pressure, density, and humidity and decreases with the increase of temperature.It can be found in Figures 12 and 13 that the gas recovery increases with the increase of λ f and λ r values, which indicate that gas recovery will be enhanced at a higher thermal conductivity of gas or rock.Figure 14 shows the shale gas recovery when heat capacity of rock is set to be 500 J/kg/K, 1000 J/kg/K, and 1500 J/kg/ K, respectively.Heat capacity represents the required energy when the temperature of 1 kilogram of rock rises by 1 Kelvin, which implies that a small heat capacity value will result in a rapid increase of temperature.As shown in Figure 14, the gas recovery will decrease with the increase of heat capacity of rock and the reduction of gas recovery will become insignificant when heat capacity of rock C pr value is higher than 1000 J/kg/K in this case.

Effect of Thermal Stimulation Temperature and
Production Pressure on Shale Gas Recovery.Thermal stimulation temperature is an important control factor to enhance shale gas recovery in the process of thermal stimulation.Figure 15 shows the effect of thermal stimulation temperature on shale gas recovery.Thermal stimulation temperature along primary fractures and secondary fractures is set to be 393.15K with no heating source, 520 K, and 620 K, respectively.It can be observed that shale gas recovery can be enhanced due to the increase of the thermal stimulation temperature.Desorption rate increases with the increase of thermal stimulation temperature and gas recovery increases by 10% after 25 years as the result of increase of thermal stimulation temperature from 393.15 K with no heating source to 520 K.However, the increase rate of gas recovery decreases when the thermal stimulation temperature increases from 520 K to 620 K.The reason is that the temperature diffusion has reached its maximum coverage area and thus desorption capability has reached its upper limit.Therefore, thermal stimulation is an effective method to enhance shale gas recovery but the thermal stimulation temperature has to be reasonably selected in the application of these techniques.
Production pressure is another important control factor to enhance shale gas recovery.Figure 16 shows the effect of production pressure on shale gas recovery.Production pressures are set to be 3 MPa and 10 MPa, and thermal stimulation temperatures associated with each pressure value are set to be 393.15K with no heating source and 520 K with  heating sources.Simulations results show that shale gas recovery decreases with increase of production pressure when the temperature are the same.Here, the gas recovery is reduced by 22% and 18% for nonheating case and heating case, respectively.Both depressurization and heating up are helpful to enhance shale gas recovery.Thermal stimulation is more effective with high temperature in terms of the improvement level of gas recovery.This implies that the thermal stimulation temperature will play a dominant role in the development of shale gas since a certain level of production pressure has to be maintained to ensure the reasonable production of shale gas during the process of thermal stimulation.

Conclusions
In this study, a novel conceptual model is developed to couple the nonisothermal adsorption with the shale gas flow model in a complex hydraulic fracture network system.To investigate the effect of thermal stimulation on the shale gas recovery, it is necessary to use the nonisothermal model, which can describe the desorption process influenced by temperature and pressure.Due to the ultralow permeability, gas flow in shale has to account for Knudsen diffusion, slippage effect, and non-Darcy flow.All factors are coupled into the mathematical model, including matrix flow, primary and secondary fracture flow, and heat transfer in shale reservoir.Nonisothermal heat transfer is also considered in the proposed model.The coupled model is solved by using finite element method.The proposed model can be validated through the comparison with an analytical proposed in the literature.Pressure distribution and temperature distribution are simulated during the process of thermal recovery.The temperature between two adjacent stages of hydraulic fractures tends to decrease as it moves away from the hydraulic 11 Geofluids fractures.After 20 years in the test case, the temperature change caused by thermal stimulation spreads to the boundary of the SRV region.The thermal stimulation mainly affects the pressure distribution in the SRV region in the early stage of the thermal-stimulated shale gas production, and the thermal effect on pressure change in the NSRV region becomes observable after 10 years of shale gas production in the test case.
Thermal properties and reservoir properties can significantly affect the shale gas recovery.These thermal properties can be characterized by volumetric heat transfer coefficient, gas thermal conductivity, rock thermal conductivity, heat capacity of rock, thermal stimulation temperature, and production pressure.Gas recovery decreases with the increase of volumetric heat transfer coefficient, and there exists a limit of the effect of volumetric heat transfer coefficient on gas recovery.In this case, the threshold value of volumetric heat transfer coefficient is 0.01 W/m 3 /K.Gas recovery increases with the increase of gas and rock thermal conductivity.It decreases with the increase of heat capacity of rock and the reduction of gas recovery will become insignificant when heat capacity of rock is higher than 1000 J/kg/K in this case.
Thermal stimulation is an effective method to enhance shale gas recovery, but the thermal stimulation temperature has to be reasonably selected in the application of these techniques since the increase rate of gas recovery decreases when the thermal stimulation temperature increases from 520 K to 620 K. Shale gas recovery decreases with increase of production pressure when the temperature remains the same.Here, 12 Geofluids the gas recovery is reduced by 22% and 18% for nonheating case and heating case.The thermal stimulation temperature is a more important factor in the development of shale gas, since a reasonably high production pressure is required to ensure the successful production of shale gas.Through the proposed model, the thermal stimulation process in shale gas reservoir can be more accurately described, which can be beneficial to predict and optimize the enhancement of shale gas recovery with the thermal stimulation.The average fluid properties, average reservoir properties, and uniform SRV are used in this model, the effects of the heterogeneous fluid and reservoir properties and the dynamic changing of thermal property parameters on the thermal stimulation process will need to be studied in the future.

Figure 1 :
Figure 1: Schematic diagram of fracture networks in shale reservoir after hydraulic fracturing.

Figure 2 :
Figure 2: Adsorption capacity change of shale gas with temperature under different pressure.

3. 2 .
Temperature and Pressure Distribution during the Process of Thermal Recovery.The model parameters used in the numerical tests are listed in Table

Figure 3 :
Figure 3: Unstructured triangular mesh grid used in finite element numerical computation with the refinement around hydraulic fractures and horizontal well.

Figure 4 :
Figure 4: Workflow of the proposed thermal-stimulated gas recovery analysis model.

Figure 5 :
Figure 5: Validation of the proposed model through the comparison against analytical solution.

Figure 6 :Figure 7 :
Figure 6: Distribution of temperature during the process of thermal stimulation (a) at 1 year, (b) at 5 years, (c) at 10 years, and (d) at 20 years.

Figure 8 :
Figure 8: Distribution of formation pressure during the process of thermal stimulation (a) at 1 year, (b) at 5 years, (c) at 10 years, and (d) at 20 years.

Figure 9 :
Figure 9: Pressure distribution profile in the thermal-stimulated shale reservoir.

Figure 10 :
Figure 10: Effect of thermal stimulation temperature on formation pressure during the process of thermal stimulation (a) at 1 year, (b) at 5 years, (c) at 10 years, and (d) at 20 years.

Figure 15 :
Figure 15: Effect of thermal stimulation temperature on shale gas recovery.

Figure 16 :
Figure 16: Effect of production pressure on thermal recovery.

Table 1 :
Setting of the model parameters in the accuracy validation case.

Table 2 :
Setting of model parameters in the numerical tests.