The Hanyi Wang

A Numerical Study of Thermal-Hydraulic-Mechanical (THM) Simulation with the Application of Thermal Recovery in Fractured Shale Gas Reservoirs*

Keywords: Shale Gas, Thermal Recovery, Thermal-Hydraulic-Mechanical (THM) Modeling, Numerical Simulation, Fractured Reservoirs, Discrete Fracture Network (DFN)

Shale gas is playing an important role in transforming global energy markets with increasing demands for cleaner energy in the future. One major difference in shale gas reservoirs is that a considerable amount of gas is adsorbed. Up to 85% of the total gas within shale may be found adsorbed on clay and kerogen. How much of the adsorbed gas can be produced has a significant impact on ultimate recovery. Even with improving fracturing and horizontal well technologies, average gas recovery factors in U.S. shale plays is only ~30% with primary depletion. Adsorbed gas can be desorbed by lowering pressure and raising temperature, reservoir flow capacity can be also influenced by temperature, so there is a big prize to be claimed using thermal stimulation techniques to enhance recovery. To date, not much work has been done on thermal stimulation of gas shale reservoirs.

In this article, we present general formulations to simulate gas production in fractured shale gas reservoirs for the first time, with fully coupled thermal-hydraulic-mechanical (THM) properties. The unified shale gas reservoir model developed in this study enable us to investigate multi-physics phenomena in shale gas formations. Thermal stimulation of fractured gas reservoirs by heating propped fractures is proposed and investigated. This study provides some fundamental insight into real gas flow in nano-pore space and gas adsorption/desorption behavior in fractured gas shales under various in-situ conditions and sets a foundation for future research efforts in the area of enhanced recovery of shale gas reservoirs.

We find that thermal stimulation of shale gas reservoirs has the potential to enhance recovery significantly by enhancing the overall flow capacity and releasing adsorbed gas that cannot be recovered by depletion, but the process may be hampered by the low rate of purely conductive heat transfer, if only the surfaces of hydraulic fractures are heated.

Introduction

Unconventional resources-shale gas and liquids, coalbed methane, tight gas and heavy oil-are called "unconventionals" because in order to economically access and produce these hydrocarbons, unconventional methods and expertise are required. Better reservoir knowledge and increasingly sophisticated technologies (especially horizontal drilling and hydraulic fracturing) make the production of unconventional resources economically viable and more efficient. This efficiency is bringing shale reservoirs, tight gas and oil, and coalbed methane into the reach of more companies around the world. With ever increasing demanding for cleaner energy, unconventional gas reservoirs are expected to play a vital role in satisfying the global needs for gas in the future. The major component of unconventional gas reservoirs comprises of shale gas. Shales and silts are the most abundant sedimentary rocks in the earth's crust and it is evident from the recent year's activities in shale gas plays that in the future shale gas will constitute the largest component in gas production globally. According to GSGI, there are more than 688 shales worldwide in 142 basins and 48 major shale basins are located in 32 countries [START_REF] Newell | Shale Gas and the Outlook for U.S. Natural Gas Markets and Global Gas Resources[END_REF]. Shale gas exploitation is no longer an uneconomic venture with the availability of improved technology, as the demand and preference for cleaner form of hydrocarbon are in ever greater demands, especially in country like China, where the main energy resource still comes from coal. Unlike conventional gas reservoirs, shale gas reservoirs have very low permeability and are economical only when hydraulically fractured. The key techniques that allow extracting shale gas commercially such as horizontal drilling and hydraulic fracturing, are expected to improve with time; however as better stimulation techniques are becoming attainable, it is important to have a better understanding of shale gas reservoir behavior in order to apply these techniques in an efficient fashion. One important aspect of shale gas reservoirs which needs special consideration is the adsorption/desorption phenomenon.

In organic porous media, gas can be stored as compressed fluid inside the pores or it can be adsorbed by the solid matrix. Similar to surface tension, adsorption is a consequence of surface energy [START_REF] Gregg | Adsorption, surface area and porosity[END_REF], which causes gas molecules to get bonded to the surface of the rock grains. The gas adsorption in the shale-gas system is primarily controlled by the presence of organic matter and the gas adsorption capacity, which depends on TOC (Total Organic Carbon), organic matter type, thermal maturity and clay minerals [START_REF] Ambrose | New pore-scale considerations for shale gas in place calculations[END_REF][START_REF] Passey | From Oil-Prone Source Rock to Gas-Producing Shale Reservoir-Geologic and Petrophysical Characterization of Unconventional Shale Gas Reservoirs[END_REF]. Generally, the higher the TOC content, the greater the gas adsorption capacity. In addition, a large number of nanopores lead to significant nanoporosity in shale formations, which increases the gas adsorption surface area substantially. The amount of adsorbed gas varies from 35-58% (Barnett Shale, USA) up to 60-85% (Lewis Shale, USA) of total gas initial in-place [START_REF] Darishchev | On Simulation of Flow in Tight and Shale Gas Reservoirs[END_REF]. Presently, the only method for accurately determining the adsorbed gas in a formation is through core sampling and analysis. However, understanding the effects that initial adsorption, and moreover, desorption has on gas production will increase the effectiveness of reservoir management in these challenging environments.

Besides horizontal drilling and hydraulic fracturing, which ignited the momentum of "shale revolution", thermal stimulation methods have been widely used for unconventional reservoirs, such as heavy oil and shale oil. During the last two decades, the development of thermal stimulation technologies, such as in-situ combustion, cyclic steam injection and SAGD (steam assisted gravity drainage), have played a major role in the implementation of different concepts of oil production from unconventional oil reserves. Today, about 60% of world oil production attributed to methods of enhanced oil recovery (EOR) comes from thermal stimulation (Chekhonin et al., 2012). As an alternative to traditional thermal methods to reduce oil viscosity and enhance recovery, electric/electromagnetic heating method has been proposed and tested in heavy oil and oil shale reservoirs [START_REF] Sahni | Electromagnetic heating methods for heavy oil reservoirs[END_REF][START_REF] Hascakir | Investigation of productivity of heavy oil carbonates and oil shales using electrical heating methods[END_REF], to convert oil shale to producible oil and gas through heating the oil shale in situ by hydraulically fracturing the oil shale and filling the fracture with electrically conductive material to form a heating element [START_REF] Symington | ExxonMobil's Electrofrac Process for In Situ Oil Shale Conversion[END_REF], as shown in Fig. 1. Already, microwave applications in oil sands bitumen and shale oil production and in petroleum upgrading are gaining considerable interest in recent years. Energy companies and petroleum researchers have been working on a variety of unconventional technologies such as microwave and radio frequency (RF) energies to recover viscous oil from shale [START_REF] Mutyala | Microwave applications to oil sands and petroleum: A review[END_REF], and it is observed that microwaves can heat up the formation much faster than conventional steam heating. In addition, nanoparticles in the form of nanofluids have been investigated for enhanced oil recovery applications and it has shown that due to absorption of electromagnetic waves by the cobalt ferrite nanoparticles, oil viscosity can be reduced, resulting in a significant increase of oil recovery [START_REF] Yahya | Cobalt ferrite nanoparticles: An innovative approach for enhanced oil recovery application[END_REF]. Even though numerous studies have investigated how to improve heavy oil/shale oil recovery by increasing formation temperature, limited studies have explored the possibility of enhancing gas recovery with similar thermal stimulation techniques. Salmachi et al. (2012) investigated the feasibility of enhancing gas recovery from coal seam gas reservoirs using geothermal resources. The results show that hot water injection for a period of 2 years into the coal seam with an area of 40 acres successfully increases average reservoir temperature by 30 ℃ and dramatically increases gas recovery by 58% during 12 years of production. Wang et al. (2015a) also concluded that thermal stimulation has the potential to enhance CBM recovery substantially by liberating a significant amount of residual adsorption gas. [START_REF] Wang | Conceptual Study of Thermal Stimulation in Shale Gas Formations[END_REF] investigated the application of thermal stimulation in hydraulically fractured shale gas formations by altering gas adsorption/desorption behavior in multiple transverse hydraulic fractures. Their work shows the efficiency of thermal treatment in shale gas formations largely depends on fracture spacing, operation conditions, gas adsorption and rock properties. [START_REF] Chapiro | Thermal Well Stimulation in Gas Shales Through Oxygen Injection and Combustion[END_REF] investigated the possibility of in-situ combustion to improve permeability in shale gas formations. They concluded that if kerogen present in sufficient quantities, methane combustion can generate enough heat to enhance the permeability. [START_REF] Yue | Increasing Shale Gas Recovery through Thermal Stimulation: Analysis and an Experimental Study[END_REF] conducted laboratory experiment on shale samples, by measuring shale gas adsorption capacity with different temperatures, their work indicates that large amount of adsorption gas can be expelled by elevating temperature and the sensitivity of gas adsorption capacity to temperature depends on rock mineralogy properties.

