Influence of temperature-viscosity behaviors of Karamay oil sand bitumen on the geomechanics in the SAGD process

Viscosity-temperature behavior is important for the thermal recovery of bitumen: firstly, it is a criterion for evaluating the evolution of mobility with temperature, which affects oil production; and secondly, it determines the role of bitumen in the geomechanics in oil sand reservoir. To quantitatively evaluate the effect of temperature on dynamic viscosity, some experimental and modeling investigations were conducted. The viscosity of Karamay oil sand bitumen over a temperature range of 20 ~ 350 °C was measured, and the viscosity-temperature behavior was simulated by a modified three-parameter model using the least square method. It was found that the viscosity of Karamay bitumen drops sharply with the temperature increase. The effect of bitumen viscosity on the geomechanics in the SAGD process was discussed in terms of reservoir deformation, fluid flow, and heat transfer behaviors. The reservoir deformation, fluid flow, and heat transfer behaviors at varying bitumen viscosity show significant differences in the drained, partially drained, and undrained geomechanical zones. The structure’s bulk modulus, effective stress, and volumetric strain in the drained zone are lower than those in the undrained zone; while the oil mobility and Peclet number show the opposite tendency. The changes in the structure’s bulk modulus, effective stress, volumetric strain, oil mobility, and Peclet number due to the phase change of bitumen increase with the increase of oil saturation. This study can provide the field engineers with guidance for the design of a proper oil recovery scheme before its implementation, and with a useful relation for the coupled thermal-flow-structure analyses.


Introduction
Heavy oil including extra-heavy oil or bitumen, mostly located in western Canada and eastern Venezuela Basin, is one of the most widely known unconventional oil resources in the world (Doan et al. 2012). China also has very rich heavy oil resources, most of which are deposited in the shallow oil sand reservoirs of the northwest margin of Junggar Basin, Xinjiang province, with an estimated resource of about 360 million tons (Yuan et al. 2013). Oil sands are these unconsolidated sandstones that contain an extremely viscous hydrocarbon that is not recoverable in its natural state by conventional oil well production methods including currently used enhanced recovery techniques (Speight 2014). This organic matrix in oil sands known as oil sand bitumen fills these pores and crevices of sandstone (Hoiberg 1964). The extremely viscous oil sand bitumen occurs in a solid or near-solid state, and is typically incapable of mobility under reservoir conditions. For example, the oil sand bitumen located in Alberta, Canada, is not mobile in the deposit and requires extreme methods of recovery to recover the bitumen (Speight 2014). By comparison, the dramatically high viscosity of Karamay oil sand bitumen ranges from 10 5 to 10 6 mPa·s under reservoir conditions, and it cannot flow at all when the reservoir temperature is lower than 80 °C (Yuan et al. 2013). Up to date, the steam-assisted gravity drainage (SAGD) technique has been implemented at a large scale in the development of oil sand bitumen reservoirs in Karamay city, and a lot of field practices have been conducted in the Xinjiang Fengcheng Oilfield, the largest extra-heavy oilfield in China.

