Modeling of reservoir temperature upon preheating in SAGD wells considering phase change of bitumen

International Journal of Heat and Mass Transfer 144 (2019) 118650

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Modeling of reservoir temperature upon preheating in SAGD wells considering phase change of bitumen Yanfang Gao, Mian Chen ⇑, Botao Lin, Yan Jin State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China

a r t i c l e

i n f o

Article history: Received 5 April 2018 Received in revised form 17 August 2019 Accepted 27 August 2019

Keywords: Oil sand Bitumen Temperature Preheating Steam circulation Phase change

a b s t r a c t Steam circulation prior to production is usually employed to establish inter-well thermal communication and create an expected initial steam chamber in the process of bitumen recovery using steam assisted gravity drainage (SAGD). To realistically evaluate the efficiency of steam circulation, and to find the indicator for terminating the steam circulation and switching to SAGD production, modeling of the temperature propagation considering the bitumen-rich oil sand reservoir properties and oilfield operations is required. This study proposed two new modified methodologies, including an analytical model as well as a numerical model, in which the latent heat effect caused by the phase change of bitumen was involved for thermo-phase change coupled analyses. The modeling was conducted for a real SAGD steam circulation project in an oilfield located in Karamay, northwest of China, using oilfield operation parameters and thermo-dynamical properties obtained from laboratory tests. The models can predict the temperature propagation and bitumen phase change behavior in the domain encompassing dual wells. In general, the predicted results obtained from both analytical and numerical approaches agree well with those in published articles in terms of the cooling time derived from temperature falloff tests and inter-well midpoint temperature. The modeling draws major conclusions that the temperature drops rapidly with the distance to wellbore in the fusion zone, but decreases slowly in the solid zone. The heat flux on wellbore exhibits an exponential reduction in the beginning, but ultimately tends to be a constant. The moving speed of phase change interface decreases with injection time, and finally increases up to only a short distance. The inter-well midpoint temperature is heated to exceed the threshold temperature, and the whole inter-well area is ultimately occupied by the fusion zone. These findings not only can provide a criterion for evaluating the efficiency of steam circulation and an indicator for terminating the steam circulation and switching to production, but also the approaches proposed may be used to analyze these heat injection and phase change problems occurring in other types of oil and gas-rich reservoirs like the natural gas hydrate formation. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction China is very rich in its ultra-heavy oil (also known as oil sand bitumen) resources. The largest ultra-heavy oil deposit is located in the northwestern part of the Junggar Basin in Xinjiang province, and possesses about 360 million tons of oil sand resources. These immobile heavy oils that have dramatic viscosities as high as 106–107 mPa s in initial reservoir conditions, are not recoverable in its natural state by conventional oil well production methods including currently used enhanced recovery techniques. Aware of this issue, the steam assisted gravity drainage (SAGD) technique has been extensively implemented in the development of Xinjiang

⇑ Corresponding author. E-mail address: [email protected] (M. Chen). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118650 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

ultra-heavy oils (bitumen) since 2012. In general, the application of SAGD consists of two stages including startup stage and production stage. The end purpose of the startup stage is just to form a desired uniform stable steam chamber and to find when and how to switch to production, which are both obstacles of SAGD technique. In the Karamay Fengcheng oilfield (the largest ultra-heavy oilfield in China), conventional startup only contains the preheating stage, during which both horizontal injector (I well) and producer (P well) in a SAGD well pair are put under steam circulation to establish inter-well communication and to create an initial steam chamber, which impacts the efficiency of bitumen recovery tremendously [1]. In the past decade, encouraged by the success in Alberta oil sands where geomechanics is proactively utilized for the in-situ oil sand recovery, dilation technique under water injection (also known as micro fracturing) has been gradually

2

Y. Gao et al. / International Journal of Heat and Mass Transfer 144 (2019) 118650

employed in the Fengcheng oilfield for shortening the subsequent cyclic preheating time [2–4], so now it is also an essential part of startup. The prediction of reservoir temperature during the startup phase is one of the most important issues for field engineers, because it can evaluate the efficiency of steam circulation and estimate the optimal opportunity of switching to SAGD production from circulation. However, it is not easy to make a good prediction due to the unique properties of Karamay oil sand bitumen. As shown in Fig. 1, the previously immobile bitumen on preheating gradually becomes significantly less viscous and undergoes a change from a solid to fluid phase. In this regard, this heat transfer problem should consider the considerable latent heat effect caused by phase change of bitumen containing a large proportion of solid waxes in the Karamay oil sand reservoir, especially in the domain encompassing the dual wells [5,6]. There is no exaggeration that this is a challenging task. Up to now, many researchers have studied the heat transfer problem in the process of SAGD. To sum up, there are two main heat transfer mechanisms in SAGD process [7–9]. Heat conduction is usually considered as the dominant mechanism [10–12], and only a few people held that the convection caused by mobile condensate flows may not be ignored [13–16]. Nevertheless, it is widely acknowledged that heat conduction is the only heat transfer mechanism in the preheating stage, because the steam can hardly infiltrate into the reservoir and therefore no mobile condensate flow exists before switching to production under a small adequate differential pressure between I well and P well [17–19]. At present, the analytical solution and numerical simulation are two frequently-used methods for the prediction of temperature propagation in the reservoir near a well pair while preheating. For the first method, published analytical models are only suitable either for the production stage where the steam chamber already exists and the considerable heat convection induced by condensate flows occurs [20–23], or for the preheating stage where the heat conduction is only considered incorporating these factors like the reservoir heterogeneity and startup strategies [17,18]. For the second method, the numerical simulation by CMG STARS upon preheating can be regarded as a preliminary procedure prior to production, and there is no substantive difference between two procedures [24–27]. In this paper, the two methods above will not be directly used in solving the heat transfer problem in the Karamay oil sand reservoir for their shortcomings. Instead, a new analytical approach and a finite element method (FEM) by ABAQUS will be employed in this field for its logical, accurate and non-empirical calculations [28]. Recently, some people have conducted preliminary studies on the heat transfer during the preheating stage, but they can’t solve the problem considering latent heat effect of bitumen and the size effect of wellbore. The aforementioned analytical and numerical models concentrated on heat conduction problems applicable in

