Sustainable Production from Shale Gas Resources through Heat-Assisted Depletion

: Advancements in drilling and production technologies have made exploiting resources, which for long time were labeled unproducible such as shales, as economically feasible. In particular, lateral drilling coupled with hydraulic fracturing has created means for hydrocarbons to be transported from the shale matrix through the stimulated network of microcracks, natural fractures, and hydraulic fractures to the wellbore. Because of the degree of conﬁnement, the ultimate recovery is just a small fraction of the total hydrocarbons in place. Our aim was to investigate how augmented pressure gradient through hydraulic fracturing when coupled with another derive mechanism such as heating can improve the overall recovery for more sustainable exploitation of unconventional resources. Knowledge on how hydrocarbons are stored and transported within the shale matrix is uncertain. Shale matrix, which consists of organic and inorganic constituents, have pore sizes of few nanometers, a degree of conﬁnement at which our typical reservoir engineering models break down. These intricacies hinder any thorough investigations of hydrocarbon production from shale matrix under the inﬂuence of pressure and thermal gradients. Kerogen, which represents the solid part of the organic materials in shales, serves as form of nanoporous media, where hydrocarbons are stored and then expelled after shale stimulation procedure. In this work, a computational representation of a kerogen–hydrocarbon system was replicated to study the depletion process under coupled mechanisms of pressure and temperature. The extent of production enhancement because of increasing temperature was shown. Moreover, heating requirements to achieve the enhancement at reservoir scale was also presented to assess the sustainability of the proposed method.


Introduction
Thermal enhanced recovery of conventional resources through steam injection, electrical heating, and in situ combustions have shown significant improvement on oil recovery. The process of heating reservoirs is primarily meant to enhance the oil mobility which makes the approach almost exclusively applied for heavy oil. However, heating can play role in altering some other properties for the favor of enhanced production. In unconventional resources of shale gas and oil, a large amount of hydrocarbons is stored in the organic constituents of the shale matrix known as kerogen [1]. Hydrocarbons are produced from those confinements by the continuous expansion of free compressed fluids and desorption of sorbed molecules through non-Darcian form of transport [2][3][4][5][6]. The modeling approach usually incorporates mechanisms such as adsorption, slip flow, and molecular and surface diffusions. These are pressure-and temperature-dependent quantities. Unconventional resources are produced under the influence of augmented pressure gradient through hydraulic fracturing under isothermal conditions (i.e., reservoir temperature being constant throughout production time). However, increasing the temperature through means of adding heat to the reservoir may alter some of the aforementioned transport and storage mechanisms for enhanced production. In this work, our aim was to investigate how coupling pressure depletion with heating can impact hydrocarbon recovery from organic materials of shales through a molecular computational approach. In addition, a reservoir scale design of hydraulically fractured reservoir assisted with heating was investigated for the power requirements and heat propagation.

