An integrated model for productivity prediction of cyclic steam stimulation with horizontal well

As the accurate productivity prediction of cyclic steam stimulation with horizontal well (CSSHW) is essential for its effective utilization in heavy oil reservoirs, a novel integrated model for predicting productivity of CSSHW was proposed in this paper. Firstly, governing equations to predict thermophysical parameters of injected steam in the horizontal well were established. On the basis, the heating model for radius of hot fluid zone and steam zone were derived respectively according to a more realistic temperature distribution of formation. Then, coupled with heating model, the productivity prediction model was established according to the principle of equivalent percolation resistance. Then, comparisons have been made among the new model results, CMG STARS results, and measured data to verify its accuracy. Finally, after the new proposed model is validated, the effects of steam quality and injection rate on heating radius and oil production were analyzed in detail. The results indicate that the result of new model is in good agreement with that of CMG STARS and measured data, verifying the correctness of the model. In addition, the radius of hot fluid zone and steam zone both decreases with horizontal well length because of the heat losses as the steam flows along the horizontal wellbore. Influential factor analysis shows that high steam quality and steam injection rate can improve heating radius and cumulative oil production.

For the above task, some classic researches have been conducted. As for the thermophysical parameter distribution in the wellbore, Ramey 22 firstly established an analytical model for temperature as a function of well depth and injection time based on the assumption of heat transfer in the wellbore being steady state. Satter, 23 Holst and Flock, 24 Gu et al, 25 and Sun et al 26 improved Ramey's model by making the over-all heat transfer coefficient dependent on depth and by considering the effect of condensation and friction losses. Recently, the models for estimating thermophysical properties and analyzing the performance of superheated steam injection in vertical wells and horizontal wells are established. 5,6,[27][28][29][30][31] Moreover, the flow and heat transfer characteristics of superheated steam in concentric dual-tubing wells were analyzed. [32][33][34][35] Sun et al [36][37][38][39][40] studied the flow of supercritical CO 2 coupled with superheated water in horizontal wellbores. As for the formation heating radius, Marx and Langenheim 41 firstly presented an analytical mathematical model for predicting heating area of steam injection with vertical well based on the energy conservation principle. Their model ignored the temperature decrease in the front of heating area. In addition, their model ignored the effect of steam overlay was not taken into account. 42 In recent years, He et al 43 improved the formation heating model by considering the nonisothermal distribution of formation temperature and the effect of steam overlay effect. Zhou et al 31 set up a formation heating model for superheated steam injection with vertical wells, which was further improved by He et al 44 who adopted a more realistic formation temperature distribution of superheated steam injection. As to our knowledge, there is big difference between vertical and horizontal section of wellbore. For example, the gas injection tubing usually is surrounded by insulation, casing, and cementing in the vertical wellbore; nevertheless, the gas injection tubing usually is surrounded by cementing or is bare to formation. 26,27,44 Also, there are great differences in the heating processes between vertical wells and horizontal wells. 28,32 For instance, for vertical wells, a part of heat loses to the boundary and the remaining part is used for heating the formation when the steam reaches the bottomhole. However, for the horizontal wells, the heating process is divided into several processes based on whether the frontier of heat area reaches boundary or not. 43,45 More importantly, as the steam flows along the horizontal wellbore, the steam temperature, pressure, and quality always change with horizontal well length. 46 Therefore, those above heating models of vertical wells cannot be used to directly predict the heating area of horizontal wells. 29,35,[47][48][49] Wu et al 50 established an analytical model for steam chamber by assuming that it has the shape of a cylinder and the thickness of pay zone is larger than the chamber diameter. However, it ignores the heat lost to the boundary when the frontier of heat area reaches boundary. Hamid et al 51 established a two-dimensional transient heat conduction model to address the formation heating process by CSSHW with consideration of the effect of partial penetration of the horizontal well in radial coordinate. He et al 45 built a mathematical model for heat radius of cyclic superheated steam stimulation with horizontal wellbore, and this method is adopted to calculate the heating radius of CSSHW in this paper. As for the production prediction of CSS and CSSHW, numerical thermal simulation method is usually adopted. However, as executing a thermal reservoir simulator is time-consuming and the use of it may be limited when geological features and fluid properties are not available, analytical productivity prediction model will be an attractive option for planning heavy oil field development. Boberg 50 presented an analytical model to predict the oil production of CSSHW. However, it does not consider the nonisothermal distribution characteristics of heating zone and does not consider the heat lost to the boundary when the frontier of heat area reaches boundary.
The main objective of this work was to establish an integrated model to accurately predict productivity of CSSHW in thin heavy oil reservoirs by taking the nonisothermal distribution characteristics of heating zone and the situation of the frontier of heating zone reaching boundary into account. In this article, governing equations to predict thermophysical parameters of injected steam in the horizontal wellbore were built firstly. Then, the equations for radius of hot fluid zone and steam zone were derived respectively according to a more realistic temperature distribution of formation. On the basis, the productivity prediction model was established according to the principle of equivalent percolation resistance. Finally, after the new proposed model is validated by comparison with the results of STARS, the effects of influence factors on heating area are analyzed in detail.