the reservoirs that contain heavy oils with low viscosity (<10000 mPa s) and low wax-content (<10%), as long as the phase change of heavy oils occurs with no or little latent heat. In these models, the analyses of temperature propagation in the reservoir around a horizontal well pair upon preheating were studied, given that reservoir thermo-physical properties are constant and the boundary condition is either constant heat flux or constant temperature [1,29]. Wellbores in these models were simplified as straight lines for some analytical solutions, and some empirical equations were used to eliminate the size effect of wellbore. Afterwards, some models were modified for the investigations of optimum startup strategies, operational procedures and heat requirements in the preheating period [18,19]. Recently, Zhu et al. proposed a modified model in which the reservoir is divided into a hot zone and a cold zone where their thermophysical parameters differ, so as to the hot zone size and shape can be estimated during steam circulation [17]. In general, all these models above are some modified forms derived from basic equations of heat conduction in the ultra-heavy oil reservoir [7,30,31]. However, it is unfortunate that latent heat effects occurring in bitumen reservoirs have not be mentioned. All in all, no approaches are available in literatures to quantitatively describe the heat transfer behavior in oil sand reservoirs considering bitumen’s phase change upon preheating. The situation hinders the design of a proper cyclic steam injection scheme before its implementation. Therefore, a new evaluation method is desired for the SAGD engineers, and the aim and objectives of this study are just to propose a mathematical description for this complex problem and further solve it by both analytical and numerical approaches for the estimation of temperature propagation in the oil sand reservoir considering latent heat effect of bitumen and size effect of wellbore. 2. Materials The ultra-heavy oil reservoirs in the Qigu and Badaowan formations were deposited during the Jurassic period. The major reservoir areas are located at the Fengcheng oilfield, which is 130 km to the northwest of Karamay city, Xinjiang province, northwest China. The formations are formed of oil sand grains filled with bitumen and mud, overlying an overlapping pinch-out belt on the upper side of the Karamay-Wuerhe overthrust belt in the northwestern margin of the Junggar basin [3]. The in-situ viscosity of Karamay bitumen exceeds 5
106 mPa s under initial reservoir conditions. These bituminous unconsolidated sandstones named oil sands are unique engineering geomaterials for two reasons: firstly, the bitumen filling pores is essentially a solid under virgin conditions; and secondly, the sands themselves are not loosely packed beach sands [32]. Instead, they have terrestrial unconsolidated, interlocked structures like but not as dense as Alberta oil sands, which are formed as a result of deeper burial and elevated

Fig. 1. Microstructures of Karamay oil sands (a) under water injection, and (b) after preheating [4].

3

Y. Gao et al. / International Journal of Heat and Mass Transfer 144 (2019) 118650

temperature over a long geological period [33]. Tables 1 and 2 exhibit these essential thermophysical parameters of the oil sands and bitumen collected from the depth of 303–304 m in the Qigu formation of the Fengcheng oilfield. Most of parameters were determined by laboratory tests and some supplementary data was derived from literatures [2,6,30,34–37]. Considering that these experimental studies (especially experiments under elevated temperature) on the thermodynamics of Karamay oil sands are very limited, some test data of Alberta oil sands is also given for reference. As shown in Table 1, the oil sand, regarded as a kind of unconsolidated sandstone, has a lower density than other rock formations such as the tight sandstone, chalk, limestone and shale. Besides predominant oil sand grains and mud components, there are considerable bitumen and little water in the oil sands. The oil sand has a lower thermal conductivity than conventional sandstones, which is unfavorable for the preheating efficiency. But the oil sand owns a lower specific heat capacity, which benefits for the whole SAGD process. It is convinced that the thermal properties of oil sand samples used are found associated with temperature [30]. As shown in Table 1, the oil sand under a higher temperature has a lower thermal property, probably because of the fusion of bitumen. As shown in Table 2, the bitumen under initial reservoir temperature contains 56.7% saturates and aromatics, which may keep solids as a form of wax before preheating. Once the wax in bitumen is heated up to a critical temperature of 43–56 °C, it will fuse and absorb heat of 266 J/g on the phase change interface [34], which is usually called the latent heat effect.

Table 2 Properties of Karamay and Alberta bitumen. Properties

Karamay

Alberta

Hydrocarbon group mass composition w (%) Saturates 42.0 Aromatics 14.7 Resins 37.2 Asphaltenes 6.1

22.0 21.0 39.0 18.0

Viscosity l (mPas) 20 (°C) 50 (°C) 100 (°C) 150 (°C) 200 (°C)

642,000 9070 130 32 12

5,060,000 61,200 846 92 22

Table 3 Model parameters for the analytical solution. Model parameters

Data

Bottom hole vapor temperature Tw (°C) Initial temperature Ti (°C) Equivalent phase change temperature Tc (°C) Thermal diffusion coefficient of the solid zone as (
106 m2/s) Thermal diffusion coefficient of the fusion zone af (
106 m2/s) Thermal conductivity for the solid zone ks (W/m°C) Thermal conductivity for the fusion zone kf (W/m°C) Wellbore radius rw (m) Reservoir density qf (
103 kg/m3) Mass content of the solid wax fwax (%) Equivalent phase change latent heat DH (J/g) Total circulation time ttotal (days) Analysis step length Dt (day)

252 12 56 0.55 0.40 1.1 0.75 0.125 1.96 4.71 266 300 1

3. Methods 3.1. Physical model The SAGD preheating stage where the steam circulates through dual horizontal wells before production, is designed to establish inter-well hydraulic and thermal communications to create a steady steam chamber. Upon preheating, the previously immobile bitumen gradually becomes significantly less viscous and undergoes a change from a solid to fluid phase due to thermal stimulation. In the preheating process, it is acknowledged to divide the whole reservoir into different zones, characterized by these properties in petrophysics, geomechanics and thermodynamics [17,38]. With respect to the formation geometry, the studied area is divided into a rheological fusion zone and a solid zone, whose thermodynamical properties like thermal conductivity, thermal diffusivity and volumetric heat capacity, differ dramatically. As

exhibited in Fig. 2, it is reasonable to first consider a simple case where there is only a single well in the oil sand reservoir, so as to uncover some desirable methods and results with regard to this complicated problem. In Fig. 2, a new 2D physical model on the heat transfer contemplating the phase change of bitumen in an oil sand reservoir around a single well was depicted, neglecting the heat transfer interference in the direction along borehole axis. The studied reservoir, between mudstone cap rock and base rock, can be reasonably regarded as an infinite reservoir, because the heat transfer can be hardly affected by the mudstones over the whole steam circulation according to field monitoring data. Meanwhile, the steam can hardly infiltrate into the reservoir and therefore reservoir pressure

Table 1 Oil sand properties in Karamay and Alberta. Properties

Karamay

Alberta

Physical properties Density q (
103 kg/m3) Porosity u (%) Absolute permeability ka (mD)

1.96 23–37 346–5801

1.75–2.40 16–46 500–1700

Mass composition w (%) Bitumen Water Solid

8.3 0.6 91.5

10.4 3.5 86.3

Thermo-dynamical properties at 20 °C Thermal conductivity k (W/m °C) Thermal diffusivity a (
106 m2/s) Volumetric heat capacity qc (
106 J/m3 °C)

0.931 0.535 1.73

1.1 0.55 2.1

Thermo-dynamical properties at 250 °C Thermal conductivity k (W/m°C) Thermal diffusivity a (
106 m2/s) Volumetric heat capacity qc (
106 J/m3 °C)

– – –

0.75 0.4 2.0

Fig. 2. Schematic geometry for heat transfer problem in an oil sand reservoir around a single well.

4

Y. Gao et al. / International Journal of Heat and Mass Transfer 144 (2019) 118650

