In this study, a one-dimensional analytical model to describe heat and mass transfer during methane hydrate dissociation under thermal stimulation in porous media has been developed. The model is based on a similarity solution that considers a moving dissociation boundary which separates t he dissociated zone containing produced gas and water from t he un -dissociated zone containing only methane hydrate . The results of temperature distribution, pressure distribution, energy efficiency, and parametric study considering various initial and boun dary conditions as well as various reservoir properties are presented and compared with previous studies . Sensitivity analysis of gas production on reservoir properties is also presented in this paper. The dissociation boundary moves faster by increasing the heat source temperature while decreasing the heat source pressure simultaneously , but the associated e nergy efficiency decreases. Increasing the well thickness has a negative effect on the energy efficiency of the process. Among the proposed thermal properties of the system, only the thermal diffusivities and c onductivites of the reservoir as well as the p orosity of the sediment affect the dissociation. The main contribution of this work is investigating analytically the hydrate dissociation using thermal stimulation by taking into account the effect of wellbore thickness and structure. 1 Introduction Gas hyd rates are cage -like substances (clathrates ) in which a large amount of gas molecules is trapped by crystalline water molecules. They contain about 170 to 180 standard cubic meters of natural gas per cubic meter of hydrate (Kamath and Godbole, 1987) . The high pressure and low [Type here] temperature underground with sufficient methane gas induce methane hydra te (MH) formation (Davie and Buffett, 2001) . The estimation of the volume of the MH exceeds the whole volume of conventional gas resources worldwide (Englezos, 1993; Makogon, 1981) . However, considering MH as a reliable and alternative source of energy for future depends on the availability of hydrates and the cost of gas prod uction process fro m them. The MH are mainly formed beneath the permafrost in arctic regions and in marine shelf sediments according to appropriate pressure and temperature conditions (Collett, 2001; Collett, 2008) . As mentioned above, MH formation and stability require specific pressure and temperature conditions, so any change in temperature or pressure from equilibrium condition s would induce MH instability. The three main methods of MH dissociation that have been reported in previous works are as follows (Li et al., 2016; Pooladi -Darvish, 2004) : i) thermal stimulation, in which the temperature of the reservoir is increased above the equilibrium temperature; ii) depressurization by decreasing the pressure inside the MH reservoir lower than the equilibrium pressure of the hydrate, so that the MH will start dissociating; and iii) depressurization in conjunction with thermal stimulation , that can improve the efficiency of the dissociation process (Wang et al., 2014) . It should be noted that gas production aspects and potential using different methods is not yet fully investigated. The mathematical bac kgrounds in the field of MH dissociation can be divided into two main categories: analytical and numerical. Analytical solutions are fast and simple, also they can provide better mechanistic view . On the other hand, the numerical studie s are more comprehen sive with smaller number of assumptions. In 1982, one of the first numerical stud ies was developed for MH dissociation under depressurization using a three dimensional (3 D) model consisting of an MH layer and a free gas zone , in which conduction heat trans fer and gas fluid flow were considered (Holder and Angert, 1982) . This work was extended in 1986 by considering water movement, which is produced from MH dissociation (Burshears et al., 1986) . In neither of the mentioned studies , the w ater/gas convective heat transfer was considered. In 1991, a numerical model was designed to predict MH dissociation in porous media by depressurization considering three phases of water, gas, and MH with gas -water flows and neglecting the effect of heat t ransfer (Yousif et al., 1991) . Then, it was extended by considering water -gas flows and convective -conductive heat trans fer (Masuda, 1997; MASUDA, 1999) . In 2000, MH dissociation upon thermal stimulation [Type here] was simulated by assuming a moving dissociation boundary, which separates the dissociated zone from the undissociated zone, and considering different media permeability (Tsypkin, 2000) . This moving boundary layer was assumed to be impermeable, so if the pressure and temperature conditions fall below the balance point the dissociation will be halted. Anot her numerical work employing finite difference method was performed to study gas produ ction form MH dissociation using depressurization (Ahmadi et al., 2004) . They considered the effect of heat transfer and reported that the process is a function of well pressure. At the same year, gas production from MHs using both depressurization and thermal stimulation was shown to be possible by using TOUGH2 simulator, which is capable of simulation by cons idering four components and up to nine phases in either kinetic or equilibrium model s and upon different dissociation processes (Moridis, 2002) . The feasibility of gas production from MH reservoirs in in the Mack enzie Delta, Northwest Territories, Canada was also studied using TOUGH2 simulator (Moridis et al., 2004) . The results revealed the gas production using depressurization and thermal stimulation, and higher efficiency achieved by using both the methods together. They additionally showed that increasing the initial formation temperature, the well temperature, and the formation thermal conductivity increase d the amount of produced gas; while, it is not affected by the perme ability of the formation and the specific heat of the rock and hydrate. The numeric al simulations in another work showed that the kinetic reaction models should be taken into account in order to avoid under -prediction of recoverable MH, however, it require s more computational effort compared to equilibrium reaction model (Kowalsky and Moridis, 2007) . Investigation on d isperse oceanic MH reservoirs with low hydrate saturation using TOUGH -Fx/HYDRATE si mulator upon depressurization showed low gas production potential, which was not economically feasible, with high amount of water production (Moridis and Sloan, 2007) . Nowadays, researchers are performing more mathematical investigations on this field accompanying with comparison with real field or experim ental results or using parameters from real reservoirs (Chen and Hartman, 2018; Konno et al., 2017; Mardani et al., 2018; Wang et al., 2019; Wang et al., 2016a; Zhao et al., 201 6). In 1990 , Selim et al. (Selim and Sloan, 1990) designed an analytical one-dimensional (1D) study on MH dissociation upon thermal stimulation method by considering a moving dissociation boundary . They took into account the gas convection heat transfer and its flow, w hile the water by-product remain ed motionless in the sediments . In 1982, a study on hydrate dissociation process upon hot water injection and depressurization employing two models has been designed: the [Type here] frontal -sweep model , and the fracture -flow model . The y reported that the gas production from hydrates using the depressurization method is feasible compared to thermal stimulation (McGuire, 1982) . Makogon (Makogon, 1997) reported analytical expressions for the temperature and pressure profiles for the MH dissociation upon depressurization using a similar process of moving dissociation boundary including the effect of the throttling process in the energy equation. Tsypkin (Tsypkin, 2001) extended it by considering similarity solutions for temperature and pressure distribut ions and the water and gas movement. Another model was also prese nted i n 2001 by including heat conduction to Makogon’s model (Ji et al., 2001) . An analytica l work on dissociation induced by depressurization reported that the effect of the gas -water two -phase flow on MH dissociation is smaller than the effect of heat transfer and the intrinsic kinetics of MH decomposition (Hong et al., 2003) . Recently, an analytical -experimental work has been reported to study MH decomposition by depressurization, thermal stimulation, and depressurization in conjunction with thermal stimulation (Wang et al., 2015) . The analytical model in that work was based on performed experi mental conditions, for example, heat transfer from outside the experimental reservoir into the reservoir is also considered . In addition to mathematical approaches, experimental tools have also been employed to investiga te MH dissociation upon different m ethods. One of the major challenges in experimental works in this field, is the size of the setup that significantly affects the outcome of tests (Wang et al., 2016b) . For example, scale of experiment setup determines the main mechanism involving in the hydrate dissociation in porous media (Wang et al., 2016b) , which could be one of the followings: i) the intrinsic kinetics of hydrate decomposition; ii) heat transfer in the decomposing zone; or iii) the multiphase flow (i.e., gas-water flow) during gas production (Hong et al., 2003) . Tang et al. (Tang et al., 2007) reported that the intrinsic kinetics of hydrate decomposition is the determining factor in the core -scale experiments. They a lso showed that in larger scale works, such as field works, heat transfer in the decomposing zone becomes the controlling mechanism. An experimental investigation of the MH dissociation in a 3D cubic hydrate simulator (CH S) upon huff and puff thermal stimu lation showed that MH dissociation process upon thermal stimulation proceeds as a moving boundary ablation process (Li et al., 2011b) . Li et al. (Li et al., 2012b) experimentally investigated the MH dissociation upon depressurization using two hydrate simulators with differen t scales. They revealed that gas production period is longer for larger scale simulators. In another experimental work, conduction heat transfer was reported as the main mechanism for heat [Type here] transfer to the dissociating zone (Zhao et al., 2012) . Wang et al. (Wang et al., 2016b) studied the hydrate dissociation behavior below the quadruple point in the s andy sediment employing a 3D Pilot -Scale Hydrate Simulator (PHS). They revealed that MH diss ociation below the quadruple point induces ice formation in pores, which in turn increases the dissociation rate. Another work, designed and performed by Wang et al . (Wang et al., 2 018), reported that depressurization in conjunction with thermal stimulation is the optimum method to dissociate water -saturated hydrate samples using a pilot -scale hydrate simulator via the three abovementioned main methods. Studies over the past decades have provided important information on MH Dissociation and its gas production potential. Apart from the interactions happening in the reservoir upon hydrate dissociation, the wellbore structure may affect the dissociation especially the heat transfer mech anism in thermal stimulation method. Furthermore, data from several studies suggest that the wellbores in the production region mainly consisted of three layers: casing, cement, and gravel (Florez Anaya and Osorio, 2014; Pucknell and Mason, 1992; Xu et al., 2014) . On the other hand , among previous mathematical works that have been done in the field of hydrate dissociation up to now (Li et al., 2010; Li et al., 2012a; Wan et a l., 2018; Zhao et al., 2016) , no analytical study has investigated the impact of heat source structure (i.e., wellbore radius and the associated outer layers ) on dissociation upon thermal stimulation using wellbore heating , which might induce unreliable outcomes while comparing to those from experiments or field works. Existing research addresses this gap and recogniz es the critical role played by the wellbore structure in hydrate dissociation via wellbore heating. In the present work , an analytical soluti on is generated to model a semi -infinite 1D hydrate reservoir in Cartesian coordinates using a flat heat source considering three main layers: casing, cement, and gravel. The model and results are validated and compared with the previous mathematical and e xperimental studies . The r esults of this wor k will be useful in assess ing the gas production from MH reservoirs using thermal stimulation. The importance and originality of this study are that it analytically explores, for the first time, the effect of wel lbore structure on the MH dissociation and t he associated heat transfer process . Therefore, t he findings should make an important contribution to the field of gas production from MH reservoirs by providing the outcomes obtained using conditions that are closer to the real-condition tests , which make them more valuable and reliable . [Type here] 2 Modeling Figure 1 show s a schematic of hydrate dissociation using thermal stimulation method in the proposed 1D flat case. It should be mentioned that the reservoir is semi -infinite, and L denotes the overall well thickness including, cement, gravel , and casing layers . The fundamental steps of MH dissociation are as follows: in the beginning of the process, the reservoir with the porosity of  is filled wi th MH in equilibrium condition s with initial temperature T0 . At time t = 0, the temperature of the heat source at x = 0, increases to a new temperature of Ti , which is above the equilibrium temperature of the hydrate, and is kept constant during the proc ess. Subsequently, MH starts dissociating and a sharp moving boundary surface is created as shown by dashed lines in Figure 1, which indicates the rate of hydrate dissociation and separates the water and gas produced in the dissociated zone (Zone I) from the un-dissociated zone (Zone II) . The t emperature of the dissociated zone is higher than the temperature in the hydrate zone inducing heat conduction from the dissociated zone to the hydrate zone and hydrate dissociation. The rate of hydrate dissociation is determined by the speed of moving interface, which decrease s as the process continues and the dissociated zone thickens . Actually, it consumes a larger part of input heat to increase the temperature of the matrix material in the dissociated zone as well as the produced water and gas at the dissociation front as the dissociation interface moves forward . The remaining input energy is consumed for hydrate dissociation and temperature increase of the matrix materials in Zone II. As shown in Figure 2 , the temp erature and pressure distributions in the Zone I ( TI and P ) are respectively s I iT T T and isP P P , and in the Zone II, the temperature distribution is 0 II s T T T and the pressure is in the equilibrium condition of the hydrate. The dissociation interface is shown by the dashed line and the associated pressure and temperature change due to time. It should be noted that TI at the outer surface of the well changes with time , but it is always lower than Ti . [Type here] Figure 1. Schematic of hydrate dissociation in the semi -infinite flat reservoir . The dissociation interface is identified by the dashed line, and the grey region show s well thickness. Figure 2. Schematic of pressure and temperature distribution in the reservoir upon hydrate dissociation. According to the mass conservation law and the sudden change in density at the dissociation front due to gas production , the produce d gas will be st reaming towards the heat source. The assumptions considered in this study, which were in accordance with the previous analytical works (Selim and Sloan, 1990; Tsimpanogiannis and Lichtner, 2007) , are as follows: i) hydrates filled the entire pores of the media; ii) temperature and pressure at the dissociation interface are at the thermodynamic equilibr ium; iii) the pr oduced water from the dissociation remains motionless in the pores; iv) thermophysical properties of the phases are constant; v) ideal behavior of gas is applied in the equation of state of gas; vi) the produced gas is in thermal equilibriu m with the local Zone I P(x,t) TI(x,t) Zone II TII(x,t) TII → T 0 Ps(t) Ts(t) x (m) 0 L ∞ S(t) Pi [Type here] sediments (with the same temperature); vii) viscous dissipation and inertia effects are neglected. It should be noted that the t hermal properties of the system are assumed to remain constant during dissociation The fundamental set of equations and the sol ution procedu re for the flat case by considering a wellbore heat source are presented by the following expression: The continuity equation of gas in Zone I is: 0g g g v tx   +=       , ,0 L x S t   (1) where, is porosity, gis gas density (kg/m3), and gv is gas velocity (m/s). The gas velocity in Zone I is calculated by the multiphase Darcy’s law, which is as follows: gkPvx=− , ,0 L x S t   (2) where k is relative permeability of gas (md),  is gas viscosity (mPa.s), and P is gas pressure (Pa). Equations 3 and 4 show the energy balance in Zones I and II, respectively: 2 1 2g pg g I II I pIC v T TTCkt x x+=   , ,0 L x S t   (3) 2 2II II IITT tx= , ,0 S x t   (4) where I is density (kg/m3) of zone I, g is gas density (kg/m3), pIC is specific heat cap acity (J/(kg.K)) of Zone I, pgC is specific heat capacity of gas (J/(kg.K)), IT is the temperature (K) in Zone I, IIT is the temperature (K) in Zone II, II is therm al diffusivity (m2/s) of Zone II, and Ik is the thermal conductivity (W/( m.K)) of Zone I. In equation 3, the gas temperature in each location in Zone I is assumed to be equal to the local sediment temperature. This assumption has been proved by Weinbaum and Wheeler (Weinbaum and Wheeler Jr, 1949) and also employed in another work by Selim and Sloan (Selim and Sloan, 1990) . The gas behavior is assumed t o obey the ideal gas law. So, the ga s density in Zone I can be evaluated by the following equation: [Type here] g ImP RT=, ,0 L x S t   (5) where m is gas molecular mass (kg/mol), and R is the universal gas constant (J/(mol.K)). The above set of equation s represents the fundamental concept of the process. In the following, the associated initial and boundary conditions are stated. Temperature and pressure at the inner wall of the wel l (x = 0) are kept constant and respectively equal to Ti and Pi. It is assumed that there is no pressure drop in the well wall so the pressure at the outer surface of the well (x=L) is Pi as well. The assumed wellbore schematic in the production region is shown in Figure 3, in which the thickness of each of the casings is 0.7 cm, and the thickness for the cement and gravel respectively is 2.5 cm and 1.5 cm . Hence, t he heat transfer equation in the well bore could be stated by equation 6: ()I i I Iw wT T TkAxR−−= , ,0 x L t= (6) 1 g sc w w s c gL LLRA k k k = + + (7) Figure 3. Schematic of the proposed wellbore structure. Different parts are as follows: 1) casing, 2) gravel pack, 3) ca sing, and 4) cement where Aw is the wellbore area (m2), wRis thermal resistivity of the well bore (W/K) , xkand xL are respectively the thermal conductivity (J/(s.m.K)) and thickness (m) with s, c, and g subsc ripts respectively stands for steel (casing), cement, and gravel . It should be noted that this proposed [Type here] structure for wellbore is a general geometry and model, and there are many of the wellbore structures and the associated geometries in the literature th at are different from this model. The pressure at the dissociation interface is calculated from the Antoine equation , which is a thermodynamic relationship with the interface temperature: ( / )ea a sA B T sP−= , ,0 x S t= (8) where sP is the pressure (Pa) and temperature (K) at the moving interface, and aA and aB are constants. The mass and energy balances at the dissociation interface are respectively represe nted in equations 9 and 10. Also, t he heat of MH dissociation is calculated by equation 1 1 (Selim and Sloan, 1990) , and equations 1 2-14 represent the boundary conditions : ( ) 0gH H g gdSFvdt+= , ,0 x S t= (9) II I II I H HdTT dSk k Qx x dt−= , ,0 x S t= (10) hyd sQ c dT=+ (11) ()I II sT T T t== , ,0 x S t= (12) 0 IITT= , ,0 xt→  (13) 0 IITT= , 0 , 0xt  = , SL= (14) where H is the hydrate density (kg/m3), IIk is the thermal conductivity (W/( m.K)) of Zone II, HdQ is hea t of MH dissociation (J/kg ), and c and d are constants . FgH in equation 8 is a constant, representing the ratio of mass of the methane gas trapped inside the MH to the mass of hydrate. In this work, the value of 0.1265 kg CH 4/kg hydrate is chosen for FgH, which has been used in a previous study by Selim et al. (Selim and Sloan, 1990) . The following equations 1 5-17 are obtained respectively from equations 1, 3, and 8 by employing equations 2 and 5 in order to eliminate the gas velocity a nd density. 0 IIP k P P t T x T x     −=         (15) [Type here] 2 1 2pg II I pIkmC TT PC P kt R x x x +=     , ,0 L x S t   (16) ( ) 0gH H IdS kmP PFdt RT x−= , ,0 x S t= (17) To solve the above set of equations, the similarity solution introduced by N eumann (Ã–zisik et al., 1993; Carslaw and Jaege r, 1959) is employed. In this solution, t he movement of the dissociation interface is assumed to be proportional w ith the square root of time ( t1/2). This assumption satisfies the initial and boundary conditions, and the following dimensionless transforma tion, shown in equation 17, is employed to simplify and solve the abovementioned equations. 4IIx t = (18) So, on the moving dis sociation front, equation 1 8 becomes: 4IIS t = , ,0 x S t= (19) Also, on the outer surface of the well: 4os IIL t = , ,0 x L t= (20) The transformation of the above equations using equation 17 is developed and presented in the supplementary information . By considering the procedure suggested by Özışık et al. (Ã–zisik et al., 1993) , Carslaw et al.(Carslaw and Jaeger, 1959) , and Selim et al. (Selim and Sloan, 1990) , error function is employed to choose the following solutions for the temperature distributions: 0( ) ( ) ()I IIT Aerf a b Aerf b B T T Cerfc = + − + =+ (21) A, B, a, and b constants are introduced in the supplementary information. The pressure distribution (Pa) in Zone I can be calculated from equation S 7 as follows: ()1/2 2 04 osgH H II IFRP P T dkm    =+ (22) [Type here] That by replacing and integrating TI from equation 2 1, becomes: ( )1/22 0 1 1 ( ) ( )( ) ( )( ( ) ( ))os os P P J K AJ L L       = + − + − (23) where ()J , ()K , and 1()L are defined in the supplementary information. Heat flux at the boundary (J/(s.m2)) as a function of time is calculated from the following formula: I fITukx=− , xL= (24) The transf ormed version of equation 2 4 considering equation 18 is provided in the supplementary information. Cumulative he at input into the reservoir from the heat source (J/m2) is calculated by integrating equation 23 with regard to time as follows : 0tI ft ITQ k dtx=− , xL= (25) In order to calculate the t otal volume at standard temperat ure and pressu re (STP) of dry gas, equation 2 6 is employed , where fpV is the total volume of produced gas per surface area of the moving interface up to time t (m3/m2), ftnis the total moles of gas produced per s urface area of the moving interface up to time t (mole/m2), and STPT and STPP are temp erature and pressure of gas at STP conditions , respectively . ft STP fp STPn RTVP= (26) More details of total v olume calculation are provided in the supplementary information. For the efficiency evaluation of gas production during the dissociation process using thermal stimulation method, the energy e fficiency ratio is introduced as follows (Song et al., 2015) : fp g f ftVQ Q= (27) where f is the energy efficiency ratio in this case , and gQis the heating value of the gas at STP conditions (J/m3). Actually, f represents the amount of energy that would be produ ced from [Type here] combustion of the produced gas to the amount of input energy to the system from the heat source. It becomes apparent that an energy efficiency higher than one is preferre d. It should be mentioned that the s urface area of the well that is in contac t with the reservoir is equal to 1m2. The obtained solution for temperature and pressure distributions can satisfy the abovementioned basic equations and boundary conditions (equa tions 1 -17) by direct substitution. MATLAB programming software is used for a ll calculations in this work. 3 Results and discussion Equation 1 8 indicates that the location and movement of the dissociation interface are dependent on  . It should be noted that  represents both the dim ensionless inter face position and dimensionless interface velocity ( 1/2/ (4 )s IIvt ). Thus, t he higher the value of  , the faster the movement of the dissociation interface. As shown in the equation S23, the value of  is dependent on sPand sT values, whi ch are calculated using the pressure and temperature at the heat source wall (out er boundary) according to the equations in the previous section. The temperature at the outer boundary of the heat source is time dependent , inducing sP and sT at the dissociation front to change over time . Figure 4 shows the value of  considering spec ific physical and boundary conditions presented in Table 1, which are obtained and set from the previous studies (Cheng et al., 2011; Dalla Santa et al., 2017; Re mund, 1999; Selim and Sloan, 1990) .  increases from a slightly lower value at the beginning but converges to a value as the temperature at the w all surface converges to temper ature Ti of the inside of the well , which is kept constant . Table 1. Parameters u sed in the modeling . Total w ellbore thickness, L, m 0.054 Cement thickness, Lc, m 0.025 Gravel thickness, Lg, m 0.015 Casing thickness, Ls, m 0.014 Contact area of well bore with reservoir in production region , m2 1 Thermal conductivity of cement , kc, W/(m.K) 0.933 Thermal conductivity of gravel , kg, W/(m.K) 0.4 Thermal conductivity of casing (steel) , ks, W/(m.K) 43.3 Porosity,  0.3 Permeability, k, μm2 1 [Type here] Thermal diffusivity of Zone I, I , μm2/s 2.89×106 Thermal conductivity of Zone I, Ik , W/(m.K) 5.57 Thermal diffusivity of Zone II, II , μm2/s 6.97×105 Thermal conductivity of Zone II, IIk , W/(m.K) 2.73 Hydrate density, H , kg/m3 913 Heat of dissociation of hydrate, HdQ , J/kg 3446.12 10 132.638sT − Gas heat capacity, pgC , J/(kg.K) 8766 Gas viscosity,  , Pa.s 10-4 Heating value of the gas at STP conidtions, gQ, MJ/m3 37.6 Molecular mass of methane, m, g/mol 16.04 Mass ratio of the methane gas trapped inside the hydrate to the mass of hydrate, FgH 0.1265 Universal gas constant, R, J/(mol.K) 8.314 Figure 4. Dimensionless position of interface during 1000 days dissociation process at Ti=563.5 K, Pi=7.6 MPa , and T0=275K . The effect of boundary conditi ons on the interface movement is presented in Figure 5, which shows the value of  with regards to different well pressure and temperature after 100 days of dissociation . The results shown in Figure 5 have the same trend as the res ults reported in the 11.11.21.31.41.51.61.71.81.9 0 200 400 600 800 1000β Time (day) [Type here] previous work by Selim et al. (Selim and Sloan, 1990) , but with a slightly lower values of  (about 0.1 difference). This difference is due to the conduction heat transfer in the wel lbore t hickness, which causes the amount of transferred heat to the dissociation surface becomes lower compared to the previous work (Selim and Sloan, 1990) , in which a flat heat source with constant temperature without wellbore thickness is considered . In all cases, decreasing the heat source pressure while increasing its temperature induces higher values of  . Figure 6 represents the relation between the interfa ce velocity and MH temperature, T0. The higher the hydra te temperature, the higher the interface movement . These results are consistent with data obtained in a Numerical work by Liang et al. (Liang et al., 2010) , who also found that decreasing the pressure and increasing the surrounding temperature would increase the rate of hydrate dissociation upon depressurization. They also validated their results against experiments performed by Masuda (MAS UDA, 1999) . Figure 5. Dimensionless position of interface at T0=280 K and various Ti and Pi values. 1.41.51.61.71.81.922.12.22.32.4 1 10 100β Pi (MPa)Ti=640K Ti=620K Ti=600K Ti=580K Ti=560K Ti=540K Ts=280 K Ti=640K Ti=600K Ti=560K Ti=620K Ti=580K Ti=540K [Type here] Figure 6. Dimensionless position of interface at Ti=563.5 K and various T0 and Pi values . Furthermore, in F igures 5 and 6, the associated temperature and pressure of the locus for which the dissociation temperature is equal to the temperature of MH (280 K) is shown . The dissociation temperature is mainly dependent on heat source pressure , while the effect of the te mperature of the heat source and MH tempera tures on Ts is almost negligible ; higher Pi increases Ts, while lower Pi decreases Ts. This is also in agreement with the previous work by Selim et al. (Selim and Sloan, 1990) , as they showed the direct dependency of Ts on Pi and mentioned that for lower pressures of 6 MPa, Ts may reduce to the freezing temperature of water which causes the interruption of the dissociation process due to ice generation. It should be mentioned that when Ts is eq ual to T0, all the transmitted heat to the interface is consumed for the dissociation and no heat is transferred to or from the hydrate zone , and the temperature of the hy drate zone remains constant. O n the other hand, as Ts falls below T0, some part of th e required heat for dissociation is absorbed from the hydrate zone , resulting in the temperature reduction of the hydrate zone . Temperature and pressure distributions are calculated and presented in Figures 7 and 8 for the following two boundary condition s (BCs): BC 1) 450iT= K, 10iP= MPa, and 0280 T= K, and 1.41.51.61.71.81.922.12.22.3 1 10 100β Pi (MPa)T0=270 K T0=272 K T0=274 K T0=276 K T0=278 K T0=280 K Ts=280 K T0=270K T0=274K T0=278K T0=272K T0=276K T0=280K [Type here] BC 2) 563.5iT= K, 7.6iP= MPa, and 0275 T= K considering three different time frames . The temperature distribution in Zones I and II are separated by a horizontal dashed line that shows the interface temperature . However, d ue to the small temperature increment at the well surface in longer time frames (as shown in the Figure 7 ), the interface temperature changes slightly . The temperature change at the well surface also causes the interface pressure to change in longer time frames . It should be noted that as the dissociation progresses and the temperature at the well surface converges to that in side the well, the temperature and pressure at the interface respectively tend to decrease and increase gradually and converge to the associate d values reported Selim et al. ’s work (Selim and Sloan, 1990) , in which the associated value s for the interface temperature and pressure for the two proposed boundary conditions were reported to be: BC 1) Ts = 285 K, Ps = 290 Pa; BC 2) Ts 282.5 K, Ps = 560 Pa . Tsimp anogiannis et al. (Tsimpanogiannis and Lichtner, 2007) in their semi -analytical model showed that increasing the temperature of well would increase the pressure at the interface, which is also shown in our results of boundary conditions of BC 1 and BC 2. The results of temperature distribution also illustrate that the distance between the dissoci ation surface and the heat source is higher in the previous study without considering the wellbore thickness (Selim and Sloan, 1990) compared to the present model due to the negative effect of heat conduction in the well, but this difference decreases in longer time frames as the dissociated zone thickens and the effect of heat source thickness disappears gradually . The results presented in Figure 7 for temperature distribution are also in good agreement with those re ported by Li et al. (Li et al., 2011b) , who experimentally investigated MH dissociation upon thermal stimulation. They also reported that the decomposition progressed by a moving boundary, which separates dissociate zone form the undissoicated zone. [Type here] Figure 7. Temperature distribution for case 1 at differe nt time frames for two initial and boundary conditions of a) BC 1 b) BC 2. Figure 8. Pressure distribution for case 1 in the dissociated zone at different time farmes for two initial and boundary conditions of a) BC 1 b) BC 2. The volume of produced gas in STP conditions , amount of input heat, and energy efficiency upon hydrate dissociation considering the two initial and boundary conditions during 100 days of gas dissociation are shown in Figure 9 . The volume of pr oduced gas and the amount of input heat are higher considering BC 2 compared to BC 1 . Howeve r, the energy efficiency is higher considering BC 2 compared to BC 1 due to the higher difference between the amount s of input heat and relatively smaller differenc e between the associated produced gas as shown in Figure 9. The results of energy efficiency indicate that increasing the heat source temperature and decreasing its 270290310330350370390410430450 0 1 2 3 4Temperature (K) Distance (m)t=10^4 s t=10^5 s t=10^6 st=104s 270320370420470520570 0 2 4Temperature (K) Distance (m)t=10^4 s t=10^5 s t=10^6 s 050100150200 0 1 2 3P-P0(Pa) Distance (m)t=10^4 s t=10^5 s t=10^6 s 0100200300400 0 1 2 3P-P0(Pa) Distance (m)t=10^4 s t=10^5 s t=10^6 sa) b) a) b) t=105s t=106s t=105s t=104s t=106s t=105s t=104s t=106s t=105s t=106s t=104s [Type here] pressure would increase the rate of dissociation, but will not increase the total efficienc y of the process. It is somewhat surprising that the energy efficiency increases to a peak point in the beginning t hen it decreases as dissociation progresses. This could be due to the constant surface area of the dissociation interface , which is equal to the surface area of the heat source , during the dissociatio n that restrict s the amount of gas production during dissociation . On the other hand, the amount of input heat increases such that it cannot reach the dissociation front due to the sediment matrix of Zone I and dissociation products . Therefore, the speed of input heat increment exceeds the speed of the produced gas increment until the temperature at the heat source surface gradually converges to that inside the heat source, and the energy efficiency also con verges to a constant value . The same trend for energy efficiency (i.e., significant increase at the beginning then decreasing smoothly) was also reported by Li et al. (Li et al., 2014; Li et al., 2011a) and Wang et al. (Wang et al., 2013) , who investigated the MH dissociation behavior up on thermal stimulation using cubic hydrate simulator s. They claimed that the sudden increase of the energy e fficiency at the beginning is due to the rapid hydrate dissociation close to the wellbore at the beginning of the process. In another experimental w ork performed by Wang et al. (Wang et al., 2014) on the gas production form MH using thermal stimulation by hot water, the same tre nd for the energy efficiency and gas production was reported . It was claimed that the pre -depressurization b efore heat stimulation caused the sudden increase of the energy efficiency at the beginning, then it decreased and converged to lower value. They re ported an energy efficiency of around 4 after 100 minutes of dissociation , which shows a good agreement with our results . In the work of Selim et al. (Selim and Sloan, 1990) an energy efficiency between 6.4 -11.2 was reported , which is close to the results presented in the current paper and of high value because of the similar conditions considered in the two works. The energy efficiency reported in this work is consistent with that of Tang et al. (Tang et al., 2005) , who experimentally investigated the gas production from MH reservoirs upon thermal stimulation. They reported an energy ratio, which was calculated in the same way as the energy efficiency in the pr esent work, of about 5 upon conditions close to the BC 2. They also revealed that increasing the initial temperature and decreasing the init ial pressure improved the energy ratio. Bayles et al. (Bayles et al., 1986) studied analyticall y MH dissociat ion upon cyclic steam injection. They reported the same trend for the energy efficiency for one year dissociation, which was converged between 4 to 9.6. Again, our results are in good agreement with [Type here] theirs, and the slight difference could be due to the direct steam injection into the reservoir and the cyclic pattern of the process. Figure 9. a) volume of produced gas , b) amount of input heat, and c) energy efficiency during hydrate dissociation for two BCs. A parametric study is performed in order to investigate the effect of various parameters, stated in Table 2, on the dissociation process. The obtained results from the parametric study , shown in Figure 10, indicate the following points: a) the higher the thermal diffusivity and conductivity of Zone I, the higher the dissociation movement ; b) the lower the thermal dif fusivity and conductivity of Zone II, the higher the dissociation movement; and c) the perme ability and gas viscosity almost have no effect on the dissociation process as the interface velocity does not cha nge by various values of each. The same results we re reported by Selim et al. (Selim and Sloan, 1990) . Zhao et 050100150200250300350400450500 0 20 40 60 80 100Produced gas (m3) Time (day)BC 1 BC 2 05001000150020002500300035004000 0 50 100Input heat (MJ/s.m2) Time (day)BC 1 BC 2 024681012 0 50 100η Time (day)BC 1 BC 2a) b) c) [Type here] al. (Zhao et al., 2015) investigated the numerical model of gas production form MH using thermal stimulation and showed that increasing the thermal conductivity promote d the dissociation process , on the other hand, different relative water -gas permeabilit ies have almost no effect on the gas production . They also numerically showed that increasing sediments’ thermal conductivity would increase the rate of gas generation at the beginning of hydrate di ssociation upon depressurization (Zhao et al., 2014) . It should be me ntioned that both of their works were validated by the experimental work of Masuda (MASUDA, 1999) . Tsimpanogiannis et al. (Tsimpanogiannis and Lichtner, 2007) investigated a parametric study on the effect of different parameters on the MH dissociation upon thermal stimul ation. They revealed that increasing the thermal conductivity of the porous media induced more MH dissociation. Table 2. Range of parameters employed in the parametric study. Porosity,  0.1 to 0.5 Permeability, k, μm2 0.1 to 5 Thermal diffusivity of Zone I, I , μm2/s 1×106 to 5×106 Thermal conductivity of Zone I, Ik , W/(m.K) 3 to 7 Thermal diffusivity of Zone II, II , μm2/s 4×105 to 8×105 Thermal conduct ivity of Zone II, IIk , W/(m.K) 1 to 5 Gas viscosity,  , Pa.s 10-4 to 10-6 [Type here] Figure 10. The effect of various values of the parameters stated in Table 2 on the interface movement after 100 days dissociation. a) thermal diffusivity and thermal conductivity of Zone I, b) thermal diffusivity and thermal conductivity of Zone II, c) various values of porosity with different permeabilities, and d) various values of por osity with different gas viscosities. 4 Conclusion s In this study an analytical model to investigate MH dissociation using thermal stimulation is presented. It takes into account the effect of wellbore thickness which may include multiple layers with differ ent thermal properties. E valuat ions of the effect s of various reservoir parameters and boundary conditions on the dissociation process have been performed . The two main factors that have been considered for assessing the process are: a) moving dissociation boundary 11.52 3 4 5 6 7β KIalpha I=10^(-6) alpha I=2(10^(-6)) alpha I=3(10^(-6)) alpha I=4(10^(-6)) alpha I=5(10^(-6)) 1.31.82.3 1 2 3 4 5β KIIalpha II=4(10^(-7)) alpha II=5(10^(-7)) alpha II=6(10^(-7)) alpha II=7(10^(-7)) alpha II=8(10^(-7)) 11.21.41.61.822.2 0.1 0.2 0.3 0.4 0.5β ϕ11.21.41.61.822.2 0.1 0.2 0.3 0.4 0.5β ϕa) b) d) c) [Type here] position/velocity ( λ), and b) energy efficiency of the process ( η). These findings have significant implications for the understanding of how the wellbore structure affect the hydrate dissociation by providing an analytical solution for thermal s timulation method by taking into account heat conduction through the wellbore , which makes the models closer to the real conditions. This study appears to b e the first study that examines associations between the wellbore structure and hydrate dissociation . In the limit when wellbore thickness approaches zero the model is in good agreement with previous similar analytical study. Additionally , the results of t he present study match with those of previous mathematical and experimental investigations. Based on the results of this investigation , the following main conclusions can be made : - The rate of dissociation and efficiency of reservoir heating are dependent on both well thickness and composition (which are design parameters). Well thickness causes a reductio n in the dissociation front speed and produced gas amount compared to the cases without considering the well thickness. - Temperature of the well boundary is not constant , which is also dependent on well thickness and composition. - As a result, p ressure and temperature at the dissociation front are not constant during dissociation . - Increasing the well temperature , but decreasing its pressure , induce higher dissociation interface velocity and volume of produced gas, but reduces energy efficiency. Thus, in order to increase the energy efficiency of the process, a detailed study on the in itial and boundary conditions including well composition (number of layers, their thickness and thermal properties) of the system should be performed. - The good agreement between the results of this study and previous experimental and numerical studies provides validation of the assumption s made in the model development . Acknowledgement Financial support for this work provided by Natural Sciences and Engineering Research Counci l of Canada (NSERC) [Type here] Nomenclature A Dimensionless constant in equation S9 ftn Total moles of produced gas per surface area of the moving interface up to time t Aa Dimensionless constant in equation 8 P Pressure in Zone I Aw Wellbore area Ps Pressure at the interface a Dimensionless constant in equation S12 Pi Pressure of the heat source B Dimensionless constant in equation S10 STPP Pressure of gas at STP conditions Ba Dimensionless constant in equation 8 HdQ Heat of MH dissociation BC Boundary conditions gQ Heating value of the gas at STP conditions b Dimensionless constant in equation S13 ftQ Total input heat to the reservoir from the heat source C Dimensionless const ant in equation S11 R Universal gas constant pIC Specific heat capacity of Zone I Rw Thermal resistivity of the wellbore pgC Specific heat cap acity of gas c Dimensionless constant in equation 1 1 S Interface position D Dimensionless constant in equation S14 Ti Temperature of the heat source d Dimen sionless constant in equation 11 IT Temperature in Zone I E Dimension less constant in equation S15 IIT Temperature in Zone II F Dimensionless constant in equation S16 T0 Initial temperature of hydrate FgH Ratio of mass of the methane gas trapped inside the MH to the mass of hydrate Ts Temperature at the interface ()G Dimensionless constant in equation S20 STPT Temperature of gas at STP conditions ()H Dimensionless constant in equation S20 t Time ()I Dimensi onless constant in equation S20 uf Heat flux from the well ()J Function in equation S17 fV Volume of produced gas per surface area of the moving interface in the time fraction of “t,t-1” k Permeability fpV Total volume of produced gas per surface area of the moving interface up to time t Ik Thermal conductivity of Zone I vs Interface velocity IIk Thermal conductivity of Zone II X Axial position ck Thermal conductivity of cement g Gas density gk Thermal conductivity of gravel H Hydrate density sk Thermal conducti vity of casing  Gas viscosity ()K Function in equation S18  Dimen sionless variable in equation 18 [Type here] L Thickness of the well  Dimen sionless constant in equation 19 1()L Function in equation S19 os Dimen sionless variable in equation 20 cL Cement thickness  Porosity gL Gravel thickness I Thermal diffusivity of Zone I sL Casing thickness II Thermal diffusivity of Zone II MH Methane hydrate gv Gas velocity m Gas molecu lar mass f energy efficiency ratio fn Total moles of produced gas per surface area of the moving interface in the time fraction of “ t,t-1” 5