HE Et al.
The horizontal wellbore is evenly divided into N segments. Taking the i-th segment as an example, the mass balance equation can be expressed as 25,[55][56][57] where w i is the mass flow rate at the outlet of the i-th segment; ρ i is the average steam density of the i-th segment; I i is the volumetric outflow rate of steam at the i-th segment.
According to Gu et al. 25 , the pressure gradient of steam in the horizontal wellbore can be given as 39,58 where dp i /dl is the pressure gradient at the i-th segment, p i is the steam pressure at the outlet of the i-th segment; f w and f p are the friction factor of pipe flow and perforation roughness, respectively; v ai and v sgi are the average velocity of steam and superficial gas velocity at the i-th segment, respectively; D is the inside diameter of casing.
As for saturated steam, there is one-to-one correspondence between steam pressure and steam temperature. An empirical correlation is recommended for this relationship as follows 59 where T si is the steam temperature at the outlet of the i-th segment.
Based on the law of energy conversation, 33,34,60 it obtains where dQ i /dl is the heat flow rate from the steam to the surrounding formation at the i-th segment, whose detailed expression can be found in Gu et al 46 ; h i is the specific enthalpy of steam at the outlet of the i-th segment. Based on thermodynamic principles, the enthalpy gradient, the second term in Equation (4), can be written as 36,61 where h si and h wi are the specific enthalpy of dry steam and the specific enthalpy of saturated water at the outlet of the i-th segment, respectively; x i is the steam quality at the outlet of the i-th segment.
Incorporating Equations (4) and (5), the steam quality distribution along the horizontal wellbore can be expressed as x 0 is the steam quality at the heel position of the horizontal wellbore; l i is the wellbore length at the outlet of the i-th segment.
Since the steam pressure and steam quality depend on each other, they should be solved by iterative method. The mass flow rate, steam temperature, steam pressure, and steam quality in each segment can be calculated with iterative method according to Equations (1), (2), (3), and (6) until the toe of the horizontal wellbore is reached.

| Formation heating model
After the steam inside the horizontal wellbore enters the formation, the formation is heated by hot steam. In order to accurately predict the formation heating area resulting from the injected steam, several basic assumptions are made as follows: 1. The formation thickness is relatively small, and therefore, the steam override effect is not taken into account, as shown in Figure 1A. 2. The hemispheric heating zones at the heel position and the toe position of the horizontal wellbore are ignored. 3. The formation is divided into three zones 45 : steam zone, hot fluid zone, and unheated zone. The temperature in the steam zone equals to the injected steam temperature, while the temperature in the hot fluid zone equals to the arithmetic mean value of the steam temperature and the initial reservoir temperature (as shown in Figure 1B).

4.
Heat transfers radially at the first stage; accordingly, the heated zone of each horizontal wellbore segment has the shape of cylinder. After the frontier of hot fluid zone and steam fluid zone arriving at the boundary, heat transfers linearly in the formation.
When the hot fluid zone frontier does not reach the boundary, all the heat carried by the injected steam is used to heat the formation. 41 Therefore, it obtains where H si is the steam zone heat injection rate of the i-th segment, H si = I i ρ i x i L vi , L vi is the steam latent heat of the i-th segment; Δl is the length of each segment; T r is the initial reservoir temperature; t is the injection time; M R is the heat capacity of reservoir.
Hence, the radius of steam zone of the i-th segment is where r si is the steam zone radius of the i-th segment. (1) Similarly, based on the energy conservation principle in the hot fluid zone, the hot fluid zone radius can be expressed as where r hi is the hot fluid zone radius of the i-th segment; H wi is the hot fluid zone heat injection rate of the i-th segment, H wi = I i i h wsi − h wri , h wsi is the water-specific enthalpy of the i-th segment at the steam temperature and h wri is the water-specific enthalpy of the i-th segment at the reservoir temperature; T hi is the hot fluid zone temperature of the i-th segment.
Based on Equations (8) and (9), the time that the hot fluid zone frontier is just reaching the boundary and the time that the steam zone frontier is just reaching the boundary can be expressed as 45 where t 1i is the time that the hot fluid zone frontier of the i-th segment reaches the boundary; t 2i is the time that the steam zone frontier of the i-th segment reaches the boundary; d is the half thickness of the formation.
When the hot fluid zone frontier reached the boundary but the steam zone frontier did not, that is, t 1i ≤ t < t 2i , the steam zone radius still can be obtained by Equation (8). Based on the energy conservation principle, a heat balance of hot fluid zone yields 41,43,45 where λ s is the thermal conduction coefficient of boundary; α s is the thermal diffusivity of boundary; δ is the instant at which the boundary becomes exposed to the hot fluid.
After Laplace transformation and inverse Laplace transformation of Equation (12), the hot fluid zone radius of the i-th segment can be obtained Similarly, when the steam zone frontier reached the boundary, the steam zone radius and the hot fluid zone radius of the i-th segment can be expressed as

