Thermal Diffusion Characteristics in Permafrost during the Exploitation of Gas Hydrate

Methane hydrate is the vast potential resources of natural gas in the permafrost and marine areas. Due to the occurrence of phase transition, the gas hydrate is dissociated into gas and water and absorbs lots of heat. The incomprehensive knowledge of endothermic reaction in permafrost sediments still restricted the production efficiency of hydrate commercial development. This endothermic reaction leads to a complex thermal diffusion in permafrost, which directly influences the phase transition in turn. In this research, the heat during the exploitation is transferred in two forms (specific heat and latent heat). Besides, the melting point is not constant but depends on the pore size of the reservoir rock. According to these features, a thermal diffusion model with phase transition is established. To calculate the governing equation, the pore size distribution is obtained by using the nuclear magnetic resonance (NMR) method. The heating tests are conducted and simulated to calibrate the coefficient (i.e., transverse surface relaxivity) of NMR. Then, the temperature field evolution of the hydrate reservoir during the exploitation is simulated by using the calibrated values. The results show that the temperature curves have a typical plateau related to the pore size distribution, which is effective to obtain the surface relaxivity. The heat transfer is remarkably limited by the endothermic effect of the phase transition. The hydrate recovery efficiency may depend largely on the heating capacity of the engineering operation and the rate of gas production. Compared to the conventional petroleum industry, it is significant to control the maximum temperature and temperature distribution in engineering operations during hydrate development. This research on the temperature behavior during onshore permafrost hydrate production could provide the theoretical support to control heat behavior of offshore hydrate production.