The potential application of thermal stimulation in shale gas formations adds extra complexity in reservoir simulation and modeling, because the rising temperature can substantially impact real gas properties, the release of adsorption gas, pressure depletion process and the associated matrix and fracture flow capacity. So it is crucial to have a reservoir simulation model that can deal with these multi-physics problems with enough confidence, however, no reservoir simulation model, that is able to encompass all the coupled physics, is currently available and no literature has attempted to do so.

In this study, we propose a fully coupled thermal-hydraulic-mechanical (THM) model for fractured shale gas reservoirs and evaluate the feasibility of using thermal stimulation in hydraulically fractured shale formations to enhance the ultimate recovery. The results of this study provide us a better understanding of the coupled behavior of fluid flow, rock deformation and heat transfer under complex reservoir conditions, and help us assess the effectiveness and potential applications of thermal stimulation methods in fractured shale gas formations. The structure of this article is as follows. First, a brief description of shale gas reservoir THM modeling will be presented. Then, this fully coupled model is applied to different thermal stimulation scenarios and the effects the thermal stimulation are assessed. Finally, conclusion remarks and discussions are presented. Detailed mathematic formulations for the proposed model are discussed in Appendix.

Thermal-Hydraulic-Mechanical Modeling in Fractured Shale Gas Reservoirs

Be able to model and simulate reservoir depletion and well performance properly in fractured formations is crucial for shale gas field development. Geological modeling and reservoir simulation provide essential information that can help reduce risk, enhance field economics, and ultimately maximize gas reserves by identifying the number of wells required, the optimal completion design and the appropriate enhanced recovery methods. The reliability of reservoir simulation results is primarily dependent on the quality of input parameters and the physical modeling of the reservoir-production systems. And the reservoir simulation model itself should be robust, representative and be able to make the most use of available data that come from various sources. Besides extremely low formation matrix permeability, some other unique features of shale gas reservoir, which can substantially impact the reservoir simulation results, should not be ignored.

In classic fluid flow mechanics where continuum theory holds, fluid velocity is assumed to be zero at the pore wall [START_REF] Sherman | The Transition From Continuum to Molecular Flow[END_REF]. This is a valid assumption for conventional reservoirs having pore radii in the range of 1 to 100 micrometer, because fluids flow as a continuous medium. Correspondingly, Darcy's equation, which models pressure-driven viscous flow, works properly for such reservoirs. However, in shale reservoirs, the ultrafine pore structure of these rocks can cause violation of the basic assumptions behind Darcy's law. Depending on a combination of pressure-temperature conditions, pore structure and gas properties, Non-Darcy flow mechanisms such as Knudsen diffusion and Gas-Slippage effects will impact the matrix apparent permeability [START_REF] Fathi | Shale Gas Correction to Klinkenberg Slip Theory[END_REF][START_REF] Michel | Parametric Investigation of Shale Gas Production Considering Nano-Scale Pore Size Distribution, Formation Factor, and Non-Darcy Flow Mechanisms[END_REF][START_REF] Swami | Non Darcy Flow in Shale Nanopores: Do We Have a Final Answer[END_REF][START_REF] Sakhaee-Pour | Gas Permeability of Shale[END_REF]. In addition, constant decreasing pore pressure during production (transient flow and pseudo-steady-state flow) can lead to reduction of thickness in gas adsorption layer and increase in the effective stress, which in turn, can impact the formation matrix microstructure and effective pore radius [START_REF] Wang | Impact of Shale Gas Apparent Permeability on Production: Combined Effects of Non-Darcy Flow/Gas-Slippage, Desorption and Geomechanics. SPE-173196-PA[END_REF]. So the overall matrix apparent permeability in shale gas reservoir is dynamic and pressure dependent, as shown in Fig. 2. In the context of hydraulic fractured shale gas formations, the local matrix permeability in the Stimulated Reservoir Volume (SRV) is space and time dependent during production. [START_REF] Wang | Impact of Shale Gas Apparent Permeability on Production: Combined Effects of Non-Darcy Flow/Gas-Slippage, Desorption and Geomechanics. SPE-173196-PA[END_REF] Besides the matrix permeability, the surface area of natural fracture networks that connected to the main hydraulic fracture and its ability to sustain conductivity are also critical for predicting long-term production in shale gas formations [START_REF] Ghassemi | Sustaining Fracture Area and Conductivity of Gas Shale Reservoirs for Enhancing Long-Term Production and Recovery[END_REF]. It is recommended that Brinell Hardness Test (BHN) and Unpropped Fracture Conductivity Test (UFCT) should be done in shales in order to determine fracture treatment types and estimate the relationship between fracture conductivity and confining stress [START_REF] Ramurthy | Surface-Area vs. Conductivity-Type Fracture Treatments in Shale Reservoirs[END_REF]. So the effects of pressure dependent matrix permeability and fracture conductivity should be included in the simulation model for any well performance and production prediction.

Fig. 2-Mechanisms that alter shale matrix apparent permeability during production

The inclusion of temperature effects add more complexity in reservoir modeling and simulation, because the changes in formation temperature not only alter gas adsorption capacity, vary real gas properties, induce thermal stress, but also impact matrix permeability and fracture conductivity in a fully coupled manner, which in turn, affect in-situ flow capacity and ultimate recovery. Fig. 3 shows laboratory measurement and theoretical model prediction of gas adsorption capacity at different temperatures from a shale sample [START_REF] Yue | Increasing Shale Gas Recovery through Thermal Stimulation: Analysis and an Experimental Study[END_REF]. By examining the adsorption curves, we can deduce that the reservoir pressure must be sufficiently low to liberate the adsorbed gas and the ultimately recoverable amount of gas is largely a function of the adsorbed gas that can be released (desorbed). Because most adsorbed gas can only be released at low reservoir pressure, due to the extremely low permeability in shale matrix, even with hydraulic fracturing, it would take considerable production time for the average pressure within the drainage area to drop to a level where most of the adsorbed gas can be liberated, and the production rate may have already reached the economical shut-in limits by then. However, if the reservoir temperature can be elevated, a significant amount of adsorption gas can be released, even at high pressure (golden line in Fig. 3). Thus, thermal stimulation techniques can be utilized as a potential method to enhance ultimate recovery from shale gas reservoir by altering shale gas desorption behavior. With the elevated formation temperature, thermal stimulation induced rock mechanical behavior can alter rock properties and hence impact the fluid flow. Numerous studies [START_REF] Burghignoli | A laboratory study on the thermomechanical behaviour of clayey soils[END_REF][START_REF] Cekerevac | Experimental Study of Thermal Effects on the Mechanical Behaviour of a Clay[END_REF][START_REF] Laloui | Thermo-Plasticity of Clays: An Isotropic Yield Mechanism[END_REF][START_REF] Wong | Structural Integrity of Casing and Cement Annulus in a Thermal Well Under Steam Stimulation[END_REF][START_REF] Xu | Thermal Impact on Shale Deformation/Failure Behaviours -Laboratory Studies[END_REF][START_REF] Yuan | Design of Caprock Integrity in Thermal Stimulation of Shallow Oil-Sands Reservoirs[END_REF] have investigated the behavior of different rocks during the process of thermal stimulation, indicating that the rocks can be weakened or strengthened depending on the factors such as initial porosity, applied confining pressure, heating rate and rock composition. However, the relevance of these results to field enhanced recovery project using thermal stimulation remains poorly understood. So in this article, we do not consider dynamic rock damage mechanisms during thermal stimulation, but effects of induced thermal stress on associated, coupled physics are well captured. A comprehensive discussion on THM modeling and mathematical formulations are presented in Appendix.

Model Description

In this section, the application of the proposed THM model is presented and discussed. Even though it is possible to simulate an entire section of a horizontal well with multiple transverse fractures, it is more efficient to simulate a unit pattern and apply symmetric boundary conditions along boundary of SRV that contains one main hydraulic fracture and natural fractures, as shown in Fig. 4 (the red line depicts horizontal wellbore, the blue line, and the distributed black line segments represent hydraulic fracture and natural fractures respectively).

Fig.4-A plane view of SRV containing hydraulic fracture and natural fractures

In practice, the primary hydraulic fracture and secondary fracture networks within each SRV unit, can be created statistically by using seismic, well and core data [START_REF] Elfeel | Effect of DFN Upscaling on History Matching and Prediction of Naturally Fractured Reservoirs[END_REF][START_REF] Cornette | Refining Discrete Fracture Networks With Surface Microseismic Mechanism Inversion and Mechanism-Driven Event Location[END_REF] and by matching hydraulic fracture propagation models. Hydraulic fracture models that commonly used in industry relies on linear elastic fracture mechanics (LEFM) (Johri and Zoback 2013;[START_REF] Weng | Modeling of Complex Hydraulic Fractures in Naturally Fractured Formation[END_REF]), but it is only valid for brittle rocks [START_REF] Wang | SPE-168600-PA. Poroelastic and Poroplastic Modeling of Hydraulic Fracturing in Brittle and Ductile Formations[END_REF], recent studies demonstrate that using energy based cohesive zone model can resolve the issue of singularity at the fracture tip and these cohesive zone based hydraulic fracture models not only can be used in both brittle and ductile formations, but also can capture complex fracture evolution, such as natural fracture interactions [START_REF] Guo | Numerical simulation of interaction of hydraulic fracture and natural fracture based on the cohesive zone finite element method[END_REF], fracture reorientation from different perforation angles [START_REF] Wang | Numerical Modeling of Non-Planar Hydraulic Fracture Propagation in Brittle and Ductile Rocks using XFEM with Cohesive Zone Method[END_REF], producing well interference (Wang 2016a), the effects of fracturing spacing and sequencing on fracture interaction/coalescence from single and multiple horizontal wells (Wang 2016b).