changes can be ignored, which means that heat convection hardly occurs. Based on all descriptions above, a heat conduction problem considering the phase change behavior of bitumen in an infinite oil sand reservoir was established in a cylindrical coordinate. However, in fact both I well and P well play a significant role upon preheating (Fig. 3). The high-pressure steam engines on the ground that can generate high pressure high temperature steam, and advanced completion strings that can permit the steam to be circulated through the dual wells, are specially designed for successful and efficient preheating. The working pressure is controlled to be higher than the initial reservoir pore pressure but no more than the fracturing pressure for security reasons. It is noteworthy that the SAGD preheating process is quite different from the steam flooding. On one hand, both I well and P well in the preheating stage are all intended to circulate steam instead of production in P well meanwhile injection in I well. On the other hand, there are both inlet strings (long tubes) and outlet strings (short tubes) in both I well and P well. The steam will be injected from the inlet strings, releasing its latent heat on wellbores, then the condensed water will be returned up to the ground by outlet strings. Little steam is forced into the reservoir filled with bitumen in the circulation stage [1,17,39], therefore it is reasonable to consider this process as a heat conduction problem with a fixed boundary of a time-dependent heat source on wellbores and a moving phase change boundary where the latent heat effect occurs. As used in many literatures [1,17–19,29], this coupled issue can be overcome by a 2D problem. In this model, the initial reservoir is divided into three zones, i.e., the pay zone, cap rock, and base rock. In general, the oil sand reservoir is usually treated as a porous media, which is saturated with water under water injection and saturated with oil-water mixtures after preheating. Here, the pay zone is regarded as a kind of phase-change composite material based on the theory of mixtures. In this paper, the phase change behavior of bitumen consisting of plenty of solid waxes will be quantified by these parameters including the fusion temperature and phase change latent heat. It was found in the oil sand reservoir, usually characterized by a higher porosity and permeability than common sandstones, that this immobile, solid bitumen does exist under initial reservoir conditions before preheating [4,6,33,36,37]. There are many field cases in which wax has been reported in a natural state in the reservoir [5]. In this regard, a consideration of latent heat effect in bitumen reservoirs is reasonable and essential. Future work will focus on some thermo-dynamical researches for the auxiliary

determination of these parameters used in analytical and numerical models below. In this paper, the numerical solution was mostly focused on, and the analytical solution can be found in these literatures [40–45]. 3.2. Numerical solution scheme To solve the problem above, a modified solution of temperature propagations in the reservoir was proposed, derived from a classical model that consists of critical governing equations about the heat conduction with a moving phase change boundary, which was initially applied in the metallurgical industry. First, a basic model for the prediction of the reservoir temperature around a single horizontal well considering the phase change of bitumen was established. The solution of the simple model above with a timedependent temperature boundary was obtained using a modified Paterson’s method. Second, the temperature field around the dual SAGD wells was calculated by superposition, given that the steam can hardly infiltrate into the reservoir and pore pressure changes can be ignored upon preheating, which means the phase change temperature of bitumen is constant. The basic model with only a single well was first considered, and the specific process was described as follows. The unsteady heat conduction models have been frequently adopted in the simulation of steam circulation in ultra-heavy oil sand reservoirs [1,17–19,29]. The temperature of the fusion zone at any time t and any position M(r) can be expressed as

@ 2 T f 1 @T f 1 @T f þ ¼ ; r @r @r 2 af @t

r w 6 r 6 r sðt Þ

ð1Þ

where t is the circulation time; r is the studied position in the cylindrical coordinate; Tf is the temperature at t and r; af is the thermal diffusivity of the fusion zone; rw is the borehole radius; rs(t) is the phase change interface. The initial conditions and boundary conditions in the fusion zone can be expressed as T f ðr; t Þjr ¼ rw ;t¼0 ¼ T i , T f ðr; t Þjr¼rw ;t>0 ¼ T w and T f ðr; t Þjr¼rsðtÞ ;t>0 ¼ T c . Here Ti is the initial temperature of oil sand reservoir; Tw is the wellbore temperature; Tc is the phase change temperature. Similarly, the temperature of the solid zone at any time and any position can be expressed as

@ 2 T s 1 @T s 1 @T s þ ¼ ; r P r sðt Þ r @r @r 2 as @t

ð2Þ

where Ts is the temperature at t and r; as is the thermal diffusivity of the solid zone. The initial conditions and boundary conditions in the solid zone can be expressed as T s ðr; t ÞjrPrw ;t¼0 ¼ T i , T s ðr; t Þjr¼rsðtÞ ;t>0 ¼ T c and T s ðr; tÞjr!1;t>0 ¼ T i . The energy balance on the phase change interface can be expressed by

kf

@T f @T s dsðtÞ j  ks ¼ f wax qf DH j j dt r¼rsðtÞ ;t>0 @r r¼rsðtÞ ;t>0 @r r¼rsðtÞ ;t>0

ð3Þ

where kf and ks are the thermal conductivities of the fusion zone and solid zone, respectively; qf is the density of the fusion zone; fwax is the mass content of the solid wax; DH is the latent heat caused by the phase change of bitumen. An exact solution can be obtainable if the solution of heat conduction equations is chosen as an exponential integral function [31,45]. Therefore, this problem can be indirectly expressed by a transcendental equation as follows.

Fig. 3. Schematic geometry for heat transfer problem in an oil sand reservoir around a well pair.


 q Tc  Ti k2 af exp k2  ks  2  exp  k a 4p as Ei  f as

! ¼ k2 f wax af qf DH

ð4Þ

5

Y. Gao et al. / International Journal of Heat and Mass Transfer 144 (2019) 118650

where q is the heat flux per unit length through wellbores; k is a dimensionless parameter; Ei (x) is the power integral function. k
pffiffiffiffiffiffiffi and and Ei (x) can be expressed as k ¼ sðt Þ= 2 af t R 1 u Ei ðxÞ ¼ x e =udu, respectively. Supposing the wellbore radius is small enough compared with the outer boundary of studied reservoir, the wellbores can be reasonably viewed as straight linear heat sources as described in Paterson’s method, where q is the strength of a linear heat source. However, the studied wellbore radius in this problem can’t be ignored because of a closed reservoir boundary like cap and base rocks, so here q is the heat flux per unit length through wellbores at r = rw. In a normal steam injection or circulation scenario, it is a common sense that the well is usually operated under a fixed bottom borehole pressure or constant wellbore temperature if the steam is saturated in wellbores, which means the heat influx will not be constant but decline continuously [29]. In this regard, q can be expressed as a function of time t, which means the heat influx per unit length through wellbores will change with time if the wellbore temperature is constant.

q ¼ qðt Þ

ð5Þ

It is worth noting that the boundary condition T f ðr; t Þjr¼rw ;t>0 ¼ T w must be satisfied at any t. Therefore, the bound-

effects, although its computational efficiency may be lower than the analytical approach. This part describes the basic energy balance, constitutive models, boundary conditions, finite element discretization, and time integration procedures [28]. The basic energy balance is [46]