3
The oil sand reservoir temperature has a very significant change in both temporal and spatial distributions during the SAGD process. In the SAGD rapid startup phase, known as water injection into the oil sand reservoir to increase the porosity and permeability for reducing the preheating period, the reservoir temperature in Karamay is just about at a range of 18.5 ~ 80 °C (Yuan et al. 2013). However, the reservoir temperature near the steam chamber in the SAGD production phase can reach more than 250 °C. It is a common sense that the bitumen viscosity is extremely sensitive to temperature, so the bitumen viscosity in the whole SAGD process will undergo a very great evolution because of temperature substantial increase. It is essential for field engineers to know how viscous the oil sand bitumen is under varying reservoir temperatures caused by steam circulation or injection. On one hand, if the heated bitumen has mobility in the reservoir, its viscosity will determine its flow capability and oil production just like the common light oils, which can be expressed by Darcy's law. On the other hand, if the reservoir temperature is very low, the bitumen viscosity will determine its phase state (solid/semisolid or liquid) because of the presence of a large proportion of solid waxes and asphaltenes in virgin bitumen (Speight 2014;Firoozabadi 1999). In this regard, the bitumen plays different roles in different SAGD phases. In the SAGD startup phase under water injection, the solid bitumen can be viewed as the skeleton and cement; but in the steam circulation and production phases, the melting bitumen plays a role as the pore fluid just like the water and gas. The Karamay oil sand bitumen is more viscous than Canadian bitumen, so the bitumen in Karamay oil sands is always regarded as the solid skeleton, because it cannot be excluded from the pores, just as shown in Fig. 1 (Li et al. 2015). However, it is not frequent to find the solid bitumen from microstructures of Athabasca oil sands in Canada, as shown in Fig. 1 (Dusseault and Morgenstern 1978). In the SAGD process, when the viscosity of crude oil drops to a critical value of 1000 mPa s, the drained zone forms and it is not very difficult for the pore fluids to flow under some driving forces (Li et al. 2004). Therefore, the study on the viscositytemperature behaviors is very important.
Up to date, the investigations on viscosity-temperature behaviors have been conducted by laboratory experiments and mathematical modeling. For the experiments, measurement instruments can be used to measure the viscosity under different temperatures according to the test standards, so as to obtain the viscosity-temperature data. Conventional test approaches for obtaining the viscosity-temperature data can be achieved by ASTM standards (e.g., ASTM 2010, 2015) using relatively simple instruments like the rotational viscometer and constant pressure viscometer. Considering the high temperature in the process of oil sand bitumen fusion under steam injection, these conventional methods are not suitable for their narrow test temperature ranges. In this regard, many advanced automatic high-temperature rheometers have sprung up for their convenience and high test precision. For the modeling study, these approaches include Walther's method, modified Walther's methods (Svrcek and Mehrotra 1988;Wright 1968), and other methods like the two-parameter Andrade equation (Hepler and Hsi 1989), double-logarithmic equation (Speight 2014), and empirical equation (AOSTRA 1984). ASTM D341-17 (2017) details the standard procedure for plotting viscosity-temperature charts that ascertain the kinematic viscosity of petroleum oil or liquid hydrocarbon at any temperature within a limited range, provided that the kinematic viscosities at two temperatures are known. Unfortunately, neither the experiments nor modeling has been intensively studied for the Karamay oil sand bitumen. These experimental and modeling methods above are not suitable for the Karamay oil sand bitumen for its dramatically high viscosity and special composition (high-content asphaltenes and resins).
Up to now, the viscosity-temperature behaviors for oil sand bitumen have been rarely discussed. Ukwuoma and Ademodi (1999) measured the apparent viscosities of Nigerian oil sand bitumen over a temperature range of 50 ~ 110 °C, and found that the bitumen became more Newtonian in the higher temperature region. Svrcek and Mehrotra (1988) published data on the viscosity of Marguerite lake and Athabasca bitumen which observed the viscosity decreased with increasing temperature and exhibited near-zero slopes in the high-temperature region. AOSTRA (1984) listed lots of test data for Cold Lake oil sand bitumen viscosity-temperature behaviors. Hepler and Hsi (1989) discussed the relation between viscosity and temperature of Canadian bitumen (Athabasca, Cold Lake, Peace River, and Wabasca). Bazyleva et al. (2010) studied the rheological properties including the viscosity over the temperature range of 200 ~ 410 K. Athabasca bitumen and Maya crude oil were found to be solid-like materials up to 1 3 260 ~ 280 K and 230 ~ 240 K, respectively, and they are both Newtonian at higher temperatures.
One of the most important reasons we focus on the viscosity-temperature relation is that the viscosity change of bitumen with temperature affects the reservoir deformation, fluid flow, and heat transfer behaviors. As shown in Fig. 2, the reservoir during SAGD can be divided into three geo-mechanical zones according to the viscosity of bitumen (Li et al. 2004). Firstly, the deformation mechanisms of the structure in these zones are different. In the undrained zone, the extremely viscous solid bitumen being a kind of cement is filled in the intergranular pores ). In the partially drained zone, the mobile bitumen can flow under certain pressure gradients, but it is still hard for the mixed water-bitumen fluids to be squeezed out of pores. Therefore, the mixed fluids in this zone play a role as same as that in the undrained zone. In the drained zone, the mobility of bitumen is similar to water so that the grains can be compacted to contact each other. In this case, the elastic deformation of the structure is predominant. As for the effective stress variation, the stress in the drained zone is theoretically dependent on the steam injection pressure that can lead to pore pressure increase due to pore fluid flow and thermal expansion that induces both pore pressure and total stress to increase. In the partially drained zone, the pore pressure change and effective stress variation are also dependent on the same factors like that in the drained zone, but the magnitude of the variation is not the same as that in the drained zone. While in the undrained zone, geomechanical shear in the undrained zone creates a condition of no volume change but leads to changes in pore pressure (Li et al. 2004). The effective stress variation in this zone is also affected by the total stress change induced by thermal expansion in the steam chamber, drained zone, and partially drained zone (Li et al. 2004). As shown in Fig. 3, in some studies, the deformations of the structure and the fluids are evaluated by viscoelastic constitutive models where the effective stress needs not be considered (Gbadam and Frimpong 2017;Chang and Meegoda 1997).
Secondly, the fluid flow behavior can be affected by the viscosity of bitumen, as shown in Fig. 4. In the drained zone, there are two liquid phases (gas is neglected), and the mixed bitumen and water flow very well in the connected pores. While in the undrained zone, the rich bitumen causes plenty of unconnected pores, which means the only liquid phase, water, is very difficult to flow. In the reservoir geomechanical simulation, the permeability used in the drained zone should be the permeability when bitumen being pore fluids, while the permeability used in the undrained zone should be the permeability to water when considering the bitumen as a part of the skeleton (Gao et al. 2019a). Besides, the viscosity changes in the drained zone will also affect the output of oil according to Darcy's law. Similarly, the porosity used in the drained zone should be calculated by the space of both bitumen and water, while the porosity used in the undrained zone is only determined by the space of water (Gao et al. 2019b). In the partially drained zone, the flow mechanism is almost the same as that in the drained zone. In the drained zone, the permeability is much higher prior to the pore collapse induced by the oil-water fluids being squeezed out of pores. The viscosity of fluids in the partially drained zone is higher than that in the drained zone because of the lower temperatures.
Thirdly, the heat transfer behavior can be influenced by bitumen viscosity. As shown in Fig. 5, the bitumen is one of the typical amorphous materials, whose phase change is performed in a temperature range (Batzle and Han 2004;Han et al. 2005Han et al. , 2006. The bitumen at a temperature lower than the glass transition point keeps a solid state, while the bitumen at a temperature higher than the melting point is   in a liquid state. When the temperature is between the two critical temperatures, the bitumen shows a quasi-solid state. In reality, the Karamay bitumen contains asphaltenes, resins, aromatics, and saturates, in which the latter two are the main components of waxes. So there is a latent heat effect in the process of bitumen fusion during the phase change of waxes. Apart from the latent heat, the heat transfer parameters in three geomechanical zones are much different (Gao et al. 2019c). For example, the thermal conductivity and thermal diffusivity in the drained zone are just half of those in the undrained zone (Hepler and Hsi 1989).
At present, few experimental and modeling approaches are available in the literature to exhibit and describe the viscosity-temperature behavior of Karamay oil sand bitumen under the elevated temperature as high as 250 °C corresponding to actual reservoir temperature under steam injection. The geomechanics induced by the viscosity change of bitumen in the SAGD process is also seldom mentioned. These situations hinder the design of a proper oil production scheme before its implementation and the analyses of the coupled thermal-flow-structure process. In this contribution, the viscosity of Karamay oil sand bitumen over a temperature range of 20 ~ 350 °C was measured by an advanced rheometer, and the viscosity-temperature response derived from the test data was simulated by a modified three-parameter model using the least square method. Based on the study of viscosity-temperature relation, the effect of bitumen viscosity on geomechanics in the SAGD process was discussed in terms of reservoir deformation, fluid flow, and heat transfer behaviors. Some useful quantitative models were proposed to evaluate the geomechanical responses under different viscosities of bitumen due to the thermal stimulation.

Viscosity-temperature relation for Karamay bitumen
The oil sand reservoirs studied were deposited during the Jurassic period in the Qigu and Badaowang formations. 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 Junggar Basin. The oil sand bitumen in three areas (areas A, B, and C) in the Fengcheng Oilfield was collected for the experiments. Region A possesses a porosity of 28.6%, an absolute permeability of 810 mD, and an oil saturation of 66.7%. Region B owns a porosity of 33.1%, and an absolute permeability of 1175 mD. Region C has a porosity of 30.3%, an absolute permeability of 2951 mD, and an oil saturation of 73.2%. These data above are all averages.
In areas A, B, and C, bitumen materials were collected from six wells (A1 ~ A6), four wells (B1 ~ B4), and two wells (C1 ~ C2), respectively. For eliminating the possible differences of rheological properties of crude bitumen samples even collected from one well, bitumen samples collected at different times for every well in area B were also used for the viscosity tests. The bitumen samples for wells B1, B2, B3, and B4 were collected from the same formations for six, five, two, and two times, respectively. Although these samples were collected from the same formation and well, their density, contents of wax, precipitate, water, and freezing point are somewhat different.
The viscosity of the sample was assessed on a rheometer (StressTech, Reologica Instruments, Sweden). The rheometer can measure the bitumen whose viscosity ranges from 0.5 to 10 8 mPa s, under a working temperature of − 20 ~ 350 °C and a working pressure as high as 40 MPa. The test procedures can be under the guidance of the instrument manual or other published literature which used this rheometer for the viscosity test (Petersson et al. 2008).

Test data
• Temperature-viscosity data of area A The test temperature ranges from 20 to 350 °C for the bitumen materials collected from six wells in this area, and the viscosity of each sample under about every 10 °C was tested. The test data are shown in the semi-logarithmic coordinate in Fig. 6.   • Temperature-viscosity data of area B For the bitumen materials collected from four wells in this area, the test temperature ranges from 25 to 200 °C, and the viscosity of each sample was tested at about every 10 °C. The specific test data are shown in the semilogarithmic coordinate in Fig. 7.
• Temperature-viscosity data of area C For the bitumen materials collected from two wells in this area, the test temperature ranges from 25 to 200 °C, and the viscosity of each sample was tested at about every 10 °C. The specific test data are shown in the semilogarithmic coordinate in Fig. 8 As shown in Figs. 6 7 and 8, the viscosity of Karamay bitumen drops sharply with the temperature increase. Figure 6 shows that, in area A, the viscosity of bitumen on reservoir conditions in A6 reaches about 6 × 10 4 mPa s, while the viscosity in A2 dramatically reaches as high as 2 × 10 6 mPa·s, which indicates a significant spatial distribution of bitumen viscosity in this area. The viscosity drops to a value of less than 10 mPa·s when the temperature is heated to more than 250 °C, which is favorable for oil flow and production. Figure 7 shows that, in region B, the initial viscosity of bitumen collected from different wells in area B is in a range of 10 6 ~ 5 × 10 6 mPa·s. Viscosity shows quite a similar range of variation at the low end of temperature. Figure 8 shows that, in region C, the viscosity of bitumen in C1 is higher than that in C2 at the same temperature lower than 100 °C, and there is little difference when the temperature is heated to more than 100 °C. The super high viscosity means that the Karamay bitumen cannot flow, and just stays in the initial position during SAGD startup phase under water injection at the temperature of 20 ~ 80 °C. The bitumen in the process of water injection plays a role of skeleton and cement, and it can affect the shear dilation capability to some degree (Lin et al. 2016).

Regression of viscosity-temperature test data
Equation (1) was usually used to depict the relation of viscosity and temperature for oils (Svrcek and Mehrotra 1988;Wright 1968;Butler 1991).
where μ is the dynamic viscosity, mPa s; T is the temperature, °C; m and n are two parameters representing oil properties.
The high-viscous oil sand bitumen shows very complex viscosity-temperature responses according to the test data. A modified model with three parameters was proposed in this paper to describe the relation between dynamic viscosity and temperature.
where a, b and c are three parameters representing bitumen properties.
As shown in Figs. 6 7 and 8, there are 7 ~ 17 sets of test data for each temperature-viscosity data regression, but there are only three undetermined parameters. In other words, the equation number (N) is more than the number of undetermined parameters. Therefore, the determination of three parameters a, b and c can be transformed into a problem: the solution of nonlinear overdetermined equations. This problem can be expressed as (1) log 10 log 10 ( + 0.7) = m log 10 (T + 273) + n (2) log 10 log 10 ( + a) = b log 10 (T + 273) + c The least square method was used to obtain an optimum solution, as shown in Eq. (4). Now, the question is to find optional solutions of a, b, and c to make g a minimum value. The solution steps to find these optional solutions in Eq. (4) is as follows: (1) Input these test data (T 1 , μ 1 ), (T 2 , μ 2 ), …, (T i , μ i ) in MATLAB; (2) Use the function of "fminsearch" in MATLAB, and find the solution of a, b, and c; (3) Input the solution of a, b, and c in Eq. (4), and calculate the solution of g.
By calculation, all the regressions agree well with the test data, except for wells A3, A4, and A6, whose fitting errors are more than 0.001. The regression curves are shown in Figs. 6 ~ 8. As shown in Fig. 6, the modeling effects for A3, A4, and A6 are not very good, and there is a critical temperature or inflection temperature near 100 °C. The water content of these samples was determined by Karl Fisher titration, and found the negligible water content is not the reason. Taxler (1961) said that the viscosity of bitumen ranges from Newtonian to highly non-Newtonian liquids depending on source and composition. Guthrie (1960) observed that the viscosity of oil sand bitumen is Newtonian when it contains more saturates and non-Newtonian when it contains high asphaltenes. To evidence that the bitumen has a non-Newtonian property, bitumen material collected from well A3 with the largest fitting error was used to analyze the content of asphaltenes. To know the components of the bitumen in well A3, the mass contents of asphaltenes, resins, aromatics, and saturates were tested according to ASTM D4124-09 (ASTM 2018). The test results are shown in Fig. 9.
The macromolecular substances including asphaltenes and resins of bitumen in the well A3 account for 20.87% and 27.42, respectively, which means a strong non-Newtonian property. To improve the regression accuracy, the piecewise regression was used for the viscosity-temperature response for the bitumen in wells A3, A4, and A6. The regression parameters are shown in Table 1.
It is worth noting that the magnitude of "a" in the left and right panels in Table 1 cannot really exhibit whether the overfitting occurs, because the independent variables for the two kinds of fitting equations are different. The piecewise regression curves and test data are shown in Fig. 10.
From Table 1 and Fig. 10, all the regression curves agree well with the test data for bitumen in wells A3, A4 and A6, (4) g = N ∑ i log 10 log 10 i + a − b log 10 T i + 273 − c 2 whose fitting error is no more than 0.001. The modeling effects for the three wells near inflection temperature are very good.
In the next step, the test data about the asphaltenes content of the other subset is desired to make the explanation more convincing. Ukwuoma and Ademodi (1999) used the bitumen sample with a low asphaltene content of 10.6%, and found their viscosity-temperature relation shows a line in the semi-log coordinate instead of two lines in our study.

Statistical characteristics of the viscosity-temperature relation
According to the studies of Svrcek and Mehrotra (1988) and Butler (1991), the dynamic viscosity in Eqs. (1) ~ (2) can be also replaced by the kinematic viscosity. They concluded that m and b in Eq. (1) or b and c in Eq. (2) are both linearly related, no matter the viscosity is dynamic or kinematic. Svrcek and Mehrotra (1988) concluded that m in Eq. (1) is -3.5556 when the viscosity is the kinematic viscosity. Butler (1991) analyzed the viscosity-temperature behaviors of heavy oils in Canada, USA, and South America, and found the relations of m = 0.3249-0.4106b using kinematic viscosity and of b = 0.3464-0.4127c using dynamic viscosity, just as shown in Fig. 11. In this study, all the modeling parameters of b and c for Karamay oil sand bitumen are plotted in Fig. 11, and a linear relation of b = − 0.41035c + 0.37217 can be used to fit all these parameters well. It is shown that the highly linear relation of b and c for Karamay oil sand bitumen is very similar to that for the bitumen in Canada, USA, and South America. The two lines are parallel to each other, and their intercepts have not much difference. The two lines can even be expressed by one equation if the modeling accuracy requirement is not very high. These highly linear  relations were interpreted as similar structures by Butler (1991).
In this regard, the viscosity-temperature behavior of Karamay oil sand bitumen can be described by a simplified model with only two parameters.
The parameter a was always viewed as a constant of 0.7 or 0.8 (Svrcek and Mehrotra 1988;Wright 1968;Butler 1991). In this paper, all the model parameters of a, b, and c as well as the modeling error log 10 (g) are exhibited in Fig. 12, where their average values have been calculated. From Fig. 12, the average values of a, b, and c as well as the modeling error log 10 (g) are − 0.71196, − 3.35549, 9.084057, and − 3.73806, respectively. It is noted that the average a for Karamay oil sand bitumen is about − 0.7 instead of 0.7 or 0.8 given in the literature (Svrcek and Mehrotra 1988;Wright 1968;Butler 1991). The viscosity behavior for the bitumen in Karamay city is different with that in Canada, USA, and south America. The a is a coefficient depending on bitumen materials, and it is a changeable coefficient.
According to the study of Svrcek and Mehrotra (1988) and Butler (1991), five effective figures are needed for the coefficients in the empirical equations. Let a = − 0.71196, and Eq. (5) can be further reduced into a model with a single parameter.
The average model can be used to generally predict the viscosity-temperature behavior if there is not any information for Karamay oil sand bitumen. The viscosity-temperature curve is shown in Fig. 13.

Effect of viscosity on the geomechanics in the SAGD process
As discussed above, the viscosity of Karamay bitumen undergoes dramatic changes when temperature increases in the SAGD process. Therefore, the roles of bitumen can be very different under different viscosities. In reality, the viscosity evolution must be taken into careful consideration in the analyses of geomechanics in the SAGD process, just as shown in Fig. 2. In this paper, the behaviors of reservoir deformation, fluid flow, and heat transfer at different viscosities in different geomechanical zones were quantitatively investigated.

Undrained zone
In this study, the structure for oil sand is defined as the combination of both the rock skeleton (mineral skeleton) and the immobile bitumen (non-mineral skeleton) in the undrained (7) log 10 log 10 ( − 0.71196) = −3.35549 log 10 (T + 273) + 9.084057 where ij is the total stress; K fr1 is the bulk modulus of the structure in the undrained zone; K s1 is the bulk modulus of the matrix in the undrained zone; p w is the pore pressure of water; δ ij is the Kronecker symbol, which is exhibited as In this paper, the bulk moduli of the structure and matrix were evaluated by Hashin's method (Hashin 1960), where the bulk modulus of a randomly heterogeneous material can be theoretically calculated. In this regard, the bulk modulus of the structure in the undrained zone can be written as where K s and v s are the bulk modulus and Poison's ratio of the mineral solid matrix, respectively; K b is the bulk modulus of the bitumen; w is the volume fraction of water; b is the volume fraction of the bitumen; is the total of w and b . Similarly, the bulk modulus of the matrix in the undrained zone can be written as where ′ b is the volume fraction of bitumen in the matrix, and it is expressed as Using the equations above, the effective stress applied on the structure is From Eq. (13), the ̂i j can evolve with the change in K b in the undrained zone. In this paper, the relation between the bulk modulus and the viscosity of bitumen was established to observe the effect of temperature-viscosity behaviors on the effective stress applied on the structure.
The relation between bitumen bulk modulus and temperature can be obtained by Eastwood's study (Eastwood 1993), where a linear relation was proven, as shown in Fig. 14.
According to the data in Fig. 14, the relation between the bulk modulus of bitumen and temperature can be expressed as Considering Eq. (7), the temperature can be written as Therefore, the relation between the bulk modulus and viscosity of bitumen is Figure 15 exhibits the evolutions of bitumen's bulk modulus with its viscosity. It can be known that the bulk modulus changes sharply when the viscosity increases under a low viscosity, and the change rate of bulk modulus slows down when the viscosity increasing. In this regard, substituting Eq. (16) into Eq. (13), the change in effective stress with viscosity can be calculated. Meanwhile, the deformation of the structure can also be determined. Using the volumetric deformation as an example, there is where v1 is the volumetric strain of the structure in the undrained zone; p ′ 1 is the mean effective stress applied on the structure.

Drained zone
When the bitumen's viscosity is less than 1000 mPa·s, the bitumen and water are both mobile fluids. Therefore, the oil sands in the drained zone can be viewed as a singleporosity medium saturated by two liquid phases. In this situation, the effective stress applied on the matrix can be calculated by (Chen et al. 1999;Gao et al. 2019d) where K fr2 is the structure in the drained zone; p b is the pore pressure of fluid bitumen.
Using Hashin's method, the structure in the drained zone can be expressed as Because water and bitumen are mixed and flow together in the same pores, the pore pressures of water and bitumen can be reasonably treated as the same value p, which can be exhibited as

Substituting Eqs. (19) and (20) into Eq. (18), there is
The volumetric strain of the structure in the drained zone can be calculated by where v2 is the volumetric strain of the structure in the drained zone; p ′ 2 is the mean effective stress applied on the structure.

Partially drained zone
It is assumed that the bulk modulus, effective stress, and deformation of structure in the partially drained zone can be approximately evaluated by the linear interpolation between the values in the drained and undrained zones (Li and Chalaturnyk 2003;Li et al. 2004;Gao et al. 2020b;Gao and Chen 2020b). The same methods are used in the latter discussions about the fluid flow and heat transfer behaviors in this zone. Therefore, the bulk modulus, effective stress, and deformation of the structure in this zone can be written as StaƟc test data (Eastwood, 1993) Ultrasonic test data (Eastwood, 1993) Regression where G represents the bulk modulus, effective stress, or deformation; the subscripts p, d, and u mean the partially drained zone, drained zone, and undrained zone, respectively; μ c1 = 1000 mPa·s and μ c2 = 20,000 mPa·s are the two critical viscosity of bitumen.
Using the parameters listed in Table 2, the evolutions of the bulk modulus, effective stress, and deformation of the structure with the viscosity can be depicted. Figure 16a exhibits the changes in structure's bulk modulus with the viscosity of bitumen. It can be observed that the structure's bulk modulus shows different evolutions in different geomechanical zones. In the drained zone, the structure's bulk modulus is not affected by the viscosity; while in the undrained zone, the structure's bulk modulus can increase with the increase of viscosity. In the partially drained zone, the bulk modulus changes sharply with the viscosity due to the fusion of bitumen. Figure 16b displays the evolutions of the effective stress applied on the structure with the viscosity of bitumen. The − c1 effective stress applied on the structure keeps constant in the drained zone, but drops slightly with the increase of viscosity in the undrained zone. The released effective stress occurring in the partially drained zone due to the fusion of bitumen reaches about 0.9 MPa. Figure 16c shows the relation between the volumetric strain of the structure and bitumen's viscosity. The structure's volumetric strain of the structure is constant in the drained zone, but decreases with the increase of viscosity in the undrained zone. The volumetric strain of the structure declines by about 0.003% due to the fusion of bitumen. The volumetric strain change is small, but the relative change seems considerable because of a small initial value. Figure 16 shows the changes in structure's bulk modulus, the effective stress applied on the structure, and the volumetric strain of the structure with the viscosity of bitumen. It can be observed that there are sudden changes when the viscosity changes from 20,000 mPa·s to 1000 mPa·s. For instance, the sudden change in Fig. 16b means the released effective stress due to bitumen's fusion.
In the Karamay Oilfield, the volume contents of bitumen or the oil saturations in different oil sand reservoirs are very different (Lin et al. 2016). For example, the bitumen-rich oil sand possesses a volume fraction of bitumen as high as about 30% (oil saturation is about 70%), while the mud-rich oil sand only owns a volume fraction of bitumen as low as 10% (oil saturation is about 40%). It is assumed that the volume content of water is 18% and the volume content of bitumen is from 8 to 28%. According to Eqs. (10), (13), (17) ~ (19), and (22), the fusion-induced bulk modulus of the structure, released effective stress and related volumetric strain can evolve with the change in the oil saturation. Figure 17 shows the evolutions of the bulk moduli, effective stresses, and volumetric strains with the oil saturation before and after the fusion of bitumen, respectively.
As shown in Fig. 17a, the structure's bulk moduli before and after the fusion of bitumen decline a lot with the increase Table 2 Parameters used to obtain the evolutions of the bulk modulus, effective stress, and deformation of the structure with the viscosity  Fig. 16 a Changes in structure's bulk modulus with bitumen viscosity, b evolutions of the effective stress applied on structure with bitumen viscosity, and c relations between the volumetric strain of structure and bitumen viscosity of oil saturation. This means that the bitumen-rich oil sand has a lower bulk modulus than the mud-rich oil sand. The change in structure's bulk modulus due to the fusion of bitumen can slightly increase with the improvement of oil saturation. Figure 17b exhibits that the effective stress applied on the structure rises with the increase of oil saturation before the fusion of bitumen, while it shows an opposite trend after the fusion of bitumen. The change in the effective stress applied on the structure due to the fusion of bitumen can increase with the increase of oil saturation. Figure 17c shows that the volumetric strains before and after the fusion of bitumen will go up with the increase of oil saturation. The change in the volumetric strain due to the fusion of bitumen can increase with the increase of oil saturation.

Undrained zone
The steady state is assumed throughout this study, and the transient period is assumed to be negligible. This assumption is frequently used in the literature (Chalaturnyk and Li 2004;Li and Chalaturnyk 2009). It is assumed that the effect of temperature on permeability is negligible. According to Eq. (7), when the viscosity of bitumen is 20,000 mPa·s, the temperature is 57 °C. In this zone, water is the only fluid that can flow, so the seepage process can be expressed as where V w is the velocity of water; K w1 is the effective permeability to water measured at reservoir temperature (20 °C) in laboratory (an effective permeability to water in the undrained zone), and K w1 = 2.886 mD according to the test data (Lin et al. 2016); μ w is the viscosity of water, mPa·s; p w is the pressure of water; M w is the mobility of water.
The relation between the viscosity of water and temperature can be written as (Korson et al. 1969)

Drained zone
When the temperature exceeds 94 °C, the reservoir will belong to the drained zone according to Eq. (7). In the drained zone, both bitumen and water can flow in the connected pores. Therefore, the seepage in this zone can be expressed as where V o is the velocity of bitumen; p o is the pressure of bitumen; K o and K w are the effective permeabilities to bitumen and water, respectively; μ o is the viscosity of bitumen; M o is the mobility of bitumen.
Given that the volume fractions of bitumen and water are b and w , respectively, the absolute permeability after the fusion of bitumen K a can be evaluated by the Kozeny-Carman equation (Wong and Li 2001;Gao et al. 2019b) As shown in Fig. 18, the relative permeabilities to water/bitumen under different water saturations are tested using the Karamay oil sand cores. It can be known that the irreducible water saturation S wi and the residual oil saturation S or are 0.296 and 0.306, respectively.
(25) u w (T) = −0.455 ln T + 2.3485 As shown in Eqs. (29) and (30) and Fig. 18, the relations between the relative permeabilities and water saturation can be depicted by the power-law equations according to the test data (Zhou et al. 2017;Yale et al. 2010).
where K rw and K ro are the water and oil (bitumen) relative permeabilities, respectively; K * rw and nw are two empirical parameters for describing the relation of K rw and S w (K * rw and nw are 0.0795 and 1.3987, respectively); K * ro and no are two empirical parameters for describing the relation of K ro and S w (K * ro and no are 1 and 2.3688, respectively). In this paper, mobility is used to quantitatively describe the flow capacity of the fluid. Using the parameters listed in Table 2, the relation between the mobility and the viscosity of water is shown in Fig. 19a. Water mobility in three geomechanical zones shows significantly different changes with water viscosity. In the undrained zone, water mobility declines from 5.7 to 2.2 mD/mPa s with the increase of water viscosity; and in the drained zone, water mobility also drops from 4 to 3.5 mD/mPa s when the viscosity increases. Due to the fusion of bitumen, water mobility decreases by about 40% in the partially drained zone. This result reveals that water mobility undergoes the change of "decrease-increase-decrease" from the close to the distant outside the steam chamber.
The relation between mobility and the viscosity of bitumen is shown in Fig. 19b. The bitumen in the undrained zone has no flow capability (M o = 0). The bitumen in the drained zone has small mobility because of the low effective permeability to oil (Fig. 18) and the high viscosity (Fig. 13). The oil mobility in the partially drained zone shows a small change due to bitumen fusion. Figure 20a shows the evolutions of the water mobility with oil saturation before and after the fusion of bitumen, respectively. For the oil sands with different bitumen contents (oil saturation), the water mobility changes due to bitumen fusion are different. It can be known that the change in water mobility due to bitumen fusion exhibits a trend of "decrease-increase." The minimum change occurs when the oil saturation is about 57%.
The water mobility after bitumen fusion is taken for example. As shown in Eq. (27), the water mobility M w is determined by K w /u w (T). Here, u w (T) is constant because T is the critical temperature when bitumen fusion.  Fig. 19 a Relations between water mobility and viscosity, and b relations between oil mobility and viscosity in different geomechanical zones Considering Eq. (29), the k w is determined by both the absolute permeability K a , and the relative permeability to water K rw . Equation (28) shows that the absolute permeability increases with the bitumen content ( b ) increasing, while Eq. (29) exhibits that the relative permeability to water decreases with the bitumen content (S o ) increasing. Therefore, the relationship between the water mobility and oil saturation is as shown in Fig. 20a. Figure 20b shows evolutions of the oil mobility with oil saturation before and after the fusion of bitumen, respectively. It can be known that the change in oil mobility due to bitumen fusion exhibits a sharp exponential increase with the increase of oil saturation. For instance, for the bitumen-rich oil sand, the oil mobility change due to bitumen fusion can reach about 0.11 mD/mPa·s; while for the mud-rich oil sand, the oil mobility only changes slightly.

Heat conduction behavior
The heat transfer mechanisms in the three geomechanical zones contain both heat conduction and heat convection related to the flow of water and oil. Firstly, the effective conductivity coefficient in these zones k eff can be expressed as (Kamath and Godbole 1987) where s b and s w are the saturations of oil (bitumen) and water, respectively; k s is the thermal conductivity of solid structure, W/m °C; is the reservoir porosity including bitumen; δ b and δ w are ratios of the thermal conductivity of oil and water, respectively, to the thermal conductivity of the solid structure, and they can be written as From Eq. (31), the evolutions of both δ b and δ w will affect k eff . In this paper, both δ b and δ w are written as functions of their viscosities. The relation between the viscosity and the conductivity of water is (Incropera et al. 2007) AOSTRA (1984) recommended a value of 0.151 W/ (m·K) as the average thermal conductivity of liquid heavy oils and bitumen, and the thermal conductivity of solid bitumen is 0.75 W m −1 K −1 (Ramos-Pallares 2017). Given that the thermal conductivity of bitumen has a linear relation with the viscosity in a range of 10 ~ 10 7 mPa s, there is In this regard, the changes in the effective thermal conductivity with the evolutions of water and bitumen viscosities at a temperature of 10 ~ 100 °C can be calculated according to Eq. (31), just as shown in Fig. 21. It can be seen that the effective thermal conductivity decreases with the increase of both water and bitumen viscosities. When the viscosities are higher, the changes are greater. In reality, the k eff is always around 1 W/(m·K), and the total change is only 0.12 W/(m·K). In this study, viscosity is temperaturedependent. When the temperature rises, the movement of water molecules will be faster. So, in reality, temperature is the fundamental cause that affects the effective thermal conductivity of water as well as reservoir.

Heat convection behavior
It is a common sense that the heat conduction is very slow, while the heat convection related to the flow behaviors of oil and water is relatively fast. Therefore, the heat convection effect is more beneficial to the temperature propagation than the heat conduction in the three zones.
In this paper, the Peclet number is used to quantitatively describe the advantage of heat convection effect over the heat conduction (Rapp 2017). Considering the flow behaviors in the three zones are different, the Peclet number in the three zones should also be separately discussed.
The Peclet number is a dimensionless number which is often used in the study of continuous transmission. The physical meaning of the Peclet number of water (P we ) is the ratio of water convection rate to water diffusion rate. In the three zones, there is both water flow. The Peclet number of water P we is defined as (Rapp 2017) where ρ w is the density of water; C w is the specific heat of water; Δp w is the differential pressure of water. Considering that the changes in the density and specific heat of water at a temperature of 10 ~ 100 °C are slight, the ρ w and C w are regarded as constant values of 995.68 kg/m 3 and 4.178 kJ/ (kg °C), respectively. (36) Considering Eq. (36) and Fig. 19a, the changes in the Peclet number of water at unit pressure (1 MPa) with water viscosity can be calculated, just as shown in Fig. 22a. It can be known that the Peclet number of water in both the drained and undrained zones declines with the increase of water viscosity. There is a sudden change due to the fusion of bitumen, and the Peclet number declines by more than 40%.
The physical meaning of the Peclet number of oil (P oe ) is the ratio of oil convection rate to oil diffusion rate. In the drained and partially drained zones, there is the oil flow. The Peclet number of oil P oe can be expressed as (Rapp 2017) where ρ o is the viscosity of oil; C o is the specific heat of oil; μ o is the viscosity of bitumen; Δp o is the differential pressure of oil. The density and specific heat of bitumen at a temperature range of 10 ~ 100 °C are about 1021 kg/m 3 and 920 J/(kg °C), respectively (Incropera et al. 2007).
Considering Eq. (37) and Fig. 19b, the changes in the Peclet number of oil at unit pressure (1 MPa) with bitumen viscosity can be known, just as exhibited in Fig. 22b. It can be concluded that the Peclet number of oil declines sharply with the increase of oil viscosity, and it shows little change in the partially drained zone. The Peclet number of oil in the undrained zone is zero, because there is no effective permeability to oil in this zone.
For the oil sands possessing different bitumen contents (oil saturations), the change in water/oil Peclet number due to bitumen's fusion is different. Figure 23(a) shows the evolutions of the Peclet number of water with oil saturation before and after the fusion of bitumen, respectively. It can be known that the change in water Peclet number due to bitumen fusion exhibits a trend of "decrease-increase."  The minimum change occurs when the oil saturation is about 57%. Figure 23b shows the evolutions of the Peclet number of oil with oil saturation before and after the fusion of bitumen, respectively. It can be known that the change in oil Peclet number due to bitumen fusion exhibits a sharp exponential increase with the increase of oil saturation. For instance, for the bitumen-rich oil sand, the oil Peclet number change due to bitumen fusion can reach about 0.10; while for the mud-rich oil sand, the oil Peclet number changes slightly.

Implications for the field and simulation work
In this contribution, the viscosity-temperature behavior of Karamay oil sand bitumen at high temperatures was studied. The effects of viscosity-temperature responses on the reservoir deformation, fluid flow, and heat transfer behaviors were quantitatively discussed in three geomechanical zones.
In the oil field, the engineers and operators can use these findings to evaluate the viscosity evolution of Karamay oil sand bitumen with temperatures in the SAGD process. Firstly, in the preheating stage, this study can be used to determine the critical temperature where bitumen can effectively flow (Gao et al. 2019c;Yuan et al. 2013).
To make thermal/hydraulic communications in the interwell zone, the midpoint temperature can be suggested to exceed 94 °C. Secondly, this study can determine the position and area of the three geomechanical zones. In this regard, the reservoir deformation, seepage, and heat transfer behaviors in the three geomechanical zones can be quantitatively assessed. Thirdly, given that the reservoir temperature distributions outside the steam chamber at any time are known, the effective stress applied on structure and the volumetric strain of structure in the SAGD process can be evaluated according to the viscosity-temperature relation. These evaluations can be also used to solve these safety problems such as borehole stability, caprock stability, and surface subsidence. Similarly, the fluid flow behavior evaluated by the viscosity-temperature relation can be employed to analyze these economic problems such as oil production, and the heat transfer behavior can be employed to evaluate the preheating period as well as the development speed of the steam chamber. All in all, this study can provide the engineers with a comprehensive approach about the "viscosity-temperature-deformation-seepageheat transfer" relation in these different geomechanical zones.
For the simulation work, this paper can provide researchers with these quantitative relations used in the coupled thermo-hydro-mechanical analyses. For instance, when the reservoir temperature distributions outside the steam chamber in the current analysis step are known, the viscosity distributions can be determined. Furthermore, all the deformation, fluid flow, and heat transfer parameters used in the next analysis step can be refreshed by this study. Finally, the deformation, seepage, and heat transfer behaviors in the three zones can be quantitatively evaluated considering the evolutions of the positions and areas of the three geomechanical zones. In reality, the role of bitumen in the oil sand reservoir during the SAGD process appears to dominate in the mechanisms of the influence of temperature-viscosity behaviors on the geomechanics in the SAGD process. The solid bitumen functions as the structure, which can support a part of the effective stress; while the fluid bitumen acts just as water, which means an effective oil production and a higher heat transfer efficiency. The different roles of bitumen in different geomechanical zones are worth careful consideration in the simulation work.

Discussion
In this study, one of the focused issues is to predict the viscosity of Karamay oil sand bitumen at these temperatures ranging from initial reservoir temperature to steam chamber temperature, so as to evaluate the fluidity of bitumen that may affect the oil production and bitumen's role in the reservoir. This paper presents a number of data for the viscosity of Karamay oil sand bitumen as a function of temperature, which can be added to the literature on the viscosity behavior of bitumen resources. The modified three-parameter model and the least square method are used to depict these viscosity-temperature relations. However, the Karamay oil sand bitumen is too viscous compared with common oils, so some tricky problems need further studies such as the non-Newtonian behavior of Karamay oil sand bitumen, sample contamination and evaporation or reaction of the sample.
Many bitumens give apparent non-Newtonian behavior at low temperatures, either due to viscoelastic behavior or due to experimental artifacts such as non-uniform sample temperature during shear. Normal practice is to simply define the temperature range of observed linear Newtonian behavior, and exclude the suspect data. By definition, the first check of non-Newtonian behavior is in the analysis of data for strain versus shear stress. Newtonian materials give a linear relationship, characterized by the viscosity. If the relationship is not linear, then the viscosity is not a sufficient description of the rheological behavior. In reality, non-Newtonian behavior cannot be deduced from a viscosity-temperature plot, and in any case the appearance of non-Newtonian behavior rules out the use of viscosity as a parameter. Therefore, the viscosity-temperature relation exhibited in this paper can not reflect the non-Newtonian behavior of Karamay oil sand bitumen. In a further study, the research gap concerning the non-Newtonian behavior of Karamay bitumen will be focused on.
The data for some samples show a discontinuity in the viscosity-temperature plot at 100 °C, which is accounted for by being rich in asphaltenes in this study. However, this may also be possibly due to water-contaminated samples in an open rheometer, although the content of water is very low by eye survey. In further research, the water content of these samples will be further determined by Karl Fisher titration, so as to exclude results from emulsions. An alternative can be used to cool a sample and rerun the experiment. Besides, the evaporation or reaction of the sample is possible to occur in relatively high temperatures. Sample measurements at up to 350 °C would be possible to evaporate a portion of the sample bitumen, depending on the exact design of the instrument. Some bitumens may also undergo slow thermal reactions at temperatures over 300 °C. In the next study, such artifacts will be investigated, for example, by running a sample to 350 °C, then cooling it down and rerunning the measurements to verify that no change in composition has occurred.
In this paper, a lot of empirical relationships from various studies are used. These relations are relatively universal, and the underlying assumptions and field applicability of these relations are similar. The reliability of the evaluation of elastic properties for a porous medium using Hashin's theory is a question worthy of consideration, because these actual complex pore structures are relatively different from these simplified pores in theoretical models. The empirical correlations may fit well with the experimental data, but they cannot explain the mechanism clearly.
Thermo-hydro-mechanical coupled behavior in the process of SAGD has become increasingly important to hydrocarbon operations including the preliminary reservoir stimulation through geo-mechanical methods and subsequent production by the propagation of steam chamber in which temperature-pressure coupling is focused on (e.g., Chalaturnyk and Li 2004;Li and Chalaturnyk 2009;Yale et al. 2010). In this study, some critical model parameters to describe the deformation and fluid flow behaviors based on the viscosity (or temperature) and the oil saturation distributions are discussed, but the effect of the latent heat effect on the thermo-hydro-mechanical coupled behavior in the partially drained zone is not involved in detail. At present, hydro-mechanical coupled behaviors under water and steam injection in Karamay oil sand reservoirs were studied (e.g., Yuan et al. 2013;Lin et al. 2016Lin et al. , 2017Gao and Chen 2020a;Gao et al. 2020a), which can provide a consultation for this study, if the phase change and latent heat effect would be included (it is a thermohydro-mechanical-phase change coupled process). Future work for this challenging problem will be discussed by the numerical method.
This study can be extended to apply to this whole oilfield and other heavy oilfields where the bitumen viscosity exceeds 20,000 mPa·s on the virgin reservoir condition. Li et al. (2004) pointed out that the ultra-heavy reservoir during SAGD process can be divided into three geo-mechanical zones according to bitumen viscosity. In their theory, the reservoir where the bitumen viscosity exceeds 20,000 mPa·s was defined as the undrained zone.

Summary and conclusion
In this contribution, the viscosity of Karamay oil sand bitumen over a temperature range of 20 ~ 350 °C was measured by an advanced rheometer, and the viscosity-temperature behaviors derived from the test data were simulated by a modified three-parameter model and the least square method. The effect of bitumen viscosity on the geomechanics in the SAGD process was discussed in terms of reservoir deformation, fluid flow, and heat transfer behaviors. Major conclusions can be drawn as follows.
The viscosity of Karamay bitumen in areas A, B, and C all drops sharply with the temperature increase. In virgin reservoir conditions, the viscosity of bitumen is dramatically high, which indicates that the Karamay bitumen keeps a solid state and cannot flow at all. The viscosity of Karamay oil drops to be a value less than 10 mPa·s when the temperature is heated to be more than 250 °C, which is favorable for the oil flow and oil production. Because of the high-content asphaltenes of bitumen for wells A3, A4, and A6 in area A, the viscosity-temperature curves show an inflection temperature near 100 °C.
The reservoir deformation, fluid flow, and heat transfer behaviors at varying bitumen viscosity show significant differences in the drained, partially drained, and undrained geomechanical zones. The structure's bulk modulus, effective stress, and volumetric strain in the drained zone are lower than those in the undrained zone; while the oil mobility and Peclet number show the opposite tendency. The changes in the deformation, fluid flow, and heat transfer due to the fusion of bitumen are very considerable for the oil sands with various oil saturation. The changes in the structure's bulk modulus, effective stress, volumetric strain, oil mobility, and Peclet number due to the phase change of bitumen increase with the increase of oil saturation. provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.