A limitation of the widely used commercial reservoir simulators is the fracture usually modeled as a narrow space with high porosity/permeability. A grid block that contains fracture segment is given higher porosity/permeability. The smallest grid block is often much thicker than typical fracture width. This requires the fracture property to be homogenized within the containing grid blocks. The orthogonal grid based mesh, used by many reservoir simulators, also imposes some limitation on fracture shape and geometries. In order to avoid extreme mesh refinement inside and around the fractures and guaranteeing solution convergence, discrete fracture networks (DFN) and unstructured mesh are implemented in our proposed model to delineate the fluid flow behavior along the fractures and discretize reservoir domain. Fig. 5 shows the unstructured mesh generated in the simulation domain that discretizes both reservoir and discrete fracture networks. The maximum and minimum element size is 0.5 m and 0.1 m, respectively and there are 96574 elements in total.

Fig.5-Unstructured mesh and DFN distribution within simulation domain

Combine hydraulic fracture modeling and microseismic data, the spatial distribution of fracture networks can be estimated by various methods (Johri and Zoback, 2013;[START_REF] Yu | Discrete Fracture Network Generation from Microseismic Data using Moment-Tensor Constrained Hough Transforms[END_REF][START_REF] Aimene | Geomechanical modeling of hydraulic fractures interacting with natural fractures -Validation with microseismic and tracer data from the Marcellus and Eagle Ford[END_REF]. And the properties of the propped fracture and un-propped fracture can be determined through conductivity test under various pressures and confining stresses [START_REF] Ghassemi | Sustaining Fracture Area and Conductivity of Gas Shale Reservoirs for Enhancing Long-Term Production and Recovery[END_REF], so that the parameters define transport capacity of fractures can be determined (e.g., fitting Eq.A-20 with laboratory data). With all these information, the special and intrinsic properties of DFNs can be estimated. How to map fracture networks across different fracturing stages based on microseismic data/hydraulic fracture modeling and fully characterize the entire horizontal completion is beyond the scope of this article. The distribution of natural fractures is randomly assigned within SRV in this study as pictured in Fig. 5.

Table 1 shows all the input parameters used for the synthetic case analysis, which includes reservoir conditions, drainage geometry, fracture conductivity, real gas parameters and rock thermal properties for a typical shale formation. The diameter of adsorption gas molecules is assumed to be 0.414 nm, which is the diameter of a methane molecule. The combination of input parameters is designed to reflect a scenario where 60% of initial gas in-place comes from adsorbed gas, this can be the case in many organic-rich shale formations [START_REF] Darishchev | On Simulation of Flow in Tight and Shale Gas Reservoirs[END_REF]. [START_REF] Wang | Conceptual Study of Thermal Stimulation in Shale Gas Formations[END_REF]. Fig. 6 shows the prediction of gas adsorption capacity at stimulation temperature (478 K) by fitting original Langmuir curve (determined by 𝑉 𝐿 and 𝑃 𝐿 ) at initial reservoir temperature (366 K) using Bi-Langmuir model.

Input

Fig.6-Shale gas adsorption capacity with different pressure and temperature

Base Case Simulation

In this section, we present our simulation results for the base case, in which the formation is under isothermal condition throughout the lifetime of reservoir depletion. In the following section, the results from base case simulation will be compared with the results when thermal stimulation is applied. Fig. 7 shows the pressure distribution in the simulated SRV unit after 1 and 20 years of production. It can be observed that the pressure disturbance has reached the boundary after 1 year's production, owing to the existence of conductive fracture networks, and pseudo-steady-state flow has already ensued from the initial transient flow, to dominate the rest of the production history. After 20 years of production, most of the area that penetrated by the main hydraulic fracture has been completely depleted due to the well-connected fracture system and effective linear flow from matrix to fracture, but the area in front the main hydraulic fracture tip has not been depleted effectively. Fig. 8 and Fig. 9 depicts the gas density and gas viscosity distribution in the simulated SRV unit after 1 and 20 years of production. It is within our expectation that under isothermal reservoir conditions, both gas density and gas viscosity are correlated to reservoir in-situ local pressure. The expansion of gas volume with decreasing gas density during depletion is the main driven mechanism of production, and the gas viscosity decreases as pore pressure declines. We can also observe that the gas viscosity is reduced locally to around 50% of its initial value in low-pressure zone, so if constant gas viscosity is used in reservoir simulation or production prediction models, the inaccurate results may lead to erroneous decisions. The corresponding matrix apparent permeability distribution is shown in Fig. 11. It can be observed that the value of matrix apparent permeability is higher within the depleted zone. And as the pressure sink propagating away from the near-hydraulicfracture-region to the entire SRV, the matrix apparent permeability increases with declining pressure, due to the effects of Non-Darcy flow/Gas-Slippage and the release of adsorption layer in the nano-pore space during production. Fig. 12 shows how much adsorption gas has been produced within the SRV unit, which can be a good indicator to estimate the ultimate recovery. We can notice that after 20 years of production, nearly 58% of adsorbed gas has been produced in the depleted zone, while only close to 30% of adsorbed gas has been produced in the area in front of the main hydraulic fracture tip. From the simulation results, we can also conclude that the maximum percentage adsorption gas that can be recovered through pressure depletion is around 58%, under our simulated production conditions and gas adsorption/desorption properties that provided in Table 1.

Fig.12-Percentage of adsorption gas produced in the SRV after 1 and 20 years of production

In order to demonstrate the importance of incorporating the evolution of both fracture conductivity and matrix permeability during production, two additional scenarios are constructed with same input parameters as the base case, but assuming constant matrix permeability and constant fracture conductivity, respectively. Because fracture conductivity is the product of fracture permeability and fracture width, so by assuming the fracture width remains the same, shift fracture permeability as a function of local stress (Eq. A.20) has the same effect as shift fracture conductivity on simulation results. In addition, the fracture storage effects have negligible impact on gas production, because the driving forces of primary shale gas recovery come from the gas volume expansion and gas desorption with declining pressure. Fig. 13 shows the predicted cumulative production for three different scenarios. It can be clearly seen that overlook the pressure dependent fracture permeability/conductivity can lead to overestimation of production, while ignoring the pressure dependent matrix permeability can result in underestimation of production. In fact, the relative importance of matrix apparent permeability evolution no only depends on reservoir micro-structure, pressure-temperature conditions, but also influenced by the density, connectivity and conductivity of natural fractures. If the overall flow capacity inside the SRV is dominated by abundant fracture networks, then the general role of matrix permeability itself diminishes (Wang 2016c). Being able to predict production correctly is crucial to optimize field development and reservoir management, so a comprehensive reservoir characterization is required to be able to delineate reservoir performance with enough certainty, because the flawed physical model can render unreliable prediction and poor history match even with available production data.

Fig.13-Prediction of cumulative production with different simulation models

Case 1: Thermal Stimulation through Hydraulic Fracture

In this section, we simulated the case of directly heating the propped hydraulic fractures, assess the proposed method of using electromagnetically sensitive proppants, and investigate how the increased formation temperature can impact reservoir depletion, cumulative production, and ultimate recovery. In this case, the thermal stimulation temperature is set to be 475K along the main hydraulic fracture and remains constant throughout the production life; all the other input parameters maintain the same as the base case. It can be expected that with elevated formation temperature, more adsorption gas can be released at the same pressure, as illustrated in Fig. 6: the extra gas can be produced is equivalent to the difference of shale rock adsorption capacity between two temperature states at a specific local pressure, if do not consider the changes in formation flow capacity.

In addition, all the pressure-temperature dependent variable will be affected by the rising temperature, so does the fully coupled THM interactions.

Fig. 14 shows the temperature distribution within the simulated SRV after 1 and 20 years of production. It can be observed that the front edge of stimulation temperature has not reached the boundary yet after 1 year's production and the temperature distribution profile reflects a heat transfer process that dominated by the conduction of rock matrix. The temperature propagation speed depends on formation thermal properties and thermal stimulation temperature. Less time will be needed to heat up the formation temperature to the desired value if the formation exhibits higher thermal conductivity or higher stimulation temperature is imposed. After 20 years of production under current simulation conditions, most of the SRV that penetrated by the hydraulic fracture has reached the desired formation temperature. Compare with the pressure depletion profile in Fig. 7, we can infer that the efficiency of pressure depletion and heat transfer is substantially impacted by the length of hydraulic fracture that penetrated into the drainage volume, even with the existence of conductive natural fractures. So unlike conventional reservoirs, it is the length and spacing of hydraulic fracture that determines the actual drainage area in shale gas formations.

results from the base case, the impact of elevated formation temperature on gas transport and desorption mechanisms can be clearly observed: The average pressure in the SRV (Fig. 15-a) is higher than that in the base case (Fig. 7) at the end of 20 years production when thermal stimulation is implemented, especially within the regime in front of the hydraulic fracture tip. This can be explained by the fact that more adsorption gas is released due to increasing formation temperature, which in turn compensates parts of the pressure depletion and leads to a lower rate of reservoir pressure decline. Compare Fig. 15-b with Fig. 10, we can note that the Knudsen Number increases from 0.5451 to 0.8151 in the SRV regimes where the hydraulic fracture penetrates. This is because the temperature and pressure dictate the alteration of Knudsen Number during production and the elevated temperature contributes to higher Knudsen Number (as reflected by Eq.A.21). However, within the regime in front of the hydraulic fracture tip, the lowest Knudsen Number even drops from 0.229 (base case) to 0.1909 (with thermal stimulation). This contrast effects of thermal stimulation in different regimes is due to the fact that, unlike the regime penetrated by the hydraulic fracture, the regime far ahead of the hydraulic fracture tip has not been impacted temperature propagation yet (Fig. 14), and the pressure in this regime is higher because pressure depletion has been compensated by the release of extra adsorption gas in the regime where temperature is elevated. This opposite effects of thermal stimulation in these two divided regions also applies to the matrix apparent permeability and the percentage of adsorption gas produced. Higher Knudsen Number in the hydraulic fracture penetrated SRV indicates more severe the Non-Darcy flow/Gas-Slippage behavior. In addition, elevated temperature also lead to higher intrinsic permeability resulted from thinner adsorption layer, that's why the maximum matrix apparent permeability increases from 330-nd (base case) to 461-nd (with thermal stimulation) in the simulated SRV, as we can discern from Fig. 15-c Fig. 16 shows the predicted cumulative production for the base case and the thermal stimulated case. The cumulative production after 20 years increases 40% when thermal stimulation is applied along the hydraulic fracture, with 122K difference between the thermal stimulation temperature and initial reservoir temperature. With elevated formation temperature, more adsorption gas can be released ultimately under current well production conditions. In addition, higher temperature leads to higher matrix apparent permeability due to stronger Non-Darcy flow/Gas Slippage behavior and slower pressure decline trend, which in turn, can relax fracture permeability reduction. On top of that, more gas can be expelled from the reservoir because of gas volume expansion at a higher temperature. All these factors lead to higher production and enhanced ultimate recovery.

Fig.16-Prediction of cumulative production with and without thermal stimulation, k ∞0 =100 nd

Next, we investigate how shale matrix permeability can influence thermal stimulation efficiency. Low formation permeability adversely impact the pressure diffusion process and in such cases, the transient flow in SRV may dominate the entire production life and the average pressure in the SRV will not be able to decline to a level that allows sufficient desorption before production rates reach economical limits. Here we reduce the initial intrinsic matrix permeability from 100 nd to 10 nd and maintain other input parameters as the same as provided in Table 1. Fig. 17 and Fig. 18 depict the formation pressure and mean effective stress in the SRV after 1 and 12 months of production, respectively. We can observe that after 1 month of production, only the regimes that around natural fractures that well-connected to the main hydraulic fracture are depleted, this is due to the low permeability of the formation and the fracture networks serve as the main path for gas flow. Also, the asymmetric nature of fracture network distribution leads to the disproportion of depletion with the SRV on both sides of the hydraulic fracture after 1 year of production. With constant far field stress, the local mean effective stress evolves with pressure depletion. 19 shows the predicted cumulative production for the base case and the thermal stimulated case, with the same initial matrix intrinsic permeability of 10 nd. The cumulative production after 20 years increases 50% when thermal stimulation is applied along the hydraulic fracture. Compared to Fig. 16, where around 40% of extra gas can be produced with thermal aid, it seems that thermal stimulation is more efficient when formation permeability is lower. This can be explained by the fact that the lower the formation permeability, the less effective of pressure depletion, so temperature elevation has more impact on gas production. Even though, thermal stimulation has the potential to enhance gas recovery significantly, but the time requires elevating formation temperature to the target level is substantial, because of low efficiency of heat transfer by pure conduction (as shown in Fig. 14), so the effect of thermal stimulation on ultimate recovery can be only prominent in late times. Considering a typical shale gas production decline trend, shown in Fig. 20, where shale gas wells lose 80% of their initial production rate in the first 2 years, and the economically preferred completion designs may be more driven by the net present value derived in the first a few years of production rather than the ultimate recovery of the well. So for practical applications of thermal stimulation techniques in shale gas reservoirs, the heat transfer efficiency has to be improved. The process of proppant transport and settling in complex fracture networks is a complicated phenomenon, where the final distribution of proppant depends on a series of factors, such as fracture geometry, pump rate, proppant concentration, proppant size and injection fluid rheology. In reality, only a fraction of the fracture networks in the SRV is propped and the effective propped volume (EPV) may play a dominant role in determining long-term well performance. Fig. 21 shows the reservoir drainage prediction in SRV and EPV after 20 years of production (the fracture network is shown in black, the virgin pressure is red, SRV in blue rectangular and EPV in green rectangular). It can be noted that most drainage happens inside EPV. Despite ongoing efforts to enhance our understanding of proppant transport in complex fracture geometry [START_REF] Raymond | Estimation of the Propped Volume Through the Geomechanical Modeling of Multiple Hydraulic Fractures Interacting with Natural Fractures[END_REF][START_REF] Roy | Proppant Transport at the Fracture Scale: Simulation and Experiment[END_REF][START_REF] Sahai | Laboratory Results of Proppant Transport in Complex Fracture Systems[END_REF], how to correctly quantifying proppant distribution within the fracture network on a field scale is still poorly understood and further investigations are required. In this section, we investigate a scenario where the heat source particles can be pushed into fracture networks, so more formation area can be heated up simultaneously. To this end, we assume that all the natural fractures that well connected to the hydraulic fracture can be filled with heat source particles, which is the upper-limit case on how large the fracture surface area can be thermally stimulated at the same time. The intention of this section is not to mimic proppant distribution in a realistic field case, but to assess the impact of heat transfer efficiency on gas production and the effectiveness of thermal stimulation.

Fig. 22 shows the temperature distributions in the SRV under such conditions. Compared to Fig. 14, where only the hydraulic fracture acts as heat source, the heat transfer efficiency increases dramatically. When most of the natural fracture networks can be thermally stimulated collectively, the whole SRV regime that penetrated by the main hydraulic fracture can achieve the desired temperature level around 1 year. With such high efficiency of heat transfer, the overall thermal stimulation effects on gas production can be altered significantly, as reflected in Fig. 23. When compared to Fig. 16, where only the main hydraulic fracture is thermal stimulated, the effect of thermal stimulation becomes prominent (cumulative production almost doubled) even within the first 5 years of production. However, as time goes on, the effects of thermal stimulation gradually tapers off as most the formation have already reached the target temperature within the first a few years. This result demonstrate that if the efficiency of heat transfer process can be improved to desorb large amount of adsorption gas and enhance the overall flow capacity during the initial stage of gas producing well, then the ultimate recovery during the economic lifetime of a gas producing well (typically, first 5 years of production matters most) can be substantially enhanced.

Fig.23-Prediction of cumulative production with and without thermal stimulation, k ∞0 =100 nd

Energy Balance Calculations

Different thermal stimulation methods have different capital cost, that related to the maturity of associated technology, the supply of raw material, the efficiency of management and the environmental impact, etc. And the price of natural gas, like other commodities, has the nature of volatile and cyclicity. So, to present a comprehensive evaluation from the economic aspects of thermal stimulation project is out of the scope of this article. Nevertheless, we can still assess the feasibility of thermal stimulation in shale formations by investigating the energy needed for thermal stimulation and the corresponding extra energy output.

The heat content of natural gas, or the amount of energy released when a volume of gas is burned, varies according to gas compositions. The primary constituent of natural gas is methane, which has a heat content of 1,010 British thermal units per cubic foot (Btu/cf), equivalent to 3.763 × 10 7 J/m 3 , at standard temperature and pressure.

Fig. 24 shows the heat content of natural gas consumed in various locations in the United States, and the average heat content of natural gas in the United States was around 1,030 Btu/cf (3.838 × 10 7 J/m 3 ). In the following calculations, we use the heat content of pure methane as a reference value to determine the volume of natural gas needed to generate an equivalent amount of input energy for thermal stimulation. The energy supplied during thermal stimulation can be estimated using the following equation:

Energy Input = ∭ ρ m C m (T -T 0 ) V ( 1
)
where ρ m and C m are formation rock density and heat capacity, respectively. T 0 is the reservoir initial temperature and T the is rock current local temperature. Integrate ρ m C m (T -T 0 ) over the entire simulated volume, the total energy input can be determined.

Based on simulation results from previous sections and energy input and output calculations, the net extra natural gas produced via thermal stimulation can be estimated, as reported in hydraulic fracture is heated, we can observe that the net recovery of extra gas increases as time goes on, regardless of matrix permeability. In absolute terms, more net extra gas can be produced when permeability is higher, but in percentage terms (compared to the scenario where thermal stimulation is not applied), thermal stimulation seems more efficient when permeability is lower. One can also notice that the calculated energy inputs in case 1 are the same, even for different matrix permeability, this is because only conductive heat transfer is modeled, due to the negligible convective heat transfer effects by gas ( refer to Eq.A.11-A.14). The consequences can be quite different if all the well-connected fracture networks can be heated up simultaneous, as reflected in case 2. The results indicate that much more net extra gas can be produced when heat transfer efficiency is improved, even within 5 years of thermal stimulation. It also shows that longer stimulation time leads to more net recovery of extra gas in absolute terms, but result in less production of extra gas by percentage. This implies that if we assess the thermal stimulation efficiency using net extra gas produced by percentage as an indicator, there exists an optimal duration on how long thermal stimulation should be applied. From the above results and analysis, we can conclude that thermal stimulation has great potential to improve shale gas ultimate recovery significantly by allowing the release of residual adsorption gas that cannot be recovered by pressure depletion alone, and increasing the overall flow capacity in the formation, by increasing matrix apparent permeability and relaxing/delaying the reduction of fracture conductivity. However, the effectiveness of thermal stimulation on ultimate recovery largely hinges on how much time it takes to elevate the formation temperature to a target level. Because the elevation of formation temperature has compound effects on all the coupled underlying physics, and the impact of these effects are cumulative. The sooner the formation can be heated up to the desired temperature, the longer time the reservoir can benefit from its altered properties, and hence, the effectiveness of thermal stimulation can be enhanced significantly.

Conclusions and Discussions

In this study, a thermal stimulation method that combines transverse hydraulic fracture and electromagnetically sensitive proppants is proposed, and the possibility of enhancing shale gas recovery behavior by elevating formation temperature is explored. In order to assess the overall effects of thermal stimulation on gas production, we presented general formulations of numerically modeling and simulation of gas production in fractured shale gas reservoirs for the first time, with fully coupled thermal-hydraulic-mechanical mechanisms. The process of real gas flow in shale matrix/ fractures, rock deformation due to changing in-situ stress and heat transfer by conduction are coupled together through temperature-pressure dependent variables.

The results indicate that the gas production and ultimate recovery have the potential to be improved substantially by elevating temperature, which alters the gas adsorption/desorption behavior and the overall flow capacity in fractured shales. But the heat transfer efficiency (i.e., how long it will take to heat up the drainage volume to desired temperature) has dominant impact on the effectiveness of thermal stimulation.

It should be mentioned that this paper only focuses on the effects thermal stimulation on gas adsorption/desorption behavior and the evolution of matrix permeability and fracture conductivity during production. It would be more realistic if the shearinduced slippage, rock failure due to thermal stress, and all the other possible mechanisms (such as the impact of temperature on clay bound water, capillary bound water and the transformation of organic matter into hydrocarbons) can be incorporated into the model to assess the overall thermal stimulation effects. In addition, the applications of other thermal stimulation methods, such as microwave, nano-particle fracturing fluid, etc., need to be explored to improve the efficiency of heat transfer process. Because different shale plays have different rock, mineral, and organic matter compositions, their gas adsorption capacity have various sensitivity to temperature, so laboratory experiment is crucial in understanding to what degree the temperature can impact gas adsorption/desorption process and how to optimize thermal stimulation design correspondingly.

Considering an average shale gas recovery factor of 30% in U.S. shale plays, there is a big prize to be claimed in terms of enhanced recovery using thermal stimulation techniques in shale gas reservoirs, however, not much work has been done in this area yet. The general approaches and the proposed mathematical framework of thermal-hydraulic-mechanical (THM) modeling presented in this article set a foundation for future research efforts, to better understand fully-coupled physics and explore new enhanced recovery techniques in fractured shale gas reservoirs.

𝑘 ∞ 𝑘 ∞0 = ( 𝑟 𝑟 0 ) 𝜂 (A. 29)
where 𝜂 is the coefficient that define the sensitivity of permeability to the changes of pore radius. Different shale formations may have different nano-pore structure typology, which leads to different values of 𝛽. In this study, 𝜂 equals 2 by assuming the overall all intrinsic permeability resembles fluid flow in a capillary tube (Beskok and Karniadakis 1999). Combining Eq. (A. 28), Eq. (A. 29) and Eq. (A. 22), we have the final expression of matrix apparent permeability:

𝑘 𝑎 = 𝑘 ∞0 (𝑟 0 ( 𝜎 𝑚 𝜎 𝑚0 ) -0.5𝐶 ∅ -𝑑 𝑚 𝑃/𝑃 𝐿 1 + 𝑃/𝑃 𝐿 ) 𝜂 𝑟 0 𝜂 𝑓(𝐾 𝑛 ) (A. 30)
Compared to the apparent matrix permeability model proposed by [START_REF] Wang | Impact of Shale Gas Apparent Permeability on Production: Combined Effects of Non-Darcy Flow/Gas-Slippage, Desorption and Geomechanics. SPE-173196-PA[END_REF], Eq. (𝐀. 𝟑𝟎) is a more general formula for shale matrix apparent permeability, with additional parameter 𝛈 to relate the changes in pore radius to the alterations in intrinsic permeability. 𝛈 can be determined from laboratory experiment for different types of shales.

As the pore pressure declines and temperature increases, formation porosity, pore radius, Knudsen number and intrinsic permeability, are subject to time and space variation in SRV, and so does the matrix apparent permeability. This matrix apparent permeability model enables us to assess the evolution of matrix flow capacity under coupled effects of pressure depletion and heat transfer process during thermal stimulation. With all these pressure and temperature dependent variables substituting into the general governing equations described in the previous sections, a closed system of fully coupled thermalhydraulic-mechanic shale gas reservoir simulation model is constructed and ready to be solved by available numerical methods. Based on the formulas introduced above, a unified shale gas reservoir simulator was constructed. Newton-Raphson method [START_REF] Ypma | Historical development of the Newton-Raphson method[END_REF] and finite element analysis [START_REF] Zienkiewicz | The Finite Element Method[END_REF] are used to solve all the coupled equations numerically. = Matrix porosity at reference conditions

Nomenclature

Fig. 1 -

 1 Fig. 1-Heating hydraulic fracture with electrically conductive material (Hoda et al., 2010)

Fig. 3 -

 3 Fig. 3-Laboratory measurement and Bi-Langmuir model prediction of gas adsorption capacity at different temperature for a shale sample (Yue et al., 2015)

Fig. 7 -

 7 Fig.7-Pressure distribution in the SRV after 1 and 20 years of production

Fig. 8 -Fig. 9 -

 89 Fig.8-Gas density distribution in the SRV after 1 and 20 years of production

Fig. 10 -

 10 Fig.10-Knudsen Number in the SRV after 1 and 20 years of production

Fig. 11 -

 11 Fig.11-Shale matrix apparent permeability in the SRV after 1 and 20 years of production

Fig. 14 -

 14 Fig.14-Temperature distribution in the SRV after 1 and 20 years of production Fig. 15 shows the distribution of pressure, Knudsen Number, matrix apparent permeability and percentage of adsorption gas produced in the simulated SRV after 20 years of production, with the application thermal stimulation. Compare with previous

  and Fig. 11. More interestingly and importantly, when comparing Fig. 15-d with Fig. 12, it can be observed that the maximum recovery of the adsorption gas in the SRV increased from 58.38% to 81.33% with the help of temperature elevation. The distribution of pressure (a), Knudsen Number (b), matrix apparent permeability (c) and percentage of adsorption gas produced (d) in the SRV after 20 years of production with thermal stimulation

Fig. 17 -

 17 Fig.17-Pressure distribution in the SRV after 1 and 12 months of production

Fig. 19 -

 19 Fig.19-Prediction of cumulative production with and without thermal stimulation, k ∞0 =10 nd

Fig. 20 -

 20 Fig.20-Typical shale gas production decline trend under different scenarios (Wang, 2016c) Case 2: Thermal Stimulation through Fracture Networks

Fig. 21 -

 21 Fig.21-Simulation of reservoir drainage after 20 years of production (Maxwell, 2013)

Fig. 22 -

 22 Fig.22-Temperature distribution in SRV when well-connected fracture networks are thermal stimulated simultaneously

Fig. 24 -

 24 Fig.24-Heat content of natural gas consumed in various locations in the United States (U.S Energy Information Administration, 2014)

Table 1 . Input parameters for simulation cases

 1 Bi-Langmuir model (Lu et al. 1995) is used to depict pressure-temperature dependent gas adsorption capacity. By regression on common gas adsorption data (i.e., Langmuir volume constant 𝑉 𝐿 and Langmuir pressure constant 𝑃 𝐿 ), gas adsorption capacity at different temperatures can be extrapolated

	Parameters	Value
	Initial reservoir pressure, 𝑝 0	38[MPa]
	Wellbore pressure, 𝑝 𝑤𝑓	5 [MPa]
	Initial intrinsic matrix permeability, 𝑘 ∞0	100[nD]
	Initial formation porosity, ∅ 𝑚0	0.01
	Initial matrix pore radius (not include adsorption layer), 𝑟 0	3[nm]
	Porosity compaction parameter, 𝐶 ∅	0.035
	Initial hydraulic fracture permeability, 𝑘 𝑓0	10000[mD]
	Hydraulic fracture porosity, ∅ 𝑓	0.5
	Hydraulic fracture width, 𝑑 𝑓	0.01[m]
	Fracture compaction parameter for hydraulic fracture, B f	0.005[1/ MPa]
	Initial natural fracture permeability, 𝑘 𝑛𝑓0	400[mD]
	Natural fracture porosity, ∅ 𝑛𝑓	0.01
	Natural fracture width, 𝑑 𝑛𝑓	0.0001[m]
	Fracture compaction parameter for hydraulic fracture, B nf	0.05[1/ MPa]
	Initial reservoir temperature, 𝑇 0	366[K]
	Thermal stimulation temperature, 𝑇 𝑠	478[K]
	Formation heat capacity, 𝐶 𝑚	1,000[J/K/kg]
	Formation heat conductivity, λ 𝑚 Density of formation rock, 𝜌 𝑚 Langmuir volume constant, 𝑉 𝐿	4[W/m/K] 2600 [kg /𝑚 3 ] 0.0113[m 3 /kg]
	Langmuir pressure constant, 𝑃 𝐿	10 [MPa]
	Diameter of adsorption gas molecules, 𝑑 𝑚	0.414[nm]

Table 2 .

 2 From the results of case 1, where only the surface of

		0.4	Thermal Stimulated thorugh Fracture Networks			Base Case
	m 3 )	0.35
	(Billion	0.3
	Production	0.2 0.25
	Cumulative	0.05 0.1 0.15
		0
		0	2	4	6	8	10	12	14	16	18	20
						Production Time (Year)

Table 2 . Energy balance calculation at the 5 th and 15 th benchmark years

 2

	Case 1	Case 2

  𝑎 𝑖 = Pre-exponential constant for ith type of adsorption, L𝑡 2 / m, 1/𝑃𝑎 𝐵 = Parameter for fracture pressure dependent permeability, L𝑡 2 / m, 1/𝑃𝑎 𝐶 𝑚 = Heat capacity of rock matrix, 𝐿 2 /𝑡 2 /𝑇, J/K/kg 𝐶 𝑝,𝑔 = Heat capacity of gas phase at constant pressure, 𝐿 2 /𝑡 2 /𝑇, J/K/kg 𝐶 ∅ = Material constant for pressure dependent porosity 𝑑 𝑓 = Fracture width, L, 𝑚 𝑑 𝑚 = Diameter of absorbed gas molecules, L, 𝑚 𝐸 = Young's modulus, m/L𝑡 2 , 𝑃𝑎 𝑓(𝐾 𝑛 ) = Non-Darcy flow correction term 𝑓 𝑖 = Fraction of ith type of adsorption 𝐹 𝑖 = Net body force along i direction, m/L𝑡 2 , 𝑃𝑎 𝐺 = Shear modulus, m/L𝑡 2 , 𝑃𝑎 ℎ 𝑒𝑞 = Equivalent heat source term, 𝑚/𝑡 3 /𝐿, J/𝑚 3 /s ℎ 𝑔 = Heat source term in gas phase, 𝑚/𝑡 3 /𝐿, J/𝑚 3 /s ℎ 𝑚 = Heat source term in rock matrix, 𝑚/𝑡 3 /𝐿, J/𝑚 3 /s 𝐽 𝑖 = Characteristic adsorption energy for ith type of adsorption, m𝐿 2 /𝑡 2 /𝑛, J/𝑚𝑜𝑙 𝑘 𝑎 = Apparent gas permeability, 𝐿 2 , 𝑚 2 𝒌 𝒂 = Gas apparent permeability tensor, 𝐿 2 , 𝑚 2 𝑘 𝑓 = Fracture permeability, 𝐿 2 , 𝑚 2 𝒌 𝒇 = Fracture permeability tensor, 𝐿 2 , 𝑚 2 𝑘 𝑓,𝑖 = Fracture permeability at initial reservoir conditions, 𝐿 2 , 𝑚 2 𝑘 ∞ = Matrix intrinsic permeability, 𝐿 2 , 𝑚 2 𝑘 ∞0 = Matrix intrinsic permeability at reference conditions, 𝐿 2 , 𝑚 2 𝐾 𝑛 = Knudsen number 𝐿 = Characteristic length of flow path, L, 𝑚 𝑚 = Total gas content, m/𝐿 3 , 𝑘𝑔/𝑚 3 𝑚 𝑎𝑑 = Gas adsorption mass per unit volume, m/𝐿 3 , 𝑘𝑔/𝑚 3 𝑀 = Molecular weight, m/n, 𝑘𝑔/mol 𝒏 = Normal vector to fracture surface 𝑃 = Reservoir pressure, m/L𝑡 2 , 𝑃𝑎 𝑃 𝑒 = Effective confined pressure, m/L𝑡 2 , 𝑃𝑎 𝑃 𝐿 = Langmuir pressure, m/L𝑡 2 , 𝑃𝑎 𝑃 𝑝𝑟 = Pseudo-reduced pressure 𝑃 0 = Pressure at reference condition, m/L𝑡 2 , 𝑃𝑎 𝒒 𝒇 = Flow rate vector in the fracture per unit length, 𝐿 2 /t, 𝑚 3 /𝑠/𝑚 𝒒 𝒈 = Velocity vector of gas phase, L/t, 𝑚/s 𝑄 𝑓 = Mass source term in fracture, m/𝐿 3 /𝑡, 𝑘𝑔/𝑚 3 /𝑠 𝑄 𝑚 = Mass source term in matrix, m/𝐿 3 /𝑡, 𝑘𝑔/𝑚 3 /𝑠 𝑟 = Effective pore radius, L, 𝑚 𝑟 0 = Effective pore radius at reference conditions, L, 𝑚 𝑅 = Universal gas constant, m𝐿 2 /𝑡 2 /𝑛/𝑇, 8.3145 𝐽/𝑚𝑜𝑙/𝐾 𝑇 = Reservoir temperature, T, 𝐾 𝑇 𝑔 = Temperature of gas phase, T, K Coefficient for pore radius dependent intrinsic permeability 𝜆 𝑒𝑞 = Equivalent thermal conductivity, m𝐿/𝑡 3 /𝑇, W/m/K 𝜆 𝑔 = Thermal conductivity of gas phase, m𝐿/𝑡 3 /𝑇, W/m/K 𝜆 𝑚 = Thermal conductivity of rock matrix, m𝐿/𝑡 3 /𝑇, W/m/K 𝜇 𝑔 = Gas viscosity, m/Lt, Pa• s 𝜈 = Poisson's ratio 𝜌 𝑔 = Gas density, m/𝐿 3 , 𝑘𝑔/𝑚 3 𝜌 𝑔𝑠𝑡 = Gas density at standard condition, m/𝐿 3 , 𝑘𝑔/𝑚 3 𝜌 𝑚 = Matrix density, m/𝐿 3 , 𝑘𝑔/𝑚 3 (𝜌𝐶 𝑝 ) 𝑒𝑞 = Equivalent volumetric heat capacity, m/𝑡 2 /𝑇/𝐿, J/K/𝑚 3 𝜎 𝑖𝑗 = Effective stress, m/L𝑡 2 , 𝑃𝑎 𝜎 𝑚 = Mean effective stress, m/L𝑡 2 , 𝑃𝑎 𝜎 𝑚0 = Mean effective stress at reference conditions, m/L𝑡 2 , 𝑃𝑎 𝜙 𝑓 = Fracture porosity 𝜙 𝑚 = Matrix in situ porosity 𝜙 𝑚0

	𝑇 𝑚	= Temperature of rock matrix, T, K
	𝑇 𝑝𝑟	= Pseudo-reduced temperature
	𝑢 𝑖,𝑗	= Component of displacement, L, m
	𝑉 𝐿	= Langmuir volume, 𝐿 3 /m, 𝑚 3 /𝑘𝑔
	𝑍	= Gas deviation factor
	𝛼	= Biot's coefficient
	𝛼 1	= permeability enhancement coefficient
	𝛽	= Thermal expansion coefficient, 1/T, 1/K
	𝛿	= Thickness of gas adsorption layer, L, 𝑚
	𝛿 𝑖,𝑗	= Kronecker delta
	𝜀 𝑖𝑗	= Elastic strain
	𝜀 𝑘𝑘	= Volumetric strain
	𝜂	=

Appendix Mathematical Formulations of Thermal-Hydraulic-Mechanical Modeling in Fractured Shale Gas Reservoirs

The fully coupled THM physical process in shale gas reservoirs involves fluid flow within the formation matrix and fractures, shale gas adsorption and desorption, real gas properties and in-situ stress that affected by both pressure and temperature locally. In this section, the general governing equations of the coupled fluid flow, geomechanics and heat transfer process will be presented first, and then the pressure-temperature dependent variables will be discussed separately. All the equations presented in this section do not have any unit conversion parameters, any unit system can be applied (e.g., SI unit or Oil Field unit), as long as they are consistent.

Gas Flow in Formation Matrix and Fractures

Due to the extremely small pore radii in shale gas formations, Non-Darcy flow/Gas-Slippage can have a significant impact on gas transport in porous medium. Ideally, Non-Darcy flow mechanism can be modeled accurately using molecular physics by capturing the interaction between molecules and pore walls on nanoscale, however, this technique is not practical for modeling flow through shale on reservoirs scale. To make simulation feasible, it is desirable to integrate molecular flow behavior with the standard Darcy's equation so that these mechanisms can be captured without extremely intensive computational efforts. Non-Darcy flow based on matrix apparent permeability can be used to describe the gas flow rate within the shale matrix:

where 𝒒 𝑔 is the gas velocity vector, 𝜇 𝑔 is the gas viscosity, and 𝒌 𝑎 is the matrix apparent permeability tensor, which is pressure and temperature dependent. The continuity equation within shale gas formation can be written as:

where 𝑚 is the gas content per unit volume, 𝜌 𝑔 is gas density , 𝑄 𝑚 is the source term and 𝑡 is the generic time. The gas content 𝑚 is obtained from two contributions:

where ∅ 𝑚 is matrix porosity, 𝜌 𝑔 ∅ 𝑚 is the free gas mass in the shale pore space per unit volume of formation, while 𝑚 𝑎𝑑 is the adsorbed gas mass per unit volume of formation .

Normally, gas flow rate inside fractures in shale reservoirs are much smaller when compared to the conventional reservoir, due to extremely small shale matrix permeability, and turbulence flow is unlikely to happen under this circumstances, so Darcy's law is adequate enough to describe the flow behavior inside the fractures. Discrete Fracture Network (DFN) with tangential derivatives can be used to define the flow along the interior boundary representing fractures within the porous medium.

where 𝒒 𝑓 is the gas volumetric flow rate vector per unit length in the fracture, 𝒌 𝑓 is the fracture permeability tensor, 𝑑 𝑓 is fracture width and ∇ T P is the pressure gradient tangent to the fracture surface. The continuity equation along the fracture reflects the generic material balance within the fracture:

where ∅ f is the fracture porosity, and 𝑄 𝑓 is the mass source term, which can be calculated by adding the mass flow rate per unit volume from two fracture walls: where 𝒏 is the vector perpendicular to fracture surface.

DFN is treated as an internal boundary condition, and it has one less dimension than the simulation domain (e.g., if reservoir domain is 2D, DFN is modeled as 1D internal line-boundary; if reservoir domain is 3D, DFN is modeled as 2D internal faceboundary). This reflects the reality for fractures with small apertures, where flow in the width direction within the fracture is negligible.

Stress Equilibrium

The porous medium is assumed to be perfectly elastic so that no plastic deformation occurs. The constitutive equation can be expressed in terms of effective stress ( 𝜎 𝑖𝑗 , ) strain ( 𝜀 𝑖𝑗 ), pore pressure (P) and temperature (T):

where G is the shear modulus, E is the Young's modulus, 𝜈 is the Poisson's ratio, 𝜀 𝑘𝑘 represents the volumetric strain, 𝛿 𝑖,𝑗 is the Kronecker delta defined as 1 for 𝑖 = 𝑗 and 0 for 𝑖 ≠ 𝑗, 𝛽 is thermal expansion coefficient and 𝛼 is the Biot's effective stress coefficient, which is assumed to be 1 in this study. The strain-displacement relationship and equation of equilibrium are defined as:

where 𝑢 𝑖 and 𝐹 𝑖 are the components of displacement and net body force in the i-direction. Combining Eq. A.7-A.9, we have a modified Navier equation in terms of displacement under a combination of applied stress, pore pressure and temperature variations:

𝑢 𝑗,𝑖𝑖 -𝛼𝑃𝛿 𝑖,𝑗 + 𝛽ΔTE𝛿 𝑖,𝑗 + 𝐹 𝑖 = 0 (A. 10)

Heat Transfer in Porous Medium

Heat transfer process is governed by a thermal diffusion equation in an isotropic porous medium and the radiative effects, viscous dissipation, and work done by pressure changes are negligible. Considering an elemental volume of a porous medium we have, for the matrix:

and for the gas phase:

where C m is the heat capacity of the rock matrix, C p,g is heat capacity at constant pressure of the gas phase, λ is thermal conductivity, ℎ is the heat source term, ρ is the density and the subscripts m and g refer to formation matrix and gas phases respectively. Setting T m = T g = T by assuming there is local thermal equilibrium at the wall of pores, and adding Eq.( A.11) and (A.12) we have

where . 13 -4) (ρC p ) eq is the equivalent volumetric heat capacity at constant pressure, λ eq is the equivalent thermal conductivity and ℎ 𝑒𝑞 is the overall equivalent heat source. Typically, the heat capacity of most formation rock is at least hundred times larger than that of gas phase and formation porosity is normally less than 3% in shale formation, so (𝜌𝐶 𝑝 ) 𝑒𝑞 , λ 𝑒𝑞 and ℎ 𝑒𝑞 are dominated by the term of (1 -∅ 𝑚 )𝜌 𝑚 𝐶 𝑚 , (1 -∅ 𝑚 )λ 𝑚 and (1 -∅ 𝑚 )ℎ 𝑚 , respectively. In addition, the gas flow rate is constrained by extremely low matrix permeability and with low gas density, the term 𝜌 𝑔 𝐶 𝑝,𝑔 𝒒 𝑔 • ∇𝑇 diminishes, so the influence of heat transfer by convection is negligible, and the heat conduction mechanism in the formation matrix dominates the entire heat transfer process. Consequently Eq. (A.13-1) can be simplified as the following:

Overall, with appropriate initial and boundary conditions, Eq. (A.2), (A.5), (A.10) and (A.14) can well describe the fully coupled THM system though the pressure and temperature dependent variables in those equations.

Pressure-Temperature Dependent Variables

In the following sections, the pressure and temperature dependent variables that needed for coupling process in the above governing equations will be determined, using either theoretical modeling or well-established correlations.

Real Gas Properties

The in-situ gas density can be calculated by real gas law:

where M is the average molecular weight of mixed gas, R is the universal gas constant. The Z-factor can be estimated by solving Equation of State (EOS) or using correlations for the gas mixtures. In this study, the Z-factor is calculated using an explicit correlation [START_REF] Mahmoud | Development of a New Correlation of Gas Compressibility Factor (Z-Factor) for High Pressure Gas Reservoirs[END_REF] based on the pseudo-reduced pressure (p pr ) and pseudo-reduced temperature (T pr ):

The advantage of using explicit correlation is to avoid solving higher order equations respect to Z-factor, which leads to multiple solutions and increases computation efforts.

Among all gas properties, gas viscosity, which is usually characterized by some available correlations, still remains uncertain. [START_REF] Lee | Adsorption Studies of Natural Gas Storage in Devonian Shales[END_REF] proposed a correlation to calculate gas viscosity at temperatures from 310 to 445 K and pressure from 0.69 to 55.16 MPa. [START_REF] Viswanathan | Viscosities of Natural Gases at High Pressures and High Temperatures[END_REF] measured the viscosity of pure methane at pressure from 34.47 to 206.84 MPa and temperatures from 310 to 478 K, and modified the gas viscosity correlation by [START_REF] Lee | Adsorption Studies of Natural Gas Storage in Devonian Shales[END_REF]:

where 𝐾 = (5.0512 -0.2888𝑀)𝑇 1.832 -443.8 + 12.9𝑀 + 𝑇 (A. 17 -2) 𝑋 = -6.1166 + ( 3084.9437 𝑇 ) + 0.3938𝑀 (A. 17 -3) 𝑌 = 0.5893 + 0.1563𝑋 (A. 17 -4)

Gas Adsorption Capacity

Langmuir isotherms (1916) are widely used to model adsorption gas content in shales as a function of pressure, which can be established for a prospective area of shale basin using available data on TOC and on thermal maturity to establish the Langmuir volume (V L ) and the Langmuir pressure (P L ). And adsorbed gas-in-place is then calculated using the formula below:

where ρ gst is gas density at standard condition. However, the temperature effect on gas adsorption capacity is not included in Langmuir isotherms. Lu et al. (1995) investigated temperature dependent adsorption curves on shale samples and proposed a Bi-Langmuir model that accounts for gas adsorption on both clay minerals and kerogen:

where 𝑓 𝑖 is defined as the ratio of the amount of the ith type of adsorption at monolayer coverage to the total amount adsorbed at monolayer coverage and each adsorption site is assumed to follow the Langmuir equation. For two types of adsorption sites (clay minerals and kerogen) we have: . 19 -4) where 𝑎 1 and 𝑎 2 is a pre-exponential constant independent of temperature, 𝐽 1 and 𝐽 2 are the characteristic adsorption energy. The five unknown independent parameters 𝑓 1 , 𝑎 1 , 𝑎 2 , 𝐽 1 and 𝐽 2 in Eq.( A.19) can be determined from a Langmuir isotherm curve at provided temperature condition by non-linear regression. Thus, the temperature effects can be included to describe shale gas adsorption capacity as a function of both pressure and temperature (as shown in Fig. 6), with available Langmuir volume and Langmuir pressure values as input parameters [START_REF] Wang | Conceptual Study of Thermal Stimulation in Shale Gas Formations[END_REF], or directly fitting and validated against experiment data within possible temperature ranges [START_REF] Yue | Increasing Shale Gas Recovery through Thermal Stimulation: Analysis and an Experimental Study[END_REF].