| Productivity prediction model
After a few days of soaking to heat the formation, the horizontal well will open for production. The basic assumptions for the productivity prediction model include the following: (a) The bottom-hole flow pressure remains constant during production and the change of pressure along the horizontal wellbore is not considered; (b) each segment is equivalent to a vertical well and the flow conforms to steady flow at small interval. The seepage process in the horizontal segment is the coupling of linear flow and radial flow, 62,63 as shown in Figure 2.
According to the principle of equivalent percolation resistance and the above formation heating model, the oil production rate and the water production rate can be expressed as 61 where Q o and Q w are oil production rate and water production rate, respectively; J oi and J wi are the oil productivity index and water productivity index of the i-th segment, respectively; p i is the average pressure of the i-th segment; p wf is the bottomhole flow pressure; R (o) 1i and R (o) 2i are the linear flow resistance and radial flow resistance of oil phase, respectively; R (w) 1i and R (w) 2i are the linear flow resistance and radial flow resistance of water phase, respectively.
When the hot fluid zone frontier does not reach the boundary, R (o) 1i and R (o) 2i can be expressed as; 63,64 When the hot fluid zone frontier reached the boundary but the steam zone frontier did not, R (o) 1i and R (o) 2i can be expressed as When the steam zone frontier reached the boundary, R (o) 1i and R (o) 2i can be expressed as where μ os , μ oh , and μ oc are the oil phase viscosity of steam zone, hot fluid zone, and unheated zone, respectively; K ros , K roh , and K roc are the oil phase relative permeability of steam zone, hot fluid zone, and unheated zone, respectively; r w is the wellbore radius; K is the absolute permeability of formation; r e is the half drainage radius. Similarly, replacing μ os , μ oh , and μ oc with μ ws , μ wh , and μ wc and replacing K ros , K roh , and K roc with K rws , K rwh , and K rwc in Equations (18)(19)(20), the linear flow resistance and radial flow resistance of water phase can be obtained. μ ws , μ wh , and μ wc are the water phase viscosity of steam zone, hot fluid zone, and unheated zone, respectively; K rws , K rwh , and K rwc are the water phase relative permeability of steam zone, hot fluid zone, and unheated zone, respectively.
Since the viscosity of oil reduces rapidly with temperature, the acquisition of reservoir temperature is essential to the calculation of oil viscosity. During the soaking stage, (14)  the temperature of steam zone and hot fluid zone decreases gradually due to radial and vertical thermal conductivity. The average temperature of steam zone and hot fluid zone at the end of soaking can be expressed as 52,54 where T asi and T ahi are the average temperature of steam zone and hot fluid zone at the end of soaking, respectively; V rsi and V rhi are the radial heat loss influence factor of steam zone and hot fluid zone, respectively; V zsi and V zhi are the vertical heat loss influence factor of steam zone and hot fluid zone, respectively.
During the production stage, the produced liquid also carries some heat, which further reduces the temperature of steam zone and hot fluid zone. 54,65 The average temperature of steam zone and hot fluid zone at the production stage can be expressed as where T asi and T ahi are the average temperature of steam zone and hot fluid zone at the production stage, respectively; δ si and δ hi are the correction coefficient of steam zone and hot fluid zone due to production, respectively. V rsi , V rhi , V zsi , V zhi , δ si , and δ hi can be calculated by the method in reference. 54 The injection of steam will result in the average reservoir pressure at the end of soaking being higher than the original formation pressure. 52,54 According to the volume balance principle, the average reservoir pressure at the end of soaking can be expressed as where P r is the initial reservoir pressure; p 1i is the average reservoir pressure of the i-th segment at the end of soaking; B w and B o are the volume factor of water and oil, respectively; C e is the total compression coefficient; β e the total thermal expansion coefficient; N i is the oil and water in place of the i-th segment; N osi and N ohi are the oil in place of steam zone and hot fluid zone of the i-th segment, respectively; G i is the cumulative injected steam volume of the i-th segment at reservoir condition. Similarly, after the well is open for production, the average formation pressure will decrease continuously due to the production of oil and water. The average reservoir pressure at the production stage can be expressed as where p i is the average reservoir pressure of the i-th segment at the production stage; N wi is the cumulative water output of the i-th segment; N oi is the cumulative oil output of the i-th segment.
According to the mass conservation principle, the oil saturation at the end of soaking and at the production stage can be obtained 4 where S i w is the water saturation of the i-th segment; S wi is the initial water saturation of the i-th segment; φ is the porosity; V si and V hi are the steam zone and hot fluid zone volume of the i-th segment, respectively.
At the end of each huff and puff cycle, some residual heat will exist in the reservoir. The residual heat from the previous cycle should be added to the next cycle. The (25)  residual heat of steam zone and hot fluid zone can be expressed as 52,54 where E rsi and E rwi are the residual heat of steam zone and hot fluid zone of the i-th segment, respectively.