Z

V


 qðt i Þ Tc  Ti k2 af exp k2  ks  2  exp  k af 4p as E  i

T f ðr w ; ti Þ ¼ T c þ

¼ k2 f wax af qf DH

ð6Þ

as

 


 q r2  Ei k2 ¼ T w Ei  w 4pkf 4af ti

ð7Þ

Similarly, k is also time-dependent. In this regard, there will be only two undetermined coefficients including q and k if t is given. In this study, a new idea was proposed that q can be viewed as a local constant in a small timestep near any time t if t is given. Therefore, the temperature distribution at a given time ti can be obtained.

 


 q r2  Ei k2 T f ðr; t i Þ ¼ T c þ Ei  4pkf 4af t i T s ðr; t i Þ ¼ T i 

  Ti  Tc r2   Ei  4as t i Ei k2 aasf

ð8Þ

cðhÞ ¼

Z

qdS þ S

ð11Þ

rdV V

dU dh

ð12Þ

Heat conduction is assumed to be governed by the Fourier law.

f ¼ k

@h @x

ð13Þ

where k is the conductivity matrix, k = k(h); f is the heat flux; and x is position. The conductivity k can be fully anisotropic, orthotropic or isotropic. Boundary conditions can be usually specified as prescribed temperature h = h(x,t) and prescribed surface heat flux q = q(x,t) per area. The Galerkin approach was used to give the spatial discretization by following equation.

Z

_ þ NN qUdV

V

Z V

@NN @h dV ¼ k @x @x

Z

Z NN rdV þ

V

NN qdS; Sq

N ¼ 1; 2;   

ð14Þ

The time integration form is shown as follows.

h R 1

i R N M NN q dU j NM dV þ V @N  kjtþDt  @N dV c M dh tþDt @x @x R R R R ¼ V NN rdV þ Sq NN qdS  D1t V NN qðU tþDt  U t ÞdV  V Dt

ð9Þ

Z

where V is a volume of solid material, with surface area S; q is the density of the material; U_ is the material time rate of the internal energy, which is the latent heat effect in this study; q is the heat flux per unit area of the body, flowing into the body; and r is the heat supplied internally into the body per unit volume. This relationship shown as Eq. (11) written in terms of a specific heat is usually used, except for latent heat effects caused by phase changes, which are given separately in terms of solidus and liquidus temperatures (the lower and upper temperature bounds of the phase change range) and the total internal energy associated with the phase change, called the latent heat. When the latent heat is given, it is assumed to be in addition to the specific heat effect (Fig. 4).

ary condition above and Eq. (4) can be combined to obtain q and k at a given time ti

!

_ ¼ qUdV

V

@NN @x

@h  k  @x dV



with hNtþDt;iþ1 ¼ hNtþDt;i þ c M ; i ¼ iteration number

When a well pair is considered, the superposition principle is used to calculate the temperature increase at M(C) induced by both injector and producer, respectively, and then the total temperature under the circulation can be obtained. When the injection time is t, this process can be expressed as

where Dt is the time analysis step length. More detained information can be found in the ABAQUS theory guide [28].

T t ¼ T i þ DT I;t þ DT P;t

4. Model verification

ð10Þ

where Tt is the total temperature at the position M(C); DTI,t, DTI,t are the temperature increases induced by I well and P well, respectively. The one-dimensionless cylindrical coordinate should also be transformed to the two-dimensional Cartesian coordinate for a convenient description and an intuitive exhibition. Therefore, the r at any point M(r) in the reservoir can be written as r = (x2 + y2)1/2. The numerical method can not only validate the analytical results, but also deal with this strong nonlinear problem involving the complicated geometry, time-dependent initial/boundary conditions as well as temperature-induced thermal parameters. The ABAQUS capability for heat transfer analysis can be intended to model solid body heat conduction with temperature-dependent conductivity/diffusivity and internal energy including latent heat

ð15Þ

It is impossible to directly measure the midpoint temperature in inter-well area by underground measurement experiments. At present, temperature falloff test has been extensively employed for an in-situ heating process to determine the heating characteristics along the horizontal wellbores [17,29]. The temperature falloff test is normally conducted after steam has been injected or circulated into the wellbore for a long period of time, and its focused object is the wellbore instead of inter-well midpoint. Duong [29] discussed the phenomena of extreme high calculated cooling time and calculated thermal diffusivity when he conducted full length assessment using his radial conductive heating model for a SAGD well pair. He attributed these phenomena to the development of steam chamber in the start-up stage. Zhu et al. [17] used

6

Y. Gao et al. / International Journal of Heat and Mass Transfer 144 (2019) 118650

Fig. 4. Schematic curve for (a) specific heat, and (b) latent heat definition.

transient temperature analysis to evaluate steam circulation in the SAGD startup process, and conducted an analysis for the temperature falloff after steam chamber forms, so as to find out whether there is a good match between his model and Duong’s. In this study, the temperature falloff simulation by ABAQUS was conducted for model verification. After steam circulation for 300 days, the I and P wells were shut in for 1000 days. Fig. 5 shows the reservoir temperature distributions after shut-in for 1000 days. During the shut-in period, the temperature falloff curves with the shut-in time were observed for analyses. Fig. 6 is the temperature evolution at the inter-well midpoint and wellbore with shutin time. Both the curves in Fig. 6 can reflect and determine the reservoir heating characteristics, but only the wellbore temperature falloff curve is usually adopted, because the temperatures can only be recorded along the horizontal wellbores by placing the temperature recorders inside the wellbores such as thermocouples or fiber-optic sensors [29]. Moreover, the temperature falloff curve with extreme slope is difficult for analyses, so the semi-log coordinate was used to obtain an approximate linear relation, which is used to determine the zero-intercept and cooling time proposed in Duong’s paper (Fig. 7). The dimensionless time t* is equal to (Tw  T)/(Tw  Ti), where T is the reservoir temperature. It is revealed that the model used in this study has a good fit with those proposed by Duong [29] and Zhu et al. [17]. The cooling time in the studied area is about 700–800 days. The solid bitumen/fusing oil interface is a

Fig. 6. Temperature falloff curves of the midpoint and wellbore with shut-in time.

Fig. 7. Semi-log plot for temperature falloff tests conducted on the wellbores.

Fig. 5. Reservoir temperature distributions after shut-in for 1000 days.

reasonable interpretation of the heat scenario at the position corresponding to shut-in for 100 days. At present, most models describing the preheating process are concerned about the thermal convection caused by steam infiltrations into the porous reservoir and about the subsequent condensate/steam interface analyses [17,29]. In their models, the circulation stage involves the period in which terminating the steam circulation and switching to SAGD production are

7

Y. Gao et al. / International Journal of Heat and Mass Transfer 144 (2019) 118650

