© H. Liu et al., published by IFP Energies nouvelles, 2017

Nomenclature

A : Area of heated zone, m²

(t): Solvent diffusion area at time t, m²

c_P : Heat capacity at constant pressure within the reservoir, J/(kg °C)

c_cap : Heat capacity of understrata and overburden rock

C_si : Solvent solubility at the steam front

C_s : Solvent concentration at ε

: Average solvent concentration in the solvent penetration depth

D_t : Effective diffusion coefficient of solvent at time t, m²/d

D : Solvent molecular diffusion coefficient, m²/s

dr/dt : Radial velocity of steam front, m/s

dx/dt : Longitudinal velocity of steam front, m/s

dδ_c /dt : Solvent diffusion rate, m/s

f_w : Water fractional flow at specific water saturation

: Slope of tangent of the water fractional flow curve at specific water saturation

h : Reservoir thickness, m

h_w : Enthalpy of saturated water, kJ/kg

H_o : Heat injection rate per unit well length, kJ/md

L_f : Horizontal displacement of cold water front in the second stage, m

L_e : Effective operating thickness of solvent, m

L_v : Latent heat of vaporization of steam, kJ/kg

l : Length of well, m

q_hloss : Heat loss rate to understrata and overburden per area, kJ/m²d

q_heff : Required heat rate associated with the steam zone expansion per unit well length, kJ/md

S_or : Residual oil saturation

S_wc : Connate water saturation

S_wt : Water saturation at the end of operating thickness

S_wp : Water saturation at the production well

T : Temperature, °C

T_r : Initial reservoir temperature, °C

t_c : Time when the first stage ends, d

r : Radius of steam zone area, m

: Mass rate of cold steam condensate per unit well length, kg/md

X : Steam quality

x : Horizontal displacement of steam front, m

Greek letters

α_r : Longitudinal dispersity factor

δ_c : Solvent penetration depth, m

ε : Distance beyond the edge of the steam front, m

ϕ : Reservoir porosity

λ_cap : Thermal conduction coefficient of understrata and overburden rock, kJ/md °C

μ : Viscosity, mPa s

ρ : Density, kg/m³

Γ(·): Gamma function

Subscripts

o: Oil

s: Steam

w: Hot water

r: Reservoir rock matrix

cap: Understrata and overburden rock

Introduction

Heavy oil resources have an important role in crude oil reserve replacement to meet the world’s future energy demands [1]. However, a considerable portion of the heavy oil resides in thin formations and in some countries more than 80% lies in reservoirs with a less than a 10 m pay zone; this has posed a major challenge for efficient development [2]. To date, the predominant in situ recovery methods for reservoirs with thin pay zone are steam injection by horizontal wells. Among these, the steam-flooding strategy has been regarded as a significantly superior technology compared to Steam-Assisted Gravity Drainage (SAGD) operations in terms of energy efficiency [3] and has been applied widely in thin heavy oil reservoirs in China [4-6]. However, significant heat losses to the understrata and overburden can render the steam-flooding process uneconomic because of the small thickness of the reservoir [7]. Solvent-Enhanced Steam Flooding (SESF) requires substantially less steam usage and net injected energy and thus is considered as a highly effective method for improving heavy oil recovery in thin reservoirs [7, 8].

Solvent properties and injection strategies are important factors that control the performance of SESF and thus determine the economic viability of the process. Therefore, it is necessary to fully investigate the effects of these factors on production performance before implementing in the field. By conducting 63 experiments on a three-dimensional elemental physical model, Redford [9] first demonstrated the effects of solvent properties on the production performance and indicated that naphtha with steam could significantly reduce oil viscosity. Besides, bitumen recovery could be improved further by adding CO₂ or ethane because they can provide more drive energy for the naphtha to enter the reservoir. Nasr and Pierce [10] conducted a series of experiments on steam-CO₂ injection strategies in a scaled model and revealed that steam-CO₂ continuous injection results in improved performance compared to steam-CO₂ sequential injection. Moreover, oil recovery rates are typically accelerated and improved during the initial stage of the process. Although the effects of the solvent properties and injection strategies of several solvents on steam flooding can be identified by well-designed experiments, the number of solvents tested is limited whereas the experiment requirements are large. Therefore, the experimental data is difficult to analyse.

Based on laboratory experiments, several attempts using numerical simulation have been undertaken to study the SESF process. Shutler and Boberg [11] first used numerical simulations to delineate the recovery mechanism of the solvent-steam process. The results indicate that medium-alkane solvents provide the greatest increase in total oil production, whereas heavy-alkanes solvents do not improve recovery. They also proposed the corollary that the production of the process is primarily determined by the placement of the solvent in the reservoir, which in turn is controlled by the steam movement and solvent volatility. Zhao et al. [8] reported reservoir simulation results on the application of SESF for a bitumen reservoir with a net pay of 4 m. The results present optimum steam and steam-solvent flooding injection strategies by comparing the cumulative Energy injected to produced Oil Ratio (cEOR) between different cases. Hence, numerical simulation was proven to be a powerful method to calculate SESF production. However, the calculation of the numerical simulation is not only time-consuming but also relies heavily on a large amount of experimental data, which limits the application of the numerical simulation method.