Limitations of Macroscopic Equations
Assuming monolayer adsorption and adsorbing sites being organized in perfectly flat planes, the Langmuir isotherm is probably the most commonly utilized theoretical model to study the sorption behavior of hydrocarbons and carbon dioxide in shale or coal bed methane reservoirs, respectively [7][8][9]. The equilibrium for a single adsorbate and the adsorbent system is given as where V abs (scf/ton) constitutes the total or absolute amount of gas adsorbed in the adsorption space; V L describes the Langmuir volume (scf/ton) corresponding to the maximum theoretical sorption capacity achievable under an infinite pressure P (psi) assumption; and P L (psi) is the Langmuir or critical desorption pressure, defined as the pressure where the volume of adsorbed gas equals fifty percent of the Langmuir volume [10,11]. It has to be pointed out that, from an experimental point of view, and particularly for high pressure experiments, determination of the absolute adsorbent amount is problematic. A more pragmatic approach is to communicate the incremental change in the form of, for instance, the surface excess adsorbed amount [12][13][14]. Given our particular interest in multicomponent mixtures, an extended form of the Langmuir isotherm model can be defined as follows [15][16][17][18]: In Equation (2), p i represents the partial pressure of the i-th component. Consequently, the total amount of gas adsorbed is given with The extended Langmuir model, as shown in Equation (3), seeks to overcome limitations of the original model but is unable to adequately account for adsorbate size effects resulting from dissimilarities of saturation capacities for pure species Langmuir isotherms [19]. This is particularly true for multicomponent sorption studies with respect to shales, where the sorption of hydrocarbons on kerogen, a complex and heterogeneous carbonaceous material, is subject to ongoing research [20][21][22][23]. Furthermore, both classical and extended Langmuir formulations cannot accommodate sorption phenomena at varying temperatures, as encountered during thermal stimulation-or gas injection-based enhanced hydrocarbon recovery techniques [24].
Fluids confined in nanoporous media behave differently. Luo et al. [25] conducted an experimental work to study the phase behavior of octane-decane binary mixture. The bubble point in their work was found to be lower than its typical values. The extent of that change in the phase behavior was found to be proportional to the degree of confinement. Walton et al. [26] used the grand canonical Monte Carlo simulation (GCMC) to study the condensation of confined fluids and found it to occur at a pressure lower than what was expected. Wongkoblap et al. [27] followed a similar approach and reported similar observations with some pore size-dependent hysteresis in the condensation/evaporation processes.
Jiang et al. [28] studied the phase transitions of hydrocarbons in nanotubes and found their critical temperature and pressure to be different than their bulk phase values. Singh et al. [29] observed a positive deviation of saturation pressure of confined fluid for a 3 nm slit-pore and negative deviation when the pore size was reduced to 1 nm.
The degree of uncertainty is even larger when more than one component is present. Didar and Akkutlu [30] and Pitakbunkate et al. [31] studied a mixture of hydrocarbons where one component had higher relative affinity toward the surfaces of nanotubes. They found the properties of the mixture to be highly correlated with the degree of confinement. Bui and Akkutlu [32] conducted a combined molecular dynamics (MD) and GCMC to study the compositional variation of a mixture of hydrocarbons in a nanotube that was in thermodynamic equilibrium with a bulk phase of constant composition. They found the composition inside the nanotube to be different than the bulk phase. Moreover, the variation of composition inside the nanotube was found to be dependent on the pressure. The fraction of heavier molecules increased as the pressure decreased. That would have some significant implications on the properties of the fluids, and hence the recovery from shale reservoirs. In this study, we considered a group of hydrocarbons confined in a representative kerogen structure undergoing a continuous depletion at (1) typical reservoir temperature and (2) increased temperature. We were interested in understanding the influence of increasing the temperature on the overall productivity.
In general, macroscopic modeling of fluid storage and transport would breakdown in nano-confinements. Typical equations are insufficient to tell apart sorbed from free fluid nor capturing selectivity of kerogen nanopores of given component over the others under given conditions of pressure and temperature. For that, a molecular computational approach was selected to investigate the system of hydrocarbons confined in a representative kerogen model for better understanding of how variables such as temperature and pressure can influence the recovery from such resources.

Molecular Simulation Approach
The build-up of pressure inside kerogen pores during hydrocarbons generation results in the formation of microcracks. This can be seen in various electron microscopy studies of shale matrix [33][34][35]. These microcracks serve as conduits of fluid transport from the organic nanopores to the larger size fractures [5]. During production, compressed hydrocarbons undergo continuous expansion and depletion through microcracks, shifting the adsorption equilibrium to the favor of desorption. The walls of the organic nanopores influence the transport mechanisms, making the continuum approach of modeling fluid transport insufficient due to the degree of confinement. Various studies have considered the transport of natural gas in nanopores and derived composite advective-diffusive-adsorptive models [3,4,36,37]. Moreover, the surfaces of kerogen can favor one component over the other, which makes the composition of a given mixture vulnerable to continuous changes during the production span.
We built our conceptual model on the basis of the repeated observations of the close associations of organic materials and microcracks. Examples of that are the SEM images shown in Figure 1. Hydrocarbons stored in the kerogen were subject to pressure gradient when the formation was hydraulically fractured. During their continuous migration to microcracks, hydrocarbons were subject to chemical and physical interactions with the kerogen surfaces. A molecular replication of kerogen was obtained, and then loaded with some arbitrary multicomponent gas mixture. After that, desorption simulations at typical reservoir temperature were performed to mimic the depletion process, as shown in Figure 2. The same procedure was repeated at elevated temperature. The goal was to closely monitor compositional variations during the depletion of shale matrix for both scenarios. This can provide some insights to how heat-assisted depletion can improve overall productivity.