conducted. Therefore, the differential pressure between two wells was usually considered in their models, while the convection is not permissible in this study. In this regard, it is impossible to directly compare these predicted results obtained by same approaches, let alone same circulation conditions and model parameters like the circulation time, working pressure, working differential pressure, borehole radius, bottom hole temperature and reservoir thermo-physical parameters. Therefore, the predicted tendency obtained by different models and physical interpretations induced by different approaches are supposed to be focused on for analyses. For instance, the predicted temperatures in literatures may be slightly higher than the modeling solutions in this study, but these temperature increases with injection time are all analogical in tendency and only differ in the slopes determined by reservoir thermo-physical parameters. For another example, the predicted temperatures in this study after being subtracted by those in literatures can be interpreted by this fact that there are both heat conduction as well as convection in the circulation process when it is operated after switching to SAGD production. In this paper, the end purpose of these mathematical descriptions in the SAGD preheating stage is to quantitatively predict the temperature distributions that are impossible to be obtained by field measurements except for a little data monitored by the sensors attached to wellbores and to provide a reliable guidance for engineers in field practices. However, there are some unavoidable differences between these predicted models and field experiments, because of the existence of a lot of essential model simplifications and hypotheses. Therefore, searching out these differences and interpreting these substantive reasons are supposed to be focused on. These factors that induce a dissatisfied modeling effect will be involved in the model in further study. For example, the densities of the fusion zone and solid zone can be different to consider the local thermal convection [12,31]. 5. Results and discussions The temperature propagation upon preheating in the oil sand reservoir around a SAGD well pair is one of the most important issues for the oilfield engineers. A good uniform propagation can not only reduce the viscosity of bitumen down to an expected viscosity at a critical temperature that makes the heavy oil self-flow under reservoir conditions (this critical temperature is normally called as pour point or fusing point in literatures [5,6]), but also benefit the generation of the initial steam chamber that impacts subsequent mature steam chamber in production stage. In this section, the temporal and spatial distributions of temperature in some typical reservoirs will be studied.

respectively. In the Fengcheng oilfield, the circulation period always maintains for about 6–12 month because of the super high viscosity of bitumen under initial reservoir conditions [2]. All key parameters used for analytical solutions are listed in Table 3. The influence of cap rock and bed rock on the heat transfer in the studied reservoir can be evaluated by the numerical simulation, and some supplement data about the reservoir shape and the thermal parameters of the cap/bed rock is shown in Table 4. In this regard, the temperature propagation in the reservoir around a single well can be obtained. First, k evolutions with injection time t are given, considering Eqs. (8) and (9) by both analytical and numerical methods. The analytical solution is slightly smaller than the numerical simulation result, but this can be neglected. The dimensionless parameter k decreasing with time t shows a sharp drop tendency in the beginning but tends to be a stable value around 0.37. The phase change interface location can be obtained if the k is given, and the analytical and numerical solutions are exhibited in Fig. 8. First, the analytical solution is very coincident with the numerical solution. Second, the interface moving speed decreases with injection time, but shows different increase rate in the whole injection process. For instance, it takes no more than 100 days for phase change interface to reach 1.5 m with the distance to wellbore axis, but finally increases up to just 2.4 m at the end of steam circulation. In a normal steam injection or circulation scenario, the well is usually operated under a fixed bottom borehole pressure or

Table 4 Model parameters for the numerical simulation. Model parameters

Data

Cap rock thickness Hc (m) Reservoir thickness HR (m) Base rock thickness Hb (m) Reservoir width W (m) Well position with the distance to top surface of base rock Hw (m) Thermal conductivity of cap rock kc (W/m °C) Thermal conductivity of bed rock kb (W/m °C) Specific heat of cap rock cc (J/gK) Specific heat of base rock cb (J/gK) Cap rock density qc (
103 kg/m3) Bed rock density qb (
103 kg/m3) Solidus temperature Ts (oC) Liquidus temperature Tl (oC)

2.5 30 2.5 30 10 2.75 2.20 0.765 0.714 2.63 2.63 55 57

5.1. Case study of I well of well pair A-1 Well pair A-1 corresponds to a recently studied reservoir in area A of the Fengcheng oilfield. The depth of the reservoir, which lies between a mud cap rock and base rock with a thickness of 2.5 m, ranges from 347 m to 377 m. The vertical depths of I and P wells are 367 m and 372 m, respectively. In this part, only I well was studied, so it is 10 m away from the top surface of base rock. The wellbore size is 0.25 m in diameter. The horizontal section of the well pair was designed to be 500 m in length. Tested results revealed that the initial reservoir temperature is 12 °C, and the bottom hole vapor temperature is kept at 252 °C. In-situ viscosity can exceed 5
106 mPa s under reservoir conditions, and the bitumen is essentially a solid under virgin conditions. Considering the phase change behavior of bitumen filling pores, the reservoir is divided into the fusion zone and solid zone, whose thermo-dynamical properties corresponds to those of 20 °C and 250 °C in Table 1,

Fig. 8. Analytical and numerical solutions of heat flux and phase change interface positions.

8

Y. Gao et al. / International Journal of Heat and Mass Transfer 144 (2019) 118650

constant wellbore temperature if steam is saturated in the wellbore. In this case, the heat influx is not constant, but declines continuously [18,29]. Therefore, the heat influx on wellbore is focused on to observe its tendency, and the results of both analytical and numerical methods reveal that the heat flux on wellbore exhibits an exponential reduction in the beginning, but tends to be a constant around 300 W/m (Fig. 8). This transient behavior is analogous to that of a well operation at a constant bottom-hole pressure as depicted by many literatures [47–49]. Fig. 9 shows the temperature distributions in the reservoir around a single well when the steam is circulated for 1, 10, 30, 60, 120 and 300 days, respectively. Here both analytical and numerical solutions are adopted for the prediction and they support each other. It is discovered that the phase change of bitumen has a significant impact on the temperature distribution around a single well. The temperature drops rapidly in the fusion zone with the distance to wellbore but decreases slowly in the solid zone. The heating area is wider and wider with time. Circulation for just 1 day has little effect on the reservoir except for the domain encompassing the well. When the preheating process is finished, all the reservoir covering a circular domain with a radius of 10 m is effectively heated, and the phase change interface position reaches about 2.4 m, which realizes the mutual validation with what is exhibited in Fig. 8. Fig. 10 exhibits how the temperatures at some given positions, including 1, 2, 4 and 8 m with the distance to the wellbore axis, respectively, evolve with injection time. The prediction results of analytical and numerical methods do validate each other. The temperature at the position of 1 m with the distance to wellbore axis rapidly grows, but slows down in the solid zone with injection time. The position at r = 2 m is under insufficient heating, and its temperature only reaches up to about 66 °C. The two positions left undergo no fusion of bitumen, and their temperatures fail to reach the phase change point. As studied previously, the reservoir temperature propagation with injection time can reflect the efficiency of steam circulation. Besides, what is usually focused on in the preheating stage for field engineers is the ultimate heating effect after circulation, as shown in Fig. 11. Temperature distributions in the reservoir after circulation are exhibited in Fig. 11(a), and it expresses same information as Fig. 9 when circulation for 300 days. Heat flux magnitudes in the reservoir after circulation are described in Fig. 11(b), which reveals similar information as Fig. 8 when r = 0.125 m. These results are obtained by the numerical method, and the analytical method