For quick calculation and an easy method for the solvent-steam process, several analytical and semi-analytical methods were developed. Sharma and Gates [12] considered mathematical models to calculate the length scales of the mass and heat transfer beyond the edge of a co-injection chamber. Gupta and Gittins [13] developed a semi-analytical approach to estimate oil drainage in the co-injection of steam and solvent. The solvent distribution beyond the steam chamber calculated by their models can be an indicator of different solvents without introducing excessive experimental data. However, their mathematical models are based on Butler’s SAGD model [14, 15] for oil drainage and do not consider the blocking effects of water on the solvent diffusion, which makes the models not applicable for SESF process.

To the best of our knowledge, applying a mathematical model to understand the SESF process for heavy oil reservoirs with thickness less than 10 m has never been undertaken. In this study, to predict SESF production performance and study the effects of solvent properties and injection strategies on the SESF process, a new semi-analytical model is proposed based on mass and energy balance principles, diffusion theory, and the theory of Buckley-Leverett. First, the mathematical model is solved in a semi-analytical method. Next, comparisons are made between the calculated model and numerical simulations to verify the proposed model. Finally, with the validated model, the effects of the solvent properties and injection strategies on the production performance of SESF are studied in detail.

1 Mathematical Model

A simplified diagram of the SESF process is presented in Figures 1 and 2. To decouple the primary mechanisms from the other complexities that accompany the SESF process, important assumptions and simplifications are provided below:

1. based on the position of the steam front, the entire production process is divided into three stages: before the steam front reaches the upper and lower boundary (steam front rising stage) (Fig. 1b), after the steam front reaches the boundary without steam condensate reaching the production well (steam front spreading stage) (Fig. 2b), after the cold steam condensate reaches the production well (steam condensate-producing stage). The pressure inside the steam zone is constant during all stages [16, 17];

2. an interface model proposed by van Lookeren [18] has been used widely to describe the shape of the steam front without considering heat loss. It is a reasonable model for the steam-flooding process in thick reservoirs; however, it is not applicable for thin reservoirs with considerably more severe heat losses to the understrata and overburden. Thus, according to the reservoir numerical simulation results, the geometries of the steam and water fronts are assumed to be a cosine curve and circular shaped, respectively;

3. all mobile oil is assumed to be displaced by the steam and hot water; thus, only residual oil is left in the steam and hot water zone;

4. the gas-phase solvent can move freely in the hot water zone owing to the high temperature and low oil saturation [17, 19], whereas the cold steam condensate may block the diffusion of solvent into the heavy oil. By investigating VAPor-EXtraction (VAPEX) experiments, Das and Butler [20] proposed that the contact area of the solvent and crude oil was directly related to the rate of the mass transfer, as indicated in Figure 3a. Thus, it is assumed that the diffusion coefficient is proportional to the contact area in a limited distance beyond the edge of the steam front, as illustrated in Figure 3b. The limited distance is defined as the effective operating thickness of the solvent, which is determined by the volatility of solvent;

   Figure 1 Illustration of steam front rising stage. a) Three-dimensional diagram of the stage; b) two-dimensional diagram of the stage.

   Figure 2 Illustration of steam front spreading stage. a) Three-dimensional diagram of the stage; b) two-dimensional diagram of the stage.

   Figure 3 Schematic of contact area change at interface. a) Das and Butler’s perspective; b) simplification based on the Das and Butler’s perspective.

5. medium solvents (C₆-C₁₂) are considered to provide the greatest increase in total oil production because they can both travel with the steam in a vaporized phase and condense in the cooler regions of the reservoir and diffuse in the crude oil [21]. Light solvents (CO₂, C₁-C₅), which have considerably higher volatility, are believed to transport further into the reservoir and thus increase the operating thickness when co-injecting with medium solvents. Therefore, the solvents used in this model are assumed to be medium solvents and a combination of medium and light solvents;

6. local equilibrium is assumed along the steam front such that the equilibrium state inside the steam zone is not affected by the phase behaviour at the steam front [22]. Other assumptions, such as no mutual solubility between the water and hydrocarbon components and no variation in the rock-fluid properties along the horizontal section of the well are identical to those proposed by Keshavarz et al. [23].

The mathematical model is divided into three parts based on the three stages (steam front rising stage, steam front spreading stage, steam condensate-producing stage). These are introduced in detail in the following.

1.1 Steam Front Rising Stage