Introduction
Gas hydrate is a solid crystalline material formed by digenetic crystallization and among gas resources [1][2][3][4]. It generally exists in permafrost regions and in marine sediments with the condition of high pressure and low temperature [5][6][7]. Over 95% of the gas hydrate resource is found in sediments below seafloors, and the amount of carbon bound in gas hydrates is estimated to total twice that in all known fossil fuels on earth [8]. Due to the high carbon-oxygen, gas hydrate is commonly recognized as the potential energy source for replacing traditional fossil fuels. The potentially largest natural gas resource remaining on earth, oceanic natural gas hydrate, may substantially supplement the natural gas supply far into the future. US, Japan, Canada, and China had conducted the hydrate trial tests which promoted the hydrate commercial hydrate process [9][10][11][12].
There are three main methods for gas production: depressurization, thermal stimulation, and use of inhibitors [13][14][15][16][17]. All methods are used to break down the hydrate phase and collect the natural gas. Because the occurrence environment of hydrate (high pressure and low temperature) is different from the conventional petroleum geology environment, once the balance is being destroyed, the gas hydrates become instable and decompose [18,19]. Moreover, it should be noted that ice and hydrate in the permafrost maintain the stability of the skeleton under the condition of perennial low temperature [20]. The original balance of permafrost reservoir will be broken by operation of the drilling, well completion, well cementation process, and production, like operation fluids will increase the temperature of reservoirs and melt the ice skeleton [19,[21][22][23][24][25]. The dissociation of hydrate skeleton by increased temperature and production also will reduce the strength of the soil skeleton. From the field experience of Canada, Russia, and US hydrate production in the permafrost region ( Figure 1) [9,11,12,26], the hydrate production operation will change the temperature of reservoir skeleton (rock, ice, and hydrate). This temperature change will dramatically affect the efficiency of hydrate production [27,28]. It may conduct the wellbore instability, sand production ( Figure 2), and methane leakage [26,[29][30][31]. The temperature range of hydrate layers in Messoyakha gas field, Mount Elbert gas hydrate well (Alaska), Malik gas hydrate site (Canada), and Qinghai-Tibet Plateau were about 8~11°C [32], -8~10°C [33,34], -20~-1°C [6,12,33], and -0.1~-5°C [35], respectively, which are different from the traditional petroleum temperature. Therefore, it is necessary to study the disturbance in the temperature field of the hydrate reservoirs in permafrost.
To study the disturbance in the temperature field of hydrate reservoir, the heat transfer process of the reservoir must be first understood. The heat is not only converted into specific heat, causing the variation in reservoir temperature. At the same time, the transformation of latent heat will consume a lot of heat. This change in heat, in turn, affects the occurrence of phase transitions, which is a complex coupling process. It should be noted that the process of heat transfer with phase transition must consider the influence of the pore structure of the rock. In rocks, the melting point is no longer a constant but is related to the pore size of the rock. The smaller the pore size of the rock, the lower the melting temperature of the pore solid. The pore size distribution can be obtained by the nuclear magnetic resonance (NMR) method [10]. Due to the electromagnetism of the hydrogen atom, the water and hydrocarbon in tiny pores can be detected by an exerted magnetic field. The NMR method is widely used in oil fields during recent years, including the measurement of porosity, pore size distribution [36], wettability [37], and permeability [38] as well as the recognition of gas and oil [39,40].
All mentioned above make the process of heat transfer at low temperature complicated. In this paper, the temperature evolution of permafrost rock during the exploitation of gas hydrate reservoir is studied. A governing equation of heat transfer with phase transition is established to investigate the evolution of temperature. Based on the Gibbs-Thomson equation [41,42], the relation between ice content and temperature can be assessed by the pore size distribution curve. By simulating the heating test, the surface relaxivity in the NMR method is accurately calibrated, and the pore size distribution is obtained. We finally simulate the temperature disturbance with the model established in the exploitation, providing the basis of the theory for further studying of wellbore stability.

Theory
Considering that the infinite reservoir is centrally symmetric with respect to the wellbore, the analysis is regarded as one dimension. In order to simplify the analysis process, similar to hydrate, we take water as the analysis object. The heat conduction as the only form of energy exchange in rocks, the governing equation of heat transfer in a one-dimensional model is [43,44] where T is temperature, t is time, ρ s is the density of rock (kg/m 3 ), c is specific heat (J/kg·°C), λ is thermal conductivity (W/m·°C), and ∇ is gradient operator. With the variation of temperature, the phase transition occurs in the pores. Therefore, the governing equation of heat conduction with the phase transition is where L is the latent heat of fusion per unit mass of water (3:33 × 10 5 J/kg), θ i is the volumetric fraction of ice, and ρ i is the density of ice. In this equation, heat is transited in the form of specific heat that changes the temperature and latent heat that changes the phase. Assuming the freezing/thawing process is instantaneous, Equation (2) can be rewritten as where ρ s c − Lρ i dθ i /dT is termed equivalent volumetric heat capacity, and the term dθ i /dT should be determined first.
When freezing in the pores, the melting temperature is not constant but depends on the pore size. The phase transition in porous solids is governed by the Gibbs-Thomson equation: where r i is the smallest pore access radius invaded by ice crystals, ΔS v is the entropy of fusion per unit volume (1.2 MPa/°C for ice), γ iw is the water-ice interface stress (0.04 J/m 2 ), and α is the contact angle of the water-ice interface (α = 0 for thawing process). The Gibbs-Thomson equation reveals that it exists a critical radius of pores for a given temperature. At this temperature, ice crystal only forms in the pores bigger than the corresponding pore ( Figure 3). This is also the principle of measuring the pore size distribution through the relation between ice content and temperature.
To obtain the pore size distribution of rock, the NMR method is undoubtedly one of the most advanced methods [10]. Compared with mercury intrusion porosimetry (MIP), NMR will cause less damage to the sample, and the sample can be a repeated test. Most importantly, the NMR measurement can detect the hydrogen atom in the pores. In this case, 2 Geofluids water as the research object is measured directly in this paper. Therefore, the measurement results should be more reliable.
By monitoring the decay of the signal, the T 2 distribution can be obtained. T 2 is the transverse relaxivity time, which is a set of attenuation coefficients correlating with the pore size. The larger the pore, the higher the T 2 value and vice versa.
Thus, a pore size distribution curve can be calculated based on the T 2 distribution curve. The pore size can be represented by the T 2 value, and the Amp value, which reflects the volume fraction occupied by the corresponding size of pores. Indeed, due to the assumption of the analysis object, the amplitude value is the proportion of pore water. The T 2 distribution can be converted into the pore size distribution according to the equation [45][46][47]: where ρ is the surface relaxivity (μm/s), S is the surface area of the pore (m 2 ), and V is the pore volume (m 3 ). The specific area (S/V) equals to 3/r if the pore is assumed to be spherical and equals to 2/r if a cylindrical model is applied, where r is the radius of the pore (m). For different materials and different pore fluids, the surface relaxivity is different, and the pore distribution is sensitive to the surface relaxivity. As mentioned above, the pore size distribution is related to the ice content evolution with time, which directly decides the variation of the equivalent volumetric heat capacity, influencing the whole process of heat transfer with phase transition. Hence, the surface relaxivity can be calibrated from the experiment of heat transfer.   3 Geofluids

Experiments
Heating tests were conducted to calibrate the surface relaxivity. All samples were cut into a column with a height of 50 mm and a diameter of 25 mm (Figure 4(a)). For each sample, a hole with a diameter of 6 mm and a depth of 15 mm was drilled for the placement of a temperature sensor, which automatically monitored and recorded the sample temperature every 30 s. The two ends of the sample were left for water saturation; the lateral side of the sample was covered by phenolic resin to reduce water loss from the rock surface. After being dried at 80°C in an oven for a week and vacuumized for two hours, the samples were saturated in a container at 10 MPa for 24 h. Based on the difference between the dry sample and the water-saturated sample, the porosity was measured as 21.9%. After that, the NMR method (Figure 4(b)) is applied to obtain the T 2 distribution. The T 2 distribution is shown in Figure 4(c).
After the preparation, all samples first were cooled to -40°C for 4 hours in a refrigerator and then were taken out to the ambient temperature. The temperature was measured by the temperature sensor every 30 s, and the temperature evolution with time of the sandstone is shown in Figure 5. The   4 Geofluids temperature varies rapidly with time, and it reveals that the rock skeleton and ice crystals in the pore are heated. After that, the temperature increases slowly and features a plateau. Obviously, this temperature plateau is related to the endothermic effect during the phase transition of ice. The occurrence of this temperature plateau is an interesting phenomenon, and its features will be analyzed in the next section. In the last, the temperature increases rapidly again and finally will stabilize at the ambition temperature.

Temperature Evolution with Time.
To calibrate the surface relaxivity, we assume that the pores are cylindrical, and thus, the specific area S/V equals to 2/r. As indicated before, the pore size distribution can be obtained from the T 2 distribution according to Equation (5). It should be noted that the surface relaxivity of sandstone range from 1 to 46 μm/s for sandstones [48], and the surface relaxivity has a sensitive effect on the pore size distribution (shown in Figure 6(a)).
To determine the value of surface relaxivity, the maximum and minimum surface relaxivities (i.e., 1 μm/s and 46 μm/s) are set for discussion. The pore size distribution under the two conditions is obtained from Equation (5). With the same shape of the pore structure curve, the peak value of the pore size distribution is 7 μm for 1 μm/s and 323 μm for 46 μm/s. The cumulative curve of sandstone is also obtained according to the relation between the porosity and the amplitude value. In Figure 6, two curves indicate the "upper bound and lower bound" of the pore structure of rocks; thus, the optimized pore structure can be obtained via further analysis. As for the pore size distribution in two values, the melting point of ice corresponding to the different pores is given by the Gibbs-Thomson equation (Equation (4)). Meanwhile, the amplitude value is the signal of pore water, and thus, the ice volume at a given temperature is the sum of pore vol-umes larger than the given pore radius. The relation between ice and temperature is shown in Figure 7(a). There is a variation rate of ice content with temperature for two conditions, both near 0°C. When the surface relaxivity is 46 μm/s, the ice content evolution with temperature decreases steeply from 0.219 around 0°C, while there is a slow descent of ice content between 0 and -20°C with surface relaxivity equals to 1 μm/s. It indicates that with the increase of the surface relaxivity, the corresponding pore size region increases as well, and the ice content evolution with temperature becomes faster.
With the determination of ice content evolution, the equivalent volumetric heat capacity is obtained, and the temperature curves in two conditions are shown in Figure 7(b) according to Equation (4). For both temperature curves, the temperature gradually increases with time first and increases slowly with time at around 0°C (i.e., both feature a plateau regime). It is apparent that the occurrence of temperature plateau is related to the endothermic effect of the phase transition of ice. The rock structure and pore ice are heated first, and thus, the temperature rises rapidly. For the temperature plateau, from Equation (3), it can be demonstrated that the appearance of this plateau results from the large increase in the equivalent heat capacity, more exactly, the latent heat term. The value of the latent heat term is determined by the dθ i /dT, which has been obtained by the ice content evolution with temperature before. Hence, the dramatic increase in the latent term is related to the rapid decomposition of ice (Figure 7(a)), implying a significant phase transition.
With using the different surface relaxivity, two temperature curves are different in two aspects. One is that the heating rate in large surface relaxivity is higher than that in small surface relaxivity before the temperature plateau. This phenomenon reflects that the equivalent heat capacity of the heating curve corresponding to the large surface relaxation rate is smaller. According to the ice evolution with temperature, the slow descent of ice content will increase the 5 Geofluids equivalent heat capacity in the condition of 1 μm/s, and so the heating rate is slow. In addition to the difference in the heating rate, the temperature plateau of 46 μm/s is longer than that of 1 μm/s, which can be explained by the relation between ice content and temperature. For the condition of 46 μm/s, the ice content within the small temperature range (around 0°C) varies sharply with the temperature. For the condition of 1 μm/s, although its ice content decreases obviously at low temperature, the overall slope is less than that of 46 μm/s, leading to a small equivalent heat capacity. Therefore, it concludes that with a large surface relaxivity, the pore size distribution measured by NMR corresponds to the large equivalent heat capacity, and the longer the temperature plateau is reflected in the temperature curve.
As discussed previously, the real surface relaxivity is calibrated from the experimental curve of the heating test, which is 3.5 μm/s. The experiment results and the simulation results are shown in Figure 8. The simulated temperature curve has a high degree of coincidence with the experiment, but what is interesting is that the temperature plateau of the experimental curve is higher than the simulated curve (above 0°C). The error is inevitable, but we think it is since the time cost for temperature variation during the melting of ice is not considered. In other words, because of the heat transfer process and the anisotropy of the sample, the temperature of the rock would continue to rise a little during the phase transition.
4.2. The Disturbance in Temperature Field. During the exploitation of gas hydrate, the dissociation of the hydrate phase will absorb a lot of heat, leading to the variation of temperature and even the disturbance in the temperature field of the reservoir. Considering the exploitation at a low-temperature region (below 0°C), ice crystals are the pore solid that melts during   Geofluids the phase transition. With the governing equation of heat transfer with phase transition established, the pore size distribution is obtained from the above NMR data. Considering these features, the simulation model in 2D is constructed (Figure 9). The diameter of the wellbore is 0.2 m, and the reservoir layer is homogeneous with a radius of 10 m. Being in contact with the drilling fluid, the wellbore boundary is the convective boundary (the convection coefficient is set 500 W/(m·°C)), and the outer boundary is the free boundary.
Assuming the initial temperature of the permafrost layer is -4°C, and the temperature of the well fluid is 20°C. Finite element analysis is used to simulate the variation of the temperature in the reservoir within 30 days. The spatial distribution of the temperature field in the reservoir at different times is shown in Figure 10. It can be seen that the temperature field at different moments decreases exponentially with the increase of distance. At 1 hour, due to the phase transition near the wellbore, the heat transfer is very slow, and the reservoir temperature field barely changes. With the decompose of the solid phase, the heat is gradually transferred from the wellbore to the surrounding area, and the temperature is transferred as far as 5 m at 30 days. It reveals that the reservoir has a multistage process of phase transition during the exploitation of hydrate, and only when the phase transition is completed, the heat can transfer to a further place without thermal limitation. In fact, it should be noted that the decomposition of hydrate will produce gas, and the transmission of gas must accelerate the transfer of heat.

4.
3. Discussion. The researches above have shown that once the temperature is higher than the formation temperature, the hydrate decomposes quickly, and the temperature in the permafrost diffuses fastly. From the engineering aspect, the temperature of drilling fluids, well completion fluids, and well cementation exothermic could be over 353.15 K, which is much higher than the hydrate decomposition temperature. Thus, the dramatical temperature increment of hydrate reservoirs by operation of drilling, well completion, and well cementation provided the massively generated gas and water, which are very likely to cause overpressure and gas channeling in the pore, even leading to blowouts and other accidents. After that, the massive decomposition of hydrates and the migration of gas and water will cause the lack of framework and the decrease of the bearing capacity in the permafrost, which easily leads to engineering accidents such as collapse, sand production, and subsidence. Meanwhile, there is no trap on the permafrost hydrate reservoirs; the production operation may induce the fractures and fluid passageway, which conduct the massive decomposed methane uncontrollable leaking into the earth's surface, water, and atmosphere. Therefore, the disturbance of the temperature field in permafrost will greatly influence the production efficiency of hydrate exploitation, causing engineering problems and even engineering accidents, environmental risks, and geology disasters.
From the previous experience [12], the onshore production of hydrate in permafrost is the road to the offshore production test. Although some researchers made the in situ temperature measurements and numerical simulation of temperature behavior in offshore methane hydrate production [49,50], the preliminary analysis could not fully explain the temperature distribution behavior and its affection. Because most marine gas hydrates occur in the submarine landslide zones, so the reservoir framework and reservoir stability would be affected by the maximum temperature and temperature distribution around the wellbore during operation of drilling, well completion, and well cementation in offshore hydrate production. It is significant to control the maximum temperature and temperature distribution during engineering operation which is dramatically different from the conventional petroleum industry. Thus, the investigation of temperature distribution in the permafrost during hydrate exploitation is a necessary way to realize the efficient and safe production of gas hydrate, and it also provides theoretical support for the exploitation of marine gas hydrate.

Conclusion
This paper bases on the heat transfer in gas hydrate-bearing sediments during exploitation. After establishing a governing 7 Geofluids equation of heat transfer with phase transition, the relation between ice content and temperature can be obtained from the pore size distribution according to the Gibbs-Thomson equation. The pore size distribution is obtained by the NMR method, and the surface relaxivity is calibrated from the experimental data. Accordingly, the process of heat transfer with phase transition is analyzed, and the disturbance in the temperature field of the hydrate reservoir is investigated. The main conclusions are as follows: (1) A governing equation is established describing the heat transfer in rocks at low temperature, in which both specific heat and latent heat terms are combined into an equivalent heat capacity. The evolution of ice content with temperature can be obtained from the pore size distribution. The NMR method is available to get the pore size distribution, and the surface relaxivity can be calibrated from the heating test (2) The temperature of samples remains quasiconstant for some moment. This is mainly related to the intense phase transition that absorbs a great quantity of latent heat. Temperature plateau is determined by the relation between ice content and temperature, which is assessed by the pore size distribution curve calculated from NMR data. Thus, the pore size distribution can be obtained by measuring the temperature evolution curve of material saturated water at low temperature. It also demonstrates that the position and length of this temperature plateau can be predicted by the pore size distribution data (3) The evolution of the temperature field in the reservoir during the exploitation is investigated according to the model of heat transfer established. The result shows that the heat transfer is remarkably limited by the endothermic effect of phase transition. Recovery efficiency may depend largely on the heating capacity of the drilling fluid and the rate of gas production (4) Compared to conventional petroleum engineering, it is key to control the maximum temperature and temperature distribution of onshore permafrost hydrate reservoirs, which would give theoretical support and guide for the controlling temperature behavior of marine gas hydrate development

Data Availability
Our data is from our experiments.

Conflicts of Interest
There is no conflict of interest regarding the publication of this paper.