Kerogen Models
Ungerer et al. [38] created six kerogen molecules corresponding to different levels of maturity. These can help as the building blocks of the nanoporous organic materials. In their study, type I was built for immature kerogen. Four versions of type II were created with some variations in their origins. Type III is represented by one version that corresponds to coal. These porotypes, which were created in a computational platform, were found to be within acceptable compositional agreement with experimental work. In our work, type II D, which corresponds to shale gas, was recreated and used for our study.
In this section, the molecular modeling approach used in this work is presented. The first subsection details building the molecular structures of the type II-D kerogen. The second part concerns the creation of realistic models for condensed kerogen bulk phases including the characterization and porosity estimation. The last part discusses the methodology in estimating the sorption of multicomponent systems into a condensed kerogen phase. The model building and simulation tools employed are those provided within the MedeA simulation environment. Specifically, these tools provide for the creation of realistic model configurations of condensed kerogen phases, which by nature are an amorhous ensemble of kerogen molecules. For the molecular mechanics and dynamics simulation, the MedeA environment makes use of the powerful program LAMMPS, developed by Plimpton and coworkers at Sandia National Laboratories [39]. For the Monte Carlo simulations in the isothermal-isochoric NVT ensemble, the MedeA environment employs the GIBBS program [40].
Within the MedeA software environment, input models and command scripts are prepared and submitted to LAMMPS and GIBBS, with raw data returned to MedeA for further processing including, in the present work, statistical analysis of calculated properties of kerogen phases, such as (among others) the density, fugacity, chemical potential, and amount of adsorbed and desorped molecules. All calculations used the polymer consistent forcefield PCFF+, which is based on the class II forcefields included with the LAMMPS distribution, following extensive refinement of nonbond parameters using liquid state PVT and cohesive property data for small molecules.
The building process for the model of the condense bulk phase of the kerogen type II-D encompassed several steps. Initially, five molecules of the kerogen type II-D were placed inside a cubic simulation cell, such that the density was 0.1 g/cc. Then a low energy configuration of randomly oriented kerogen molecules was created, using a Monte Carlo algorithm as implemented in the Amorphous Materials module of MedeA. To build a dense kerogen bulk phase, the resulting structure was relaxed in a sequence of several molecular dynamics (MD) simulations, such as illustrated in Figure 3. The first MD simulation was performed within the canomical (NVT) ensemble. The initial structure of randomly oriented kerogen molecules was headed up to a temperature of 900 K for 250 ps. The next four MD simulations were performed within in the isobaric (NPT) ensemble to condition the density of the bulk kerogen phase for a particular pressure. Stepwise, the temperature was decreased from 900 to 300 K at a constant pressure of 20 MPa (2900 psi), as shown in Figure 3. Within the MD simulations, temperature and pressure were controlled with the Nosé-Hoover thermostat and barostat, respectively [41]. The time step used in all calculations was 1 fs.
At the final conditions of pressure and temperature, the density was calculated to be around 1-1.2 g/cc. The porosity of the final structure was estimated by a stimulated helium pycnometry using the program GIBBS. As in the experimental helium pycnometry method, in the simulation the effective pore volume was also determined by adsorption of helium [42,43]. The calculated porosity of the final bulk phase of the kerogen type II-D was 19%.