Fig. 10. Temperature changes T(t) with time at some given distances.

can also be used to give almost identical solutions according to preceding analyses. 5.2. Case study of well pair A-1 In this part, both I well and P well of well pair A-1 were studied. All the information shown in Section 5.1 is also suitable for this section, so it only needs to supply the information of P well of well pair A-1. The P well during the preheating stage is actually used for steam circulation just like I well but not completely, because actually an adequate differential pressure from I well to P well should be tentatively applied for improving the heat convection by possibly existed oil flows and obtaining considerable preliminary oil production at the end of the preheating. Here investigations about the process of tentatively switching to SAGD were not involved, and some relevant research progress can be found in literatures [24,27]. The vertical depth of P well is recorded as 372 m, which is 5 m away from the top surface of base rock. Tested results showed that the bottom hole vapor temperature in both I well and P well are kept at about 252 °C. Considering the phase change behavior of bitumen filling pores, the reservoir around well pair

Fig. 9. Temperature distribution T(r) with distance at some given time.

Y. Gao et al. / International Journal of Heat and Mass Transfer 144 (2019) 118650

Fig. 11. Numerical solutions of (a) temperature distribution T(r), and (b) heat flux magnitude q(r) after circulation.

Fig. 12. Temperature distribution T(x,y) when circulation for (a) 30, (b) 60, (c) 120, and (d) 300 days by analytical method.

9

10

Y. Gao et al. / International Journal of Heat and Mass Transfer 144 (2019) 118650

A-1 involving I and P wells can be divided into the fusion zone and solid zone. The wellbore size of dual wells is both 0.25 m in diameter. The circulation period is kept for about 6–12 month [2]. Other information is completely same as that of the case study of I well of well pair A-1 in Section 5.1. In the analytical method, step length Dt of 1 day is chosen to analyze the steam circulation for 300 days. When the injection time ti = iDt (1  i  300) is determined, the dimensionless parameter k(ti), flux on wellbore q(ti) and phase change interface moving speed drs(t)/dt can be obtained according Eqs. (6) and (7). In this regard, the temperature distribution T(r) at t = ti can be calculated. Then the superposition principle was used to obtain the results under circulation by a well pair. The two-dimensional Cartesian coordinate replaces the one-dimensional cylindrical coordinate, so the temperature distribution at time t is expressed in the form of T(x,y) instead of T(r). The coordinate origin is set on the midpoint of the well pair for symmetrical characteristics. As for the numerical method, a heat conduction problem involving considerable latent heat effect with time-dependent thermo-physical parameters was simulated by ABAQUS software. The studied oil sand reservoir was considered as a phase change material (solid bitumen will fuse under high temperature, but oil sand grains will not), and there is no latent heat exchange behavior in the cap rock and base rock. A four-node linear heat transfer quadrilateral element type (DC2D4, [28]) was selected for this simulation. The initial temperature of 12 °C is loaded on the whole

studied area through a predefined temperature field. The temperature conditions on the wellbores are both 252 °C, which brings into correspondence with the measured bottom hole temperature data. Calculation results by both analytical and numerical methods are exhibited in Figs. 12 and 13, respectively. In Fig. 12, an area of 10 m
20 m around the well pair is studied for its temperature distribution when the injection time is 30, 60, 120 and 300 days, respectively. In these figures, these contour lines corresponding to the phase change temperature are marked for determining the location of phase change interface. When circulation for just one month, it can be found there is only a narrow area near the well pair influenced by the preheating. Two months later, the interwell region is obviously heated, but the phase change interfaces induced by injector and producer, respectively, have not met each other. This phenomenon can be greatly improved after circulation for four months, because two phase change interfaces have already met in this moment, and the temperature propagation speed in the inter-well area is faster than before, because of the higher temperature change rate in the fusion zone than that in the solid zone. When the circulation is over, the phase change interface has already moved to a position far from the wellbore, which means the whole inter-well area is covered by the fusion zone. Therefore, the inter-well region has a great possibility in a successful thermal communication during startup. In Fig. 13, the numerical solutions for a simulated area of 30 m
35 m around the well pair are shown. When the injection

Fig. 13. Temperature distribution T(x,y) when circulation for (a) 30, (b) 60, (c) 120, and (d) 300 days by numerical method.

Y. Gao et al. / International Journal of Heat and Mass Transfer 144 (2019) 118650

is over, the cap rock and base rock have not been influenced by the preheating, and their temperatures still keep 12 °C, therefore the oil sand reservoir (covering an area 30 m
30 m) involving latent heat effect induced by the phase change of bitumen can be considered as the only studied object for simplicity. It is discovered that the solutions by numerical method are highly consistent with those by previous analytical method, so there is no need for specifically interpreting the information described in Fig. 13. In next part, the midpoint temperature will be investigated for finding a quantified criterion of terminating the steam circulation and switching to SAGD production.

5.3. Midpoint temperature of well pair A-1 The midpoint temperature between injector and producer is widely employed as a criterion of inter-well thermal communication in the startup process [1,17,19]. Usually, the lowest temperature at the midpoint between the SAGD well pair is selected as a threshold temperature, which is used as an indicator for terminating the steam circulation and switching to SAGD production [27]. The threshold temperature is discrepant in different literatures, but ranges from 80 °C to 120 °C [1,17,26,27,39,50]. Analytical and numerical solutions are exhibited in Fig. 14. If 80 °C is selected as the threshold temperature, the startup should switch to SAGD production after circulation for about 175 days. When the circulation is over, the inter-well midpoint temperature will reach 95– 104 °C. It is noted that the analytical solution is smaller than the numerical one if the injection time exceeds the switching time (around 6 months), because the temperature interference caused by phase change behaviors between I well and P well is not considered in the analytical method. Nevertheless, they agree well with each other when circulation time is shorter than 6 months. Therefore, it is acceptable in terms of temperature prediction accuracy before switching operations at the threshold time. The heat convection has a considerable influence on the temperature propagation at the end of the preheating stage [1]. How to evaluate its impact is a complicated work. The degree of its effect depends on the mobile water saturation in the formation. To estimate this effect, a new concept named heating ring was introduced by Duong. In his paper, heating rings were employed to measure the heat storage in the heated bitumen [29]. In this study, steam was considered hardly to infiltrate into the reservoir, therefore the mobile water saturation is very low. But this only occurs when there is no effective pressure difference between the well pair

120 Analytical Numerical

Temperature T (oC)

100

Threshold temperature

80 60 40

Switching to SAGD producon

20

11