| CALCULATION FLOWCHART FOR THE MODEL
Based on the established model, the calculation flowchart is shown in Figure 3.

| Model verification
To verify the correctness of the proposed integrated model, a horizontal steam injector well H-470 in KMK oilfield, Aktyubinsk, northwest of Kazakhstan, is used as an example to calculate the heating radius and to predict the oil production rate. The basic parameters of formation and steam injection parameters of well H-470 are listed in Table 1. The relative permeability curve and the viscosity-temperature curve can be obtained in reference. 54 We compare the result of the new integrated model with that of CMG STARS and the measured data as shown in Figure 4. It is clearly observed from Figure 4A that the heating radius predicted by the new proposed model is in good agreement with the CMG STARS simulation result, which supports the reliability of the new model. It also shows that there are some differences between them. The reason for the difference is that we made some simplification for the sake of quick calculation. We simplify the formation temperature of hot fluid zone to the arithmetic mean value of the steam temperature and the initial reservoir temperature, but in reality the temperature of hot fluid zone tends to be exponentially. Figure 4B shows that the oil production rate predicted by the new model and by CMG STARS and the actual measured data are consistent. Therefore, this new model is proved to be accurate and efficient. It should be pointed out that the result of Wu et al's model is larger than that of the new model and measured date for the reason that it ignores the heat lost to the boundary when the frontier of heat area reaches boundary. Figure 5 shows the effect of steam quality at the horizontal wells' heel position (x 0 ) on the profiles of heating radius in the horizontal wellbore. From Figure 5A,B, it is observed that for a given steam quality, the hot fluid zone radius and the steam zone radius both decrease with horizontal well length. This is because the steam quality and the steamspecific enthalpy decreases with horizontal well length as the steam flows along the horizontal wellbore on account of the heat and mass transfer into the reservoir, thereby decreasing of the heating radius along the horizontal wellbore. Meanwhile, Figures 5 and 6 indicate that the larger the steam quality is, the greater the heating radius and cumulative oil production. For instance, when the steam quality equals to 0.4, the radius of hot fluid zone and steam zone at the well length of 100 m are 4 m and 1.8 m, respectively, and the cumulative oil production at the end of the first cycle is 1352 m 3 . However, when the steam quality increases to 1, the radius of hot fluid zone and steam zone will increase to 4.7 and 2.9 m, respectively, and the cumulative oil production at the end of the first cycle will increase to 1544 m 3 . Consequently, to ensure the heating effect of injected steam on the formation, the steam quality should not be too small. Figure 7 shows the effect of injection rate at the horizontal wells' heel position (w 0 ) on the profiles of heating radius in the horizontal wellbore. From Figure 7A,B, it is clearly found that the larger the injection rate is, the greater the radius of hot fluid zone and steam zone. The reason is that the larger the injection rate at the horizontal wells' heel position, the larger mass of steam injected into the formation per unit of time result in greater heating area. Figure 8 indicates that the larger the steam injection rate is, the greater the heating radius and cumulative oil production. Therefore, enhancing the injection rate may be a good choice to increase the heating radius during steam injection.

| CONCLUSIONS
This paper proposes an integrated model for productivity of CSSHW in thin reservoirs by coupling with the wellbore model and formation heating model. The following conclusions can be summarized from the results of this work:

1.
Since the proposed integrated model considered the change in steam's thermophysical parameters with horizontal well length, the nonisothermal distribution characteristics of heating zone and the situation of the frontier of heating zone reaching boundary, the results predicted by this new model are in good agreement with the CMG STARS results and measured data, which verify the correctness of the model.

2.
The radius of hot fluid zone and steam zone both decreases with horizontal well length on account of the heat losses as the steam flows along the horizontal wellbore. 3. High steam quality and steam injection rate can improve heating area and cumulative oil production.