Fracture Pressure Dependent Permeability

It is a well-known fact that most shale formations have massive pre-existing natural fracture networks that are generally sealed by precipitated materials weakly bonded with mineralization. Such poorly sealed natural fractures are generally reported to interact heavily with the hydraulic fractures during the injection treatments, serving as preferential paths for the growth of complex fracture network. The reactivation of such pre-existing planes of weakness (i.e., natural fractures, micro-fractures, fissures) is well documented and observed in microseismic monitoring [START_REF] Cipolla | A Practical Guide to Interpreting Microseismic Measurements[END_REF][START_REF] Zakhour | Real-Time Use of Microseismic Monitoring for Horizontal Completion Optimization Across a Major Fault in the Eagle Ford Formation[END_REF]. Ideally, the fracture conductivity on un-propped and propped fracture surfaces, as a function of confining stress, fluid type, and proppant type should be investigated for a specific reservoir in order to optimize stimulation design and field development, and the fracture pressure-dependent permeability (PDP) data should be incorporated into reservoir simulation model to better forecast the well performance under various operation conditions. The expected reduction in fracture permeability caused by the increasing effective stress during production is can be described by power law relationship [START_REF] Cho | Pressure-Dependent Natural-Fracture Permeability in Shale and Its Effect on Shale-Gas Well Production[END_REF], here, we use the following correlation:

where 𝑘 𝑓 is the fracture permeability, 𝑘 𝑓,𝑖 is the fracture permeability at initial reservoir conditions, B is a fracture compaction parameter that can be determined from experimental data, and 𝜎 𝑚 is the mean effective stress, which is influenced by declining pressure and rising temperature during stimulation.

Correction for Gas Apparent Permeability in Nano-Pores

Darcy's law cannot describe the actual gas behavior and transport phenomena in nano-porous media. In such nanopore structure, fluid flow departs from the well-understood continuum regime, in favor of other mechanisms such as slip, transition, and free molecular conditions. The Knudsen number [START_REF] Knudsen | Die Gesetze der Molekularstromung und der inneren Reibungsstromung der Gase durch Roḧren (The laws of molecular and viscous flow of gases through tubes)[END_REF] is a dimensionless parameter that can be used to differentiate flow regimes in conduits at micro and nanoscale. For conduit with radius r, Knudsen number can be estimated by

Table A1 shows how these different flow regimes, which correspond to specific flow mechanisms, can be classified by different ranges of 𝐾 𝑛 .

𝐾 𝑛 0 -10 -3 10 -3 -10 -1 10 -1 -10 The apparent permeability of shale matrix can be represented by the following general form:

where 𝑘 ∞ is the intrinsic permeability of the porous medium, which is defined as the permeability for a viscous, non-reacting ideal liquid, and it is determined by the nano-pore structure of porous medium itself. 𝑓(𝐾 𝑛 ) is the correlation term that relates the matrix apparent permeability and intrinsic permeability. Sakhaee-Pour and Bryant (2012) developed correlations which are based on lab experiments. They proposed a first-order permeability model in the slip regime and a polynomial form for the permeability enhancement in the transition regime using regression method:

𝑓(𝐾 𝑛 ) = { 1 + 𝛼 1 𝐾 𝑛 Slip Regime 0.8453 + 5.4576𝐾 𝑛 + 0.1633𝐾 𝑛 2 Transition Regime 0.1 < 𝐾 𝑛 < 0.8 (A. 23) where 𝛼 1 is permeability enhancement coefficient in slip regime. To ensure the approximation of continuity of 𝑓(𝐾 𝑛 ) at the boundary region of slip and transition regime, where no existing model available, 𝛼 1 is set to be 4. For all the simulation cases presented in this study, 𝐾 𝑛 is always larger than 0.1, so only the polynomial expression for transition regime is automatically used in our simulator.

Intuitively, Eqs. (A.21) and (A.23) predict net increases of 𝐾 𝑛 and 𝑓(𝐾 𝑛 ) with decreasing pore pressure, respectively. However, decreasing pore pressure can also lead to reduction in pore radii and, in turn, reduce the intrinsic permeability 𝑘 ∞ . It can be seen from Eq. (A.22) that the evolution of matrix apparent permeability 𝑘 𝑎 is determined by two combined mechanisms (i.e., the variations in 𝑘 ∞ and 𝑓(𝐾 𝑛 ) during production).

Even though many studies (Lunati and Lee 2014; [START_REF] Sigal | Examination of the Importance of Self Diffusion In the Transportation of Gas In Shale Gas Reservoirs[END_REF][START_REF] Kazemi | An analytical model for shale gas permeability[END_REF] have developed various theoretical models to account for the effects of permeability enhancement in nano-pore structure of shale gas, but how to quantify these effects in complex, heterogeneous shale systems on a field scale is still not well understood within scientific and industry community. In addition, all these proposed models do not account for the dynamic changes of stress and adsorption layer as pressure declines (as reflected in Fig. 2) in a fractured system, which can significantly overestimate the effects of permeability enhancement on production and production decline trend (Wang 2016c). Recent study (Wang et al 2015b) also reveals that under most real shale-gas-reservoir conditions, gas adsorption and the non-Darcy flow are dominant mechanisms in affecting apparent permeability, while the effects of surface diffusion caused by adsorbed gas can be largely overlooked. So in this article, molecule diffusion will not be discussed in the following derivations for matrix apparent permeability.

Wang and Marongiu-Porcu (2015) conducted a comprehensive literature review on gas flow behavior in shale nano-pore space and proposed a unified matrix apparent permeability model, which bridges the effects of geomechanics, non-Darcy flow and gas adsorption layer into a single mode, by considering the microstructure changes in nano-pore space. In a general porous media, the loss in cross-section area is equivalent to the loss of porosity,

where variables with a subscript 0 correspond to their value at the reference state, which can be laboratory or initial reservoir conditions. Laboratory measurements by [START_REF] Dong | Stress-Dependence of the Permeability and Porosity of Sandstone and Shale from TCDP Hole-A[END_REF] shows that the relationship between porosity and stress follows a power law relationship, and can be expressed using the concept of mean effective stress 𝜎 𝑚 [START_REF] Wang | Impact of Shale Gas Apparent Permeability on Production: Combined Effects of Non-Darcy Flow/Gas-Slippage, Desorption and Geomechanics. SPE-173196-PA[END_REF]:

𝐶 ∅ is a dimensionless material-specific constant that can be determined by lab experiments. For silty-shale samples, the values of 𝐶 ∅ range from 0.014 to 0.056 [START_REF] Dong | Stress-Dependence of the Permeability and Porosity of Sandstone and Shale from TCDP Hole-A[END_REF]. Combining Eq. (A. 24) and Eq. (A. 25) leads to the relationship between pore radius and local stress r = 𝑟 0 ( 𝜎 𝑚 𝜎 𝑚0 ) -0.5𝐶 ∅ (A. 26)

When the adsorption layer is considered, the thickness of the gas adsorption layer, δ ,can be interpolated based on a Langmuir type functional relationship: δ = 𝑑 𝑚 𝑃/𝑃 𝐿 1 + 𝑃/𝑃 𝐿 (A. 27) where 𝑑 𝑚 is the average diameter of gas molecules residing on the pore surface. And the effective pore radius (A. 26) can be modified as: r = 𝑟 0 ( 𝜎 𝑚 𝜎 𝑚0 )

-0.5𝐶 ∅ -𝑑 𝑚 𝑃/𝑃 𝐿 1 + 𝑃/𝑃 𝐿 (A. 28) The relationship between intrinsic permeability and pore radius can have the following relationship