before switching to production. Once the steam is pushed into the porous reservoir under a controlled relatively high inter-well pressure and the mobile condensed water generates, the effect of convection is absolutely supposed to be taken into account. At present, a shorter circulation period for a desired midpoint temperature can be achieved by geomechanics methods. Thermo-pressure-deformation coupled behavior in the SAGD process has become increasingly important in guiding the artificial operations in field practices. For example, the preliminary reservoir stimulation by the geo-mechanical method, i.e., water injection for permeability enhancement, is a pressure-deformation coupled process. The subsequent production by the propagation of steam chamber is a thermo-pressure-deformation coupled process [51– 57]. At present, thermo-pressure-deformation coupled behaviors under hot water injection in Karamay oil sand reservoirs were studied [2–4,50], which can provide a good consultation for further study if the latent heat effect is also included (then it is actually a thermo-pressure-deformation-phase change coupled process). Future research on this complicated problem will be discussed by numerical method. 6. Summary and conclusion This paper proposed an analytical model and a numerical model for predicting the thermal propagation upon preheating in a bitumen-rich oil sand reservoir. The latent heat effect resulting from the phase change of bitumen was involved in these models for thermo-phase change coupled analyses. Major conclusions can be drawn as follows: The analytical and numerical approaches both provided satisfactory predictions for the temperature propagations during the circulation stage for SAGD engineers. Application on well pair A-1 showed that solutions obtained from analytical and numerical methods support each other, and exhibit their respective advantages. The heat flux on wellbore exhibits an exponential reduction in the beginning, but tends to be a constant around 300 W/m. The interface moving speed decreases with injection time, and it takes no more than 100 days for phase change interface to reach 1.5 m with the distance to wellbore axis, and to increase up to just 2.4 m after circulation. The temperature drops rapidly in the fusion zone with the distance to wellbore but decreases slowly in the solid zone. The heating area is wider and wider with time. When circulation for just one month, there is only a narrow heated region. When circulation is over, the phase change interface has already moved to a far position, which means the whole inter-well area will be covered by the fusion zone. It should be switched to SAGD production when circulation for about 175 days. And when the circulation is over, the inter-well midpoint temperature will reach 95–104 °C. The temperature falloff test simulation widely acknowledged in literatures was conducted and the parameters derived from the semi-log relations between the shut-in time and temperature were abstracted for analyses, and the results showed a good match with those in literatures. The analysis solutions, for example, the midpoint temperature evolutions with injection time, can not only give a criterion for evaluating the efficiency of steam circulation but also point out an indicator for terminating the steam circulation and switching to SAGD production. Declaration of Competing Interest

0 0

50

100

150

200

250

Injection time t (days) Fig. 14. Midpoint temperature propagations with injection time.

300

We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that

12

Y. Gao et al. / International Journal of Heat and Mass Transfer 144 (2019) 118650

could be construed as influencing the position presented in, or the review of, the manuscript entitled ‘‘Modeling of reservoir temperature upon preheating in SAGD wells considering phase change of bitumen”. Acknowledgements We sincerely appreciate the funding provided by the National Natural Science Foundation of China (No. 51490651, 51404281, 51325402). We also owe gratitude to Xinjiang Oilfield Corporation for providing the field-collected oil sand samples and the field injection data. Appendix A. Supplementary material Supplementary data to this article can be found online at https://doi.org/10.1016/j.ijheatmasstransfer.2019.118650. References [1] A.N. Duong, T. Tomberlin, M. Cyrot, A new analytical model for conduction heating during the SAGD circulation phase, in: International Thermal Operations and Heavy Oil Symposium, Calgary, Alberta, 2008, SPE-117434-MS. [2] B. Lin, Y. Jin, H. Pang, A. Cerato, Experimental investigation on dilation mechanisms of land-facies Karamay oil sand reservoirs under water injection, Rock Mech. Rock Eng. 49 (4) (2016) 1425–1439. [3] B. Lin, Y. Jin, S. Chen, A criterion for evaluating the efficiency of water injection in oil sand reservoirs, J. Petrol. Sci. Eng. 149 (2017) 322–330. [4] B. Lin, S. Chen, Y. Jin, Evaluation of reservoir deformation induced by water injection in SAGD wells considering formation anisotropy, heterogeneity and thermal effect, J. Petrol. Sci. Eng. 157 (2017) 767–779. [5] A. Firoozabadi, Thermodynamics of Hydrocarbon Reservoirs, McGraw Hill Book Co., 1999. [6] J.G. Speight, The Chemistry and Technology of Petroleum, fifth ed., CRC Press, 2014. [7] R.M. Butler, Thermal Recovery of Oil and Bitumen, Prentice Hall Inc., 1991. [8] A. Al-Bahlani, T. Babadagli, A critical review of the status of SAGD: where are we and what is next, in: SPE Western Regional and Pacific Section AAPG Joint Meeting, Bakersfield, California, 2008, SPE-113283-MS. [9] A. Al-Bahlani, T. Babadagli, SAGD laboratory experimental and numerical simulation studies: a review of current status and future issues, J. Petrol. Eng. 68 (2009) 135–150. [10] N. Edmunds, On the difficult birth of SAGD, J. Can. Petrol. Technol. 38 (1) (1999) 14–24. [11] L. Liang, An analytical model for cyclic steaming of horizontal wells (M. Sc. thesis), Department of Petroleum Engineering, Standard University, 2005. [12] M. Irani, S. Ghannadi, Understanding the heat-transfer mechanism in the SAGD process and comparing the conduction and convection flux in bitumen reservoirs, SPE J. 18 (1) (2013) 134–145. [13] Y. Ito, S. Suzuki, Numerical simulation of the SAGD process in the Hangingstone oil sands reservoir, in: SPE Annual Technical Meeting, Calgary, Alberta, 1996, PETSOC-96-57. [14] Y. Ito, S. Suzuki, Numerical simulation of the SAGD process in the Hangingstone oil sands reservoir, J. Can. Petrol. Technol. 38 (9) (1999) 27–35. [15] S.M. Farouq-Ali, Is there life after SAGD, J. Can. Petrol. Technol. 36 (6) (1997) 20–23. [16] Y. Ito, Discussion of on the difficult birth of SAGD, J. Can. Petrol. Technol. 38 (5) (1999) 22–23. [17] L. Zhu, F. Zeng, G. Zhao, A. Duong, Using transient temperature analysis to evaluate steam circulation in SAGD start-up processes, J. Petrol. Sci. Eng. 100 (2012) 131–145. [18] B. Moini, N. Edmunds, Quantifying heat requirements for SAGD startup phase: steam injection, electrical heating, J. Can. Petrol. Technol 52 (2) (2013) 89–94. [19] S. Ghannadi, M. Irani, R. Chalaturnyk, Overview of performance and analytical modeling techniques for electromagnetic heating and applications to steamassisted-gravity-drainage process startup, SPE J. 21 (2) (2016) 311–333. [20] R.M. Butler, D.J. Stephens, The gravity drainage of steam-heated heavy oil to parallel horizontal wells, J. Can. Petrol. Technol. 20 (2) (1981) 90–96. [21] J.C. Reis, A steam-assisted gravity drainage model for tar sands: linear geometry, J. Can. Petrol. Technol. 31 (10) (1992) 14–20. [22] S. Akin, Mathematical modeling of steam-assisted gravity drainage, SPE Reserv. Eval. Eng. 8 (5) (2005) 372–376. [23] A. Azad, R.J. Chalaturnyk, A mathematical improvement to SAGD using geomechanical modeling, J. Can. Petrol. Technol. 49 (10) (2010) 53–64.