Results and Discussions
The nanoporous material was loaded with a mixture of C 1 -C 4 at a pressure of 2500 psi and temperature equal to 194 • F (i.e., typical reservoir temperature). The bulk phase composition was arbitrary selected to be 0.3, 0.15, 0.2, and 0.35 for the four alkanes, respectively. The initial composition of the mixture was 0.673, 0.255, 0.0073, and 0.065, respectively, as shown in Table 1. The sorption of the different components was obtained by using Grand Canonical ensemble [44]. The parameters of such ensemble are the adsorbed molecule chemical potential, volume, and temperature. The fugaicties were imposed for each component (i.e., its gas mixture with relative high pressure), which was the equivalent of using the chemical potential as follows: where µ i , µ io , R, T, f i , and f io are the chemical potential for component i, ideal chemical potential for component i, gas constant, temperature, fugacity of component i, and ideal fugacity of component i, respectively. All the adsorbent components were treated as rigid, except C 1 due to the high fugacity value. Rigid type simplified the problem because there was no need to consider the solid-solid interaction energy [40]. The kerogen type II-D molecule shown in Figure 3 was recreated at the molecular platform and the nanoporous material was created by following the description provided in the previous section. Porosity was selected as 19.6% by adding 10 molecules of heptane during the construction phase and then deleting them afterwards. The initial composition of components in the nanoporous structure was obtained as shown in Table 2. Desorption of hydrocarbons were carried out at a typical reservoir temperature of 194 • C and then at 284 • C. During the desorption stages, number of molecules of each component was monitored for further analysis. The total number of desorbed molecules for both cases is shown in Figure 4. The recovery factor RF is calculated as where n is the total number of molecules present at a given depletion stage and n i is the initial number of molecules present at initial conditions. The ultimate recovery factor at the last desorption was 36%, whereas that of the heated case was almost 70%, indicating an improved recovery by two times. The compositional variations at each pressure stage were obtained for typical reservoir and heated temperature cases, as shown in Figures 5 and 6, respectively. For the typical temperature, the lightest component which was methane was selectively desorbed, leaving other heavier molecules behind. This behavior was referred to as the fractionation effect, where continuous compositional variations were observed. It can be seen from Figure 5 that gas composition is maintained constant above certain pressure, which is in this case was around 1500 psi (this pressure is called saturation pressure, which is the maximum adsorption capacity pressure). The structure above this pressure tended to hold on to its contained gases. However, the extent of that fractionation became more pronounced below saturation pressure (i.e., maximum adsorption capacity pressure). On the other hand, the fractionation effect was less severe for the heated case, as components maintained almost constant composition during desorption stages. The influence of compositional variations on the phase behavior was investigated by obtaining a P-T diagram for given composition of both cases at two stages of depletion, as shown in Figures 7 and 8, respectively. It can be seen that the two-phase envelope was larger for the unheated case compared with the case of depletion at a higher temperature. The two-phase envelope of the unheated scenario continued to expand as kerogen pressure dropped further. This effect was attributed to the accumulation of heavier components. On the other hand, the two-phase envelope of the heat-assisted depletion almost stayed constant as all components were expelled from the structure uniformly. Hence, heating can reduce fractionation effect and potentially condensate problems to maximize the overall recovery.   . P-T diagram generated at the gas composition of typical (orange) and heated (blue) when the kerogen pressure dropped to 1000 psi. The typical reservoir temperature had a larger two-phase envelope due to accumulation of heavier components. Figure 8. P-T diagram generated at the gas composition of typical (orange) and heated (blue) when the kerogen pressure dropped to 500 psi. The two-phase envelope continued to increase for the typical reservoir temperature while maintaining its shape for the heated case.

Reservoir Scale Heat Transfer Model
In this section, the feasibility of heating an unconventional reservoir by either injecting hot gas or by introducing a heating element was investigated. The model assumed transverse fractures intersecting the wellbore lateral created by multistage hydraulic fracturing (MSF), as shown in Figure 9. It also assumed that a 30 ft spacing between hydraulic fractures along a 4500 ft long lateral was designed, producing 150 transverse fractures. The fracture half-length was 400 ft, whereas the fracture height was 150 ft. Due to symmetry, the simulation domain was half of the stimulated reservoir volume (SRV), as shown in Figure 9. A coupled heat and mass transfer model was used to study the heat propagation while injecting hot fluids. The detailed mathematical description is given in the Appendix A. The reservoir temperature was around 200 • F, whereas the gas was heated to 400 • F before injection. The thermal properties of the reservoir rock and fluids are shown in Table 3. The gas injector was simulated in this case to understand the heat propagation, given that the producer was neglected. It was found that the hot gas can advance heat to a sufficient distance only when a considerable amount was being injected for long periods. Figure 10 shows the heat propagation at a steady-state condition when injecting 0.5 MMscf/D and 1.0 MMscf/D of hot gas. The hydraulic fracture (one wing) was placed in the middle of the x-axis (0 location). As expected, most of the heating occurred within the fracture surroundings. Additionally, the larger the injection rate, the longer the heat penetration distance.  Hot gas injection was not feasible for many reasons. The process of gas injection is very expensive and could cause severe corrosion when CO 2 is the selected gas. Additionally, the practice of heating gas while injecting and separating it after production is costly and unpractical. A considerable amount of gas must be injected; despite this, the heat does not fully propagate to the tip of the fracture. A large amount of heat could be lost around the wellbore, especially if it has to travel a long distance before reaching the perforation clusters. Figure 11 shows a comparison between the initial geothermal temperature and after injection temperature around the vertical wellbore section. Notice that the wellbore was placed on the west edge of the figure. As can be observed from the temperature magnitudes, a large amount of heat was lost to the formation before the hot fluids reached the fracture. The other option would be placing a heated element around the lateral section of the wellbore. This would assume that the heated element was permanently placed at a constant temperature of 400 • F while the reservoir was under production (10 MMscf/D). The heating element was considered as a constant temperature boundary condition at the south side of the domain, as shown in Figure 12. The feasibility of the heated element arises from two facts: first, the heating propagation from the element into the reservoir does not act against the formation linear advection direction. Second, heat propagation, in this case, is transient and kept reaching new fronts as time progress (see Figure 12). As a matter of fact, heat propagation did not terminate while the heated element was operating. The mechanism at which the heat was propagated was through conduction within the formation's rock and fluids. In contrary to the previous scenario, the heat propagated almost equally through the SRV, except around the producing fractures. The reason was that the fluids were advected against the heat propagation only within the fracture.
Besides the heating time, two other factors can influence the heat propagation. First, the advection rate can significantly enhance heat transfer within the SRV, causing the heat to reach longer distances. Figure 13a shows the temperature map when the lateral was assumed to produce at 100 MMscf/D. The heat reached around 90 ft along the fracture direction, whereas only 55 ft for the 10 MMscf/D production case (see Figure 12b). Additionally, increasing the spacing between fractures can enhance heat propagation, as can be seen in Figure 13b. This is because an SRV will receive more heat flux, which improves the efficiency of heat transfer.