The process of the steam front rising stage is short because of the small thickness of the reservoir, which makes the generating of steam condensate negligible. Therefore, the velocity of the steam front can be calculated by the energy conservation equation:$$\begin{array}{l}{H}_{\mathrm{o}}=\left[{\rho }_{\mathrm{r}}{c}_{\mathrm{Pr}}\frac{(1-\varphi )}{\varphi }+{\rho }_{\mathrm{o}}{S}_{\mathrm{or}}{c}_{\mathrm{Po}}+{\rho }_{\mathrm{s}}\Delta S_{\mathrm{o}}{c}_{\mathrm{Ps}}+{\rho }_{\mathrm{w}}{S}_{\mathrm{wc}}{c}_{\mathrm{Pw}}\right]\varphi \mathrm{dl}(T_{\mathrm{s}}-T_{\mathrm{r}})\frac{\mathrm{d}A}{\mathrm{d}t}\\ =2\pi \left[{\rho }_{\mathrm{r}}{c}_{\mathrm{Pr}}\frac{(1-\varphi )}{\varphi }+{\rho }_{\mathrm{o}}{S}_{\mathrm{or}}{c}_{\mathrm{Po}}+{\rho }_{\mathrm{s}}\Delta S_{\mathrm{o}}{c}_{\mathrm{Ps}}+{\rho }_{\mathrm{w}}{S}_{\mathrm{wc}}{c}_{\mathrm{Pw}}\right]\varphi \mathrm{dl}(T_{\mathrm{s}}-T_{\mathrm{r}})\mathrm{rU}\end{array}$$(1) $$\Delta S_{\mathrm{o}}=1-S_{\mathrm{or}}-S_{\mathrm{wc}}$$(2)where H_o is the heat injection rate per unit well length; A is the area of the steam front; ρ_r, ρ_o, ρ_s, and ρ_w are the densities of the rock, crude oil, steam, and water, respectively; c_Pr, c_Po, c_Ps, and c_Pw are the heat capacities of the rock, crude oil, steam, and water, respectively; ϕ is the reservoir porosity; S_or is the residual oil saturation; S_wc is the connate water saturation; dl is the unit length of well; T_s is the steam temperature; T_r is the initial reservoir temperature; r is the radius of steam zone area; and U is the velocity of steam front.

The unsteady-state solvent concentration distribution along the steam front can be described by Equation (3) [24]:$$\frac{\mathrm{\partial }}{\mathrm{\partial }\epsilon }\left[(D+{\alpha }_{\mathrm{L}}U)\frac{\mathrm{\partial }{C}_{\mathrm{s}}}{\mathrm{\partial }\epsilon }\right]+U\frac{\mathrm{\partial }{C}_{\mathrm{s}}}{\mathrm{\partial }\epsilon }=\frac{\mathrm{\partial }{C}_{\mathrm{s}}}{\mathrm{\partial }t}$$(3)