[24] G. Parmar, L. Zhao, J. Graham, Start-up of SAGD wells: history match, wellbore design and operation, J. Can. Petrol. Technol. 48 (1) (2009) 42–48. [25] J. Yuan, McFarlane, Evaluation of steam circulation strategies for SAGD startup, J. Can. Petrol. Technol. 50 (1) (2011) 20–32. [26] M. Anderson, D. Kennedy, SAGD startup: leaving the heat in the reservoir, in: SPE Heavy Oil Conference Canada, Calgary, Alberta, 2012, SPE-157918-MS. [27] J. Zhang, Y. Fan, B. Xu, B. Yang, Y. Yuan, Y. Yu, Steam circulation strategies for SAGD wells after geomechanical dilation start-up, in: SPE Canada Heavy Oil Technical Conference, Calgary, Alberta, 2016, SPE-180705-MS. [28] ABAQUS, Abaqus theory guide, version 2016, Dassault Systemes Simulia Corp., 2016. [29] A.N. Duong, Thermal transient analysis applied to horizontal wells, in: International Thermal Operations and Heavy Oil Symposium, Calgary, Alberta, 2018, SPE-117435-MS. [30] L.G. Hepler, C. His, AOSTRA technical handbook on oil sands, bitumen and heavy oils, Alberta Oil Sands Technology and Research Authority, 1989. [31] N. Ozisik, Heat Conduction, second ed., John Wiley & Sons Inc, 2014. [32] Y. Gao, M. Chen, B. Lin, Y. Jin, Experimental investigation on compressibility of Karamay oil sands under water injection, in: 51st U.S. Rock Mechanics/ Geomechanics Symposium, San Francisco, California, 2017, ARMA-2017-0597. [33] P.M. Collins, Geomechanical effects on the SAGD process, in: SPE International Thermal Operations and Heavy Oil Symposium, Calgary, Alberta, 2005, SPE97905-MS. [34] S. Himran, A. Suwono, G.A. Mansoori, Characterization of alkanes and paraffin waxes for application as phase change energy storage medium, Energy Source 16 (1) (1994) 117–128. [35] S. Li, J. Wang, H. Tan, Z. Wu, Study of extraction and pyrolysis of Chinese oil sands, Fuel 74 (8) (1995) 1191–1193. [36] P.M. Collins, The false lucre of low-pressure SAGD, J. Can. Petrol. Technol. 46 (1) (2007) 20–27. [37] P.M. Collins, Geomechanical effects on the SAGD process, SPE Reserv. Eval. Eng. 10 (4) (2007) 367–375. [38] P. Li, R.J. Chalaturnyk, M. Polikar, Issues with reservoir geomechanical simulations of the SAGD process, J. Can. Petrol. Technol. 43 (5) (2004) 30–40. [39] L. Guindon, Heating strategies for steam-assisted-gravity-drainage startup phase, J. Can. Petrol. Technol. 54 (2) (2015) 81–84. [40] A. Lazaridis, A numerical solution of the multidimensional solidification (or melting) problem, Int. J. Heat Mass Tran. 13 (9) (1970) 1459–1477. [41] Y.K. Chuang, J. Szekely, On the use of Green’s functions for solving melting or solidification problems, Int. J. Heat Mass Tran. 14 (9) (1971) 1285–1294. [42] Y.K. Chuang, J. Szekely, The use of green’s functions for solving melting or solidification problems in the cylindrical coordinate system, Int. J. Heat Mass Tran. 15 (5) (1972) 1171–1174. [43] M.S. Selim, R.C. Seagrave, Solution of moving-boundary transport problems in finite media by integral transforms. I. Problems with a plane moving boundary, Ind. Eng. Chem. Fund. 12 (1) (1973) 1–8. [44] M.S. Selim, R.C. Seagrave, Solution of moving-boundary transport problems in finite media by integral transforms. II. Problems with a cylindrical or spherical moving boundary, Ind. Eng. Chem. Fund. 12 (1) (1973) 1171–1174. [45] S. Paterson, Propagation of a boundary of fusion, P. Glasgow Math. Assoc. 1 (1952) 42–47. [46] A.E. Green, P.M. Naghdi, A general theory of an elastic-plastic continuum, Arch. Ration. Mech. An. 18 (4) (1965) 251–281. [47] A.F. Van-Everdingen, W. Hurst, The application of the Laplace transformation to flow problem in reservoirs, J. Petrol Tech. 1 (12) (1949) 305–324. [48] C.E. Jacob, S.W. Lohman, Nonsteady flow to a well of constant drawdown in an extensive aquifer, T. Am. Geophys. Union 33 (4) (1952) 559–569. [49] M.J. Fetkovich, Decline curve analysis using type curves, J. Petrol. Tech. 32 (6) (1980) 1065–1077. [50] X. Yuan, S. Dou, J. Zhang, S. Chen, B. Xu, Consideration of geomechanics for insitu bitumen recovery in Xinjiang, China, in: SPE Heavy Oil ConferenceCanada, Calgary, Alberta, 2013, SPE-165414-MS. [51] R.J. Chalaturnyk, P. Li, When is it important to consider geomechanics in SAGD operations, J. Can. Petrol. Tech. 43 (4) (2004) 53–61. [52] P. Li, R.J. Chalaturnyk, Geomechanical model of oil sands, in: SPE International Thermal Operations and Heavy Oil Symposium, Calgary, Alberta, 2005, SPE97949-MS. [53] P. Li, R.J. Chalaturnyk, History match of the UTF phase a project with coupled reservoir geomechanical simulation, J. Can. Petrol. Tech. 48 (1) (2009) 29–35. [54] D.P. Yale, T. Mayer, J. Wang, Geomechanics of oil sand under injection, in: 44th U.S. Rock Mechanics Symposium and 5th U.S.-Canada Rock Mechanics Symposium, Salt Lake City, Utah, 2010, ARMA-10-257. [55] Y. Yuan, B. Xu, B. Yang, Geomechanics for the thermal stimulation of heavy oil reservoirs-Canadian experience, in: SPE heavy oil conference and exhibition, Kuwait City, Kuwait, 2011, SPE-150293. [56] Y. Gao, M. Chen, B. Lin, Y. Jin, An analytical model of hydraulic dilation area for Karamay oil sand reservoir under water injection in SAGD wells, J. Petrol. Sci. Eng. 179 (2019) 1090–1101. [57] Y. Xia, Y. Jin, M. Chen, K. Chen, An enriched approach for modeling multiscale discrete- fracture/matrix interaction for unconventional-reservoir simulations, SPE J. 24 (1) (2019) 349–374.