Conclusions
The common practice of producing unconventional resources through hydraulic fracturing can be enhanced when coupled with heating. The effect of temperature change on recovery was investigated for kerogen where significant amounts of hydrocarbons are stored. The computational work was carried out on a molecular replication of a kerogen sample using molecular simulations. The results showed an improved recovery by a factor of two with reduced fractionation effect in comparison with the non-heated case for the same applied pressure difference when the temperature was increased by 90 • F. Furthermore, the energy requirements to achieve that temperature rise at reservoir scale was investigated. Desired reservoir temperature for enhanced recovery could be achieved through resistive heating or fluid injection. Acknowledgments: Authors gratefully acknowledge College of Petroleum Engineering and Geoscience CPG support.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
An integrated flow and heat transfer model was developed for a system containing a wellbore, fracture, and reservoir. These were integrated through the shared boundary conditions. The mathematical description of each model is shown below. Detailed description of the model was provided by Aljawad et al. [45].
The wellbore model is transient one-dimensional (1D), and projected to estimate the temperature of the fluid incoming to the fracture from the wellbore. The wellbore model is numerically coupled with a reservoir model to estimate the heat transfer from the wellbore to the surrounding formation. To obtain the velocity profile in the wellbore, the following continuity equation is solved: where t is injection time, R w is the wellbore radius, ρ f is the injected fluid density, γ is the wellbore open ratio, v L,p is the velocity of fluids entering the fracture, and z is the direction along the wellbore length. Obtaining the velocity profile makes it possible to solve the wellbore energy balance equation. The heat transfer equation within the wellbore can be written as where T wb is the temperature of the wellbore,Ĉp f is the fluid's specific heat capacity , T f |B is the fracture temperature at the wellbore/fracture boundary, U is the overall heat transfer coefficient between the wellbore and formation, and T r|B is the temperature at the wellbore/formation (i.e., reservoir) boundary. The subscript B is used to represent the value at the boundary of the wellbore, reservoir, and fracture. Solving Equation (A2) requires information about the fracture and formation temperatures. A procedure that requires iteration was implemented to solve the complete system. The geothermal temperature was considered as the initial condition. The radial reservoir heat conduction model was coupled along the wellbore, which can be written as where ρĈp is the average reservoir rock and fluid property, k e is the effective average thermal conductivity, and r is the radial direction away from the wellbore. During injection, the fluid loss from the fracture induces fluid flow in the reservoir. The flow in the coupled fracture/reservoir system can be estimated by solving diffusivity equation for compressible fluids, written as ∇ · (k nD ∇m(p)) = ϕµc t ∂m(p) ∂t (A4) where k nD is the non-Darcy permeability, ϕ is the formation porosity, c t is the total compressibility, and m(p) is the pseudo-pressure function. That enables solving the heat transfer model in the reservoir/fracture system, which can be described as ρ fĈ p f ∂T r ∂t + −ϕβ T T r ∂p ∂t + ρ fĈ p f v · ∇T r = ∇ · k f ∇T r + (β T T r − 1)v · ∇p (A5) where β T is the thermal expansion factor, v is the velocity vector, k f is the fluid's thermal conductivity vector, and p is the formation pressure. The first term in the above equation represents the accumulation of heat in the fracture, the second term is heat convection, and the last term is heat conduction.