The corresponding boundary conditions are:$$\{\begin{array}{cc}{C}_{\mathrm{s}}=0,& t_{\mathrm{s}}=0\\ C_{\mathrm{s}}=C_{\mathrm{si}},& \epsilon =0\\ C_{\mathrm{s}}=0,& \epsilon \to \mathrm{\infty }\end{array}$$(4)where ε is the distance beyond the edge of the steam front; C_si is the solvent solubility at the steam front; C_s is the solvent concentration at ε; D is the solvent molecular diffusion coefficient; α_L is the longitudinal dispersity factor; and t_s is the solvent injection time.

Owing to the nonlinearity of Equation (3), the Heat Integral Method (HIM) is used to determine the concentration profile. This method is a powerful approximation method for solving a variety of diffusion problems [25, 26]. Rabiei Faradonbeh et al. [24] realized that the exponential form of the HIM method provides more accurate predictions of the concentration profile than the polynomial form; therefore, the exponential form of the concentration distribution is used for the calculations of the fluid mixture and oil flow rate in this study.

The concentration distribution in exponential form ahead of the interface can be expressed as:$$C_{\mathrm{s}}=C_{\mathrm{si}}{e}^{-\frac{\epsilon }{{\delta }_c}}$$(5)where δ_c is the solvent penetration depth. Transforming Equation (3) into integral form in the domain of the solvent penetration depth and substituting Equation (5) for C_s, Equation (3) can be simplified into:$$\left(D+{\alpha }_{\mathrm{L}}U\right)\frac{1}{{\delta }_c}-U=\frac{\mathrm{\partial }{\delta }_c}{\mathrm{\partial }t}, \hspace{1em}\left(t=0,\hspace{0.5em}{\delta }_c=0\right)$$(6)

According to the Butler’s theory, all mobile oil is displaced by the steam and only residual oil is left in the steam zone. Therefore, based on the mass conservation equation and assumption (6), the oil production rate in this stage can be expressed as:$$q_{\mathrm{oi}}={\rho }_{\mathrm{o}}\varphi \Delta S_{\mathrm{o}}\mathrm{dl}\frac{\mathrm{d}r}{\mathrm{d}t}+{\rho }_{\mathrm{o}}\varphi \left(1-S_{\mathrm{wc}}\right)\mathrm{dl}{\overline{C}}_{\mathrm{s}}\frac{\mathrm{d}{A}_{{\delta }_c\left(t\right)}}{\mathrm{d}t}$$(7) $${\overline{C}}_{\mathrm{s}}=\frac{{\int }_0^{{\delta }_c}{C}_{\mathrm{si}}{e}^{-\frac{\epsilon }{{\delta }_c}}\mathrm{d}\epsilon }{{\delta }_c}$$(8)where is the solvent diffusion area at time t and is the average solvent concentration in the solvent penetration depth.

As indicated in Figures 4a and 4b, the steam front rising stage can be divided into two secondary stages according to the position of the solvent boundary. They can be described by Equations (9) and (10), respectively.$$\frac{\mathrm{d}{A}_{{\delta }_c\left(t\right)}}{\mathrm{d}t}=\frac{\mathrm{d}\left\{\frac{\pi }{2}\left[{\left(r+\frac{\partial {\delta }_c}{\partial t}\right)}^2-r^2\right]\right\}}{\mathrm{d}t},\begin{array}{l}\\ \end{array}\left(r+\frac{\partial {\delta }_c}{\partial t}\le \frac{h}{2}\right)$$(9) $$\frac{\mathrm{d}{A}_{{\delta }_c\left(t\right)}}{\mathrm{d}t}=\frac{\mathrm{d}\left\{\left[h\sqrt{{\left(r+\frac{\partial {\delta }_c}{\partial t}\right)}^2-{\left(\frac{h}{2}\right)}^2}+{\left(r+\frac{\partial {\delta }_c}{\partial t}\right)}^2\left(\pi -2\arccos \frac{h}{2\left(r+\frac{\partial {\delta }_c}{\partial t}\right)}\right)\right]-\frac{\pi r^2}{2}\right\}}{\mathrm{d}t},\begin{array}{l}\\ \end{array}\left(r+\frac{\partial {\delta }_c}{\partial t}>\frac{h}{2}\right)$$(10)

1.2 Steam Front Spreading Stage

The steam front spreading stage is schematically presented in Figure 2. In this stage, the steam is already touching the understrata and overburden; therefore, the heat injected is divided into two parts based on the energy conservation equation:$$H_{\mathrm{o}}=q_{\mathrm{heff}}+q_{\mathrm{hloss}}$$(11)where q_hloss is the heat loss rate to the understrata and overburden per area and q_heff is the required heat rate associated with the steam zone expansion per unit well length.

The heat consumption rate used to expand the steam zone is expressed by the energy increase of the area. As described in assumption (2) and Figure 2b, the area of the steam zone is twice that of the hot water zone. Therefore, the required heat rate associated with the steam expansion per unit well length in this stage can be expressed as:$$\begin{array}{l}{q}_{\mathrm{heff}}=\frac{2}{3}\left[{\rho }_{\mathrm{r}}{c}_{\mathrm{Pr}}\frac{(1-\varphi )}{\varphi }+{\rho }_{\mathrm{o}}{S}_{\mathrm{or}}{c}_{\mathrm{Po}}+{\rho }_{\mathrm{s}}\Delta S_{\mathrm{o}}{c}_{\mathrm{Ps}}+{\rho }_{\mathrm{w}}{S}_{\mathrm{wc}}{c}_{\mathrm{Pw}}\right]\left(T_{\mathrm{s}}-T_{\mathrm{r}}\right)\varphi {\rho }_{\mathrm{o}}\Delta S_{\mathrm{o}}\frac{\mathrm{d}A}{\mathrm{d}t}\\ \begin{array}{ll}& \end{array}+\frac{1}{3}\left[{\rho }_{\mathrm{r}}{c}_{\mathrm{Pr}}\frac{(1-\varphi )}{\varphi }+{\rho }_{\mathrm{o}}{S}_{\mathrm{or}}{c}_{\mathrm{P}\mathrm{o}}+{\rho }_{\mathrm{s}}\Delta S_{\mathrm{o}}{c}_{\mathrm{Pw}}\right]\left(T_{\mathrm{w}}-T_{\mathrm{r}}\right)\varphi {\rho }_{\mathrm{o}}\Delta S_{\mathrm{o}}\frac{\mathrm{d}A}{\mathrm{d}t}\end{array}$$(12)

According to assumption (2), the geometry of the water front is assumed to be a diameter-constant circular shape; thus, the change of the heated area can be transformed to the change of distance:$$\frac{\mathrm{d}A}{\mathrm{d}t}=h\frac{\mathrm{d}x}{\mathrm{d}t}$$(13)

Carslaw and Jaeger [27] built and solved the heat loss model of a semi-infinite plate with a constant temperature boundary. As the hot water mainly contacts the understrata and overburden, the heat loss rate is written as:$$q_{\mathrm{hloss}}=2{\int }_{t_c}^t(T_{\mathrm{w}}-T_{\mathrm{r}})\sqrt{\frac{{\lambda }_{\mathrm{cap}}{\rho }_{\mathrm{cap}}{c}_{\mathrm{cap}}}{\pi (t-\tau )}}\frac{\mathrm{d}x}{\mathrm{d}\tau }\mathrm{d}\tau$$(14)where T_s and T_w are the temperature of the steam and hot water, respectively; ρ_cap is the density of the understrata and overburden rock; λ_cap is the thermal conduction coefficient of the understrata and overburden rock; c_cap is the heat capacity of the understrata and overburden rock; x is the horizontal displacement of the steam front; and t_c is the time when the first stage terminates.

Rabiei Faradonbeh et al. [28] also employed HIM to determine the temperature profile. Therefore, the temperature of the hot water can be expressed as:$$T_{\mathrm{w}}=\frac{{\int }_0^{{\delta }_{\mathrm{cT}}}\left(T_{\mathrm{s}}-T_{\mathrm{r}}\right)e^{-\frac{\epsilon }{{\delta }_{\mathrm{cT}}}}\mathrm{d}\epsilon }{{\delta }_{\mathrm{cT}}}+T_{\mathrm{r}}$$(15)

Based on the mass conservation law, the generating rate of the cold steam condensate per unit well length can be given as:$$v_{m_{\mathrm{c}}}=\frac{H_{\mathrm{o}}}{X{L}_{\mathrm{v}}+h_{\mathrm{w}}}-{\rho }_{\mathrm{s}}\mathrm{dl}\frac{\mathrm{d}{A}_{\mathrm{s}}}{\mathrm{d}t}-{\rho }_{\mathrm{w}}\mathrm{dl}\frac{\mathrm{d}{A}_{\mathrm{w}}}{\mathrm{d}t}$$(16)where is the mass rate of cold steam condensate per unit well length; X is the steam quality; L_v is the latent heat of steam vaporization; h_w is the enthalpy of saturated water; and A_s and A_w are the areas of the steam zone and hot water zone, respectively.

The areas of the steam zone and hot water zone in Equation (16) can be rewritten as Equation (17) and Equation (18) because of the cosine curve and circular shape of the steam and water fronts, respectively:$$A_{\mathrm{s}}=\frac{2h}{\pi }\left(\frac{h}{2}+x\right)$$(17) $$A_{\mathrm{w}}=xh+\frac{\pi h^2}{8}-A_{\mathrm{s}}$$(18)

Shutler and Boberg [11] proposed that piston-like flow through a thin layer is assumed; therefore, the isothermal two-phase flow theory of Buckley and Leverett [29] can be introduced into isothermal zones. Based on Shutler’s theory, the horizontal displacement of the cold water front can be expressed as:$$L_{\mathrm{f}}=\frac{f_{\mathrm{w}}^{\mathrm{\prime }}\left(S_{\mathrm{wf}}\right)}{\varphi h\mathrm{dl}}{\int }_{t_c}^t\frac{v_{m_{\mathrm{c}}}}{{\rho }_{\mathrm{w}}}\mathrm{d}t+x$$(19)where is the slope of the tangent of the water fractional flow curve at the cold water front. Because the solvent penetration depth is extremely thin when the cold water begins to generate, the horizontal displacement of the steam front is assumed to be not affected by the viscosity reduction by the solvent.

Cold water is assumed to block the diffusion of the solvent into the crude oil. Based on assumption (4), the effective diffusion coefficient of the solvent can be given by Figure 5 and Equations (20)-(22):$${\overline{S}}_{\mathrm{w}}=\{\begin{array}{l}{S}_{\mathrm{wc}}+\frac{1}{L_{\mathrm{e}}\varphi h\mathrm{dl}}{\int }_{t_c}^t\frac{v_{m_{\mathrm{c}}}}{{\rho }_{\mathrm{w}}}\mathrm{d}t,\begin{array}{l}\\ \end{array}(L_{\mathrm{f}}\le x_2)\\ S_{\mathrm{wt}}+\frac{1}{L_{\mathrm{e}}\varphi h\mathrm{dl}}\left[1-f_{\mathrm{w}}\left(S_{\mathrm{wt}}\right)\right]{\int }_{t_c}^t\frac{v_{m_{\mathrm{c}}}}{{\rho }_{\mathrm{w}}}\mathrm{d}t,\begin{array}{l}\\ \end{array}(L_{\mathrm{f}}>x_2)\end{array}$$(20) $$L_{\mathrm{e}}=x_2-x_1$$(21) $$D_t=D\frac{1-{\overline{S}}_{\mathrm{w}}}{1-S_{\mathrm{wc}}}$$(22)where L_e is the effective operating thickness of the solvent; is the average water saturation in L_e; x₁ and x₂ are the coordinates of the steam edge and the effective operating front, respectively; and D_t is the effective diffusion coefficient of solvent at time t.

The solvent diffusion rate and oil production rate can also be calculated by Equations (6) and (7).

1.3 Steam Condensate-Producing Stage

The steam condensate-producing stage begins when the steam condensate is produced from the production well. The water saturation profile between the steam edge and production well is determined by the oil viscosity profile, which is in turn influenced by the solvent concentration profile.

The following mixing rule is used for calculating the mixture viscosity of the oil phase [30]:$$\ln {\mu }_{\mathrm{m}}=\left(1-C_{\mathrm{si}}{e}^{-\frac{\epsilon }{{\delta }_c}}\right)\ln {\mu }_{\mathrm{o}}+C_{\mathrm{si}}{e}^{-\frac{\epsilon }{{\delta }_c}}\ln {\mu }_{\mathrm{s}}$$(23)where μ_o, μ_s, and μ_m are the viscosities of the heavy oil, solvent, and mixture of heavy oil and solvent, respectively.

By employing the theory of Buckley and Leverett, the slope of the tangent of the water fractional flow curve at the production well can be written as:$$f_{\mathrm{w}}^{\mathrm{\prime }}\left(S_{\mathrm{wp}}\right)=\frac{\epsilon \varphi h\mathrm{dl}}{{\int }_{t_c}^t\frac{v_{m_{\mathrm{c}}}}{{\rho }_{\mathrm{w}}}\mathrm{dt}}$$(24)

The oil production rate can be rewritten from Equation (7) by considering the water and solvent production:$$q_{\mathrm{oi}}=\left[{\rho }_{\mathrm{o}}\varphi \Delta S_{\mathrm{o}}\mathrm{dl}\frac{\mathrm{d}r}{\mathrm{d}t}+{\rho }_{\mathrm{o}}\varphi \left(1-S_{\mathrm{wc}}\right)\mathrm{dl}{\overline{C}}_{\mathrm{s}}\frac{\mathrm{d}{A}_{{\delta }_c\left(t\right)}}{\mathrm{d}t}\right]f_{\mathrm{w}}\left(S_{\mathrm{wp}}\right)\left(1-C_{\mathrm{s}}\right)$$(25)where S_wp is the water saturation at the production well, which is determined by and μ_m, and is the fractional flow of water at the production well, which can be acquired by considering μ_m and S_wp.

2 Calculation for Mathematical Model

2.1 Model Treatment

A semi-analytical method is used to solve the model due to the complexity of the equations. To solve the semi-analytical model conveniently, several equations must first be solved.

The equation describing the movement of the steam front can be obtained by integrating Equations (11)-(14).$$\begin{array}{l}{H}_{\mathrm{o}}-2\hspace{1em}{\int }_0^t(T_{\mathrm{w}}-T_{\mathrm{r}})\sqrt{\frac{{\lambda }_{\mathrm{cap}}{\rho }_{\mathrm{cap}}{c}_{\mathrm{cap}}}{\pi (t-\tau )}}\frac{\mathrm{d}x}{\mathrm{d}\tau }\mathrm{d}\tau =\frac{2}{3}\left[{\rho }_{\mathrm{r}}{c}_{\mathrm{Pr}}\frac{(1-\varphi )}{\varphi }+{\rho }_{\mathrm{o}}{S}_{\mathrm{or}}{c}_{\mathrm{Po}}+{\rho }_{\mathrm{s}}\Delta S_{\mathrm{o}}{c}_{\mathrm{Ps}}+{\rho }_{\mathrm{w}}{S}_{\mathrm{wc}}{c}_{\mathrm{Pw}}\right]\left(T_{\mathrm{s}}-T_{\mathrm{r}}\right)\varphi \frac{\mathrm{d}{A}_T}{\mathrm{d}t}\\ \begin{array}{llllllllllll}& & & & & & & & & & & \end{array}+ \frac{1}{3}\left[{\rho }_{\mathrm{r}}{c}_{\mathrm{Pr}}\frac{(1-\varphi )}{\varphi }+{\rho }_{\mathrm{o}}{S}_{\mathrm{or}}{c}_{\mathrm{Po}}+ {\rho }_{\mathrm{s}}\Delta S_{\mathrm{o}}{c}_{\mathrm{Pw}}\right]\left(T_{\mathrm{w}}-T_{\mathrm{r}}\right)\varphi \frac{\mathrm{d}{A}_T}{\mathrm{d}t}\end{array}$$(26)

Equation (26) belongs to the Volterra integral equations of the second kind, which makes it difficult to solve directly. In this article, Equation (26) is solved with the support of a Laplace transform and the displacement of the steam front is given as:$$x={\int }_{t_c}^t\frac{H_{\mathrm{o}}}{B}{e}^{\left\{{\left[\frac{A}{B}\Gamma \left(0.5\right)\right]}^2\tau \right\}}\mathrm{erfc}\left(\frac{A}{B}\Gamma \left(0.5\right)\sqrt{\tau }\right)\mathrm{d}\tau$$(27)where$$A=2\left(T_{\mathrm{w}}-T_{\mathrm{r}}\right)\sqrt{\frac{{\lambda }_{\mathrm{cap}}{\rho }_{\mathrm{cap}}{c}_{\mathrm{cap}}}{\pi \left(t-\tau \right)}}$$(28) $$B=\frac{2}{3}\left[{\rho }_{\mathrm{r}}{c}_{\mathrm{Pr}}\frac{(1-\varphi )}{\varphi }+{\rho }_{\mathrm{o}}{S}_{\mathrm{or}}{c}_{\mathrm{Po}}+{\rho }_{\mathrm{s}}\Delta S_{\mathrm{o}}{c}_{\mathrm{Ps}}+{\rho }_{\mathrm{w}}{S}_{\mathrm{wc}}{c}_{\mathrm{Pw}}\right]\left(T_{\mathrm{s}}-T_{\mathrm{r}}\right)\varphi +\frac{1}{3}\left[{\rho }_{\mathrm{r}}{c}_{\mathrm{Pr}}\frac{(1-\varphi )}{\varphi }+{\rho }_{\mathrm{o}}{S}_{\mathrm{or}}{c}_{\mathrm{Po}}+{\rho }_{\mathrm{s}}\Delta S_{\mathrm{o}}{c}_{\mathrm{Pw}}\right]\left(T_{\mathrm{w}}-T_{\mathrm{r}}\right)\varphi$$(29)

2.2 Calculation Flow for Mathematical Model

Because the solvent diffusion rate and steam front velocity change with time, time is first discretized into small intervals. As mentioned above, the following steps are used for the calculation of the oil production rate:

1. for a given heat injection rate and solvent component, Equations (1), (2), (5), and (6) are coupled and solved for the steam front velocity U, radius of the steam zone area r, solvent diffusion rate dδ _c /dt, and solvent penetration depth δ _c ;

2. the production rate in the first stage is calculated using Equations (7)-(10);

3. repeat steps 2 and 3 substituting δ _c calculated in step 2 into Equation (6) until the steam front reaches the understrata and overburden, i.e., when the process enters the next stage;

4. the steam front velocity dx/dt and the steam front displacement x in the second stage are calculated using Equations (15) and (27);

5. the cold steam condensate generating rate is calculated using Equations (16)-(18). Then, Equations (5), (6), and (19)-(22) are coupled and solved for the solvent dδ _c /dt, δ _c , and the cold steam condensate distribution beyond the edge of the steam front;

6. the production rate in the second stage is also calculated using Equations (7)-(10);

7. repeat steps 5-7 substituting δ _c calculated in step 6 into Equation (6) until the cold steam condensate reaches the production well, when the process enters the third stage;

8. , dδ _c /dt, and δ _c are calculated in the same manner as step 6. The water fractional flow f _w(S _wo) in the third stage is calculated using Equations (5), (23), and (24). Then, the oil production rate of the third stage can be calculated using Equation (25);

9. repeat steps 9 substituting δ _c calculated in step 9 into Equation (6) until f _w(S _wo) = 1, when the process ends.

3 Results and Discussion

3.1 Model Validation

It is difficult to address all combinations of different solvents by experiments and numerical simulation, making the validation of the SESF process far from comprehensive. Given this situation, the model with solvent concentration value of zero is validated to ensure that the proposed model is a credible base model for the SESF calculation. This method of validation is considered as reasonable in this article because the major objective of the model is to investigate the effects of the solvent properties on the production of SESF, not that of a specific solvent.

STARS is an effective thermal recovery reservoir numerical simulator developed by the CMG Company. As the results of STARS are exact numerical solutions that can closely reveal the real reservoir production process, it frequently serves as a reference to test and verify the correctness of other models. Hence, to validate the correctness of the proposed model, the results of the model are compared with those of STARS. The parameters used in the proposed model and STARS are all based on the average property parameters of one oil field in China. The well spacing between the production well and injection well is 40 m. The thickness of the reservoir is 7 m. The steam injection rate per unit length is 0.09 t/md, which is low because of the thin thickness of the reservoir. The temperature of the steam is 300 °C and the oil viscosity is 300 mPa s at the reservoir condition. Table 1 lists the reservoir properties and injection parameters used in the calculation of both the new model and STARS.

The steam zone development comparison between the proposed model and STARS is made using Figure 6. The geometries and movements of the steam front and hot water front predicted by the proposed model are similar to the numerical simulation results, which also proves the reasonability of assumption (2). Figure 7 depicts the cold steam condensate front position at the end of the steam front spreading stage. The difference of the condensate front position between the two methods is only 2 m, which also proves the feasibility of the Buckley-Leverett theory in the steam-flooding process in a thin layer.

Figure 8 displays the oil production rate as a function of the steam flooding production time. The oil production rate increases before the steam front reaches the understrata and overburden. In the second stage, the production rate remains virtually stable because the hot water zone between the steam zone and boundaries alleviates the heat loss. After the steam condensate reaches the production well, the oil production rate begins to decrease sharply; in this stage the steam zone and hot water zone stop expanding at the point when the heat loss rate to the overburden and understrata equals the heat injection rate [31, 32].

RMSRE (Rooting-Mean-Square of Relative Error) (Eq. 30) and RE (Relative Error) (Eq. 31) are employed to evaluate the accuracy of the proposed model in calculating the cumulative oil production and oil production rate [33]:$$\mathrm{RMSRE}=\sqrt{\frac{1}{\mathrm{NPTS}}\sum _{i=1}^{\mathrm{NPTS}}{\left[\frac{y_i^{\mathrm{cal}}-y_i^{\mathrm{sim}}}{y_i^{\mathrm{sim}}}\right]}^2}\times 100\%$$(30) $${\mathrm{RE}}_i=\frac{y_i^{\mathrm{cal}}-y_i^{\mathrm{sim}}}{y_i^{\mathrm{sim}}}\times 100\%$$(31)where NPTS is the number of data points; is the calculated value at point i by the proposed model; and is the calculated value at point i by STARS.

In Figure 9, the cumulative oil calculated by the proposed model and STARS are compared. The RE of the cumulative oil for each point is less than 7%; the RMSRE is only 3.5%. Thus, the proposed model is accurate in calculating the steam-flooding process. Figure 10 indicates the RE distribution of the oil production rate with time. It can be clearly observed that the RE is significantly higher in the third stage than in the first and second stages. According to our analyses, the greater difference between the two results in the third stage is largely due to the one-dimensional assumption of the Buckley-Leverett theory. The RE and absolute error (difference between the calculated oil production rate by the proposed model and STARS) in the later portion of the third stage is less than that in the middle of the same stage. As is known from a previous study, the solvent mainly affects the early and late stages of the steam injection process [10, 34], which is consistent with the low RE stage in the proposed model. Therefore, the proposed model is sufficiently accurate to calculate the SESF process.

After the validation of the proposed model, the oil production rate of the steam flooding is compared with that of SESF, as illustrated in Figure 11. The solubility of the solvent at the steam front is 0.2, the diffusion coefficient of the solvent into the heavy oil is 5.5 × 10⁻² m²/d, which is one order of magnitude greater than the Das and Butler’s results that pertain to Heleshaw laboratory models of VAPEX [35]. Gupta and Gittins [13] regard the higher diffusion coefficient as reasonable because of the higher temperature and existence of a blanket layer. The longitudinal dispersity factor is 5 × 10⁻³ m and the viscosity of solvent is 4.5 mPa s. It can be seen that oil rate increases rapidly in the first stage with the contact area increasing until the steam front reaches the boundaries when the oil rate reaches its highest point. Then, it begins to decline owing to the lower solvent concentration gradient beyond the steam front and the blocking effects of the steam condensate on the solvent diffusion. Moreover, the production time is shortened because of the improvement of the oil viscosity profile beyond the edge of the steam caused by the dissolution of the solvent into the heavy oil. These are consistent with the working stage of the solvent mentioned in other literatures [10, 34], which in turn verifies the proposed model.

3.2 Effects of Solvent Properties on SESF

Figure 12 displays the effects of different effective operating thicknesses of the solvent on both the production rate and reduced steam oil ratio. The reduced steam oil ratio is defined as the steam oil ratio difference between SESF and steam flooding, which represents the thermal efficiency increment by the solvent. We can see that the oil rate decreases slower with higher effective operating thickness, while the production time remains unaffected. The larger the effective operating thickness, the greater the thermal efficiency increment. This means that by pushing the solvent further into the reservoir, the thermal efficiency can be increased without additional consumption of solvent, thus increasing the economic benefits. This is similar to the phenomenon in the experiments of Redford [9].

Figure 13 indicates that the solvent solubility in heavy oil also has significant effects on the SESF process. It can be observed that both the oil production rate in the early stage and thermal efficiency increment increase with the solvent solubility. Moreover, according to Equations (23) and (24), for the benefit of a lower mixture viscosity beyond the steam edge, a higher solvent solubility tends to result in a shorter production time without decreasing the cumulative oil production, such that heavy oil is recovered considerably more efficiently.

However, it should be noted that a larger operating thickness always means a lower solvent solubility in terms of a single solvent. Consequently, an optimal solvent must be carefully selected for the best performance of SESF when injecting a single solvent with steam. Multi-solvent injection, which is rarely mentioned in other analytical models, is provided with a simple method to simulate in the proposed model using different values of operating thickness and solubility.

3.3 Effects of Injection Strategy on SESF

Figure 14 illustrates the effects of different injection rates on SESF production. It can be concluded that the differences between the largest oil rate and the stable oil rate are virtually the same among different injection rates. The reason is that the diffusion rate is largely determined by the concentration gradient and contact area, not the interface velocity. Further, the generating of water within the operating thickness accelerates with an increasing injection rate, which results in a sharp decline of oil production rate at an early recovery stage for a high injection rate. Moreover, Figure 14 indicates that the thermal efficiency increment increases with a decrease of the injection rate because of the longer contacting time of the solvent with the oil and thus relatively shorter production times compared with steam flooding.

Figure 15 displays the effects of different well spacing on both the production rate and reduced steam oil ratio. It is clear that the oil rates in the early stage are identical for different well spacing. Moreover, the larger the well spacing, the lower the thermal efficiency increment. This is because the diffusion rates are the same with different well spacing and thus the oil is unaffected by the solvent increases related to the well spacing. Therefore, solvents with higher operating thickness, which means wider working range, are more suitable for longer well spacing.

From Figures 14 and 15, it can also be concluded that despite the fact that a higher injection rate may result in lower thermal efficiency increment when well spacing is constant, it may be more suitable for longer well spacing owing to the improvement of saturation and oil viscosity profiles beyond the steam edge caused by the longer contact time. Thus, an optimal injection rate must be carefully selected for the best performance of SESF with a given well spacing.

Conclusion

In this study, a semi-analytical model for predicting SESF production performance was established. Then the theoretical model was validated using STARS results and error analysis. The production performances of steam flooding and SESF were also compared. Moreover, with the theoretical model, the effects of solvent properties, injection strategies on oil production rate, and thermal efficiency increment were investigated. The findings can be summarized as follows:

1. compared with the STARS results for one case based on average reservoir parameters, the proposed model proved to be sufficiently accurate in calculating the SESF production performance;

2. SESF can achieve a considerably higher oil production rate at the early recovery stage than steam flooding. Further, the production time is shortened owing to the improvement of the viscosity profile beyond the edge of the steam, which is caused by the solvent dissolution into the heavy oil;

3. the operating thickness and solubility of a single solvent are proportional to the oil production increment of SESF; however, they cannot be increased simultaneously. Thus, an optimal solvent must be carefully selected for the best performance of SESF when injecting a single solvent with steam. The proposed model also provided a simple method to simulate multi-solvent injection using different values of operating thickness and solubility;

4. a higher injection rate can result in a lower thermal efficiency increment when well spacing is constant. However, it may be suitable for longer well spacing owing to the improvement of the viscosity profile beyond the steam edge, which is caused by longer contact time between the solvent and crude oil.

Acknowledgments

The authors wish to thank the Research Institute of Petroleum Exploration & Development, Liaohe Oilfield Company, PetroChina (2013-JS-9399). This work was also supported in part by a grant from the National Science and Technology Major Projects of China (2016ZX05012-005-004) for financial support. Discussions with colleagues of an oil reservoir numerical simulation group provided great insight for our paper.

References

All Tables

All Figures