© I. Malinouskaya et al., published by IFP Energies nouvelles, 2018

1 Introduction

Over the past decade, the method of Steam Assisted Gravity Drainage (SAGD) has become increasingly important for heavy oil recovery because of the large petroleum reserves accessible using this thermal Enhanced Oil recovery (EOR) process . In addition to the heat conduction from steam to oil leading to the decrease of oil viscosity and density, this process is also associated with geomechanical effects [1]  . Increasing the pressure and temperature leads to a modification of petrophysical properties, as porosity and permeability [2]. Moreover, the modification also induces changes in the productivity index of the wells, mostly forgotten in the reservoir simulation.

There are several models using either iterative either full coupling of the reservoir simulations with geomechanics [3–6]. Most of these models are developed for the fractured media and thermal effects are neglected. An iterative coupling proposed by [6] takes into account the porosity changes according to the pressure, temperature and total stress variations. However, this model was applied only on porothermoelastic materials. Then, [7] proposed different methods for the stress and permeability approximations according to the pressure variations. They demonstrated that the approximative methods yield the results very close to the direct coupled simulations. The important temperature changes during SAGD process suggests to take them into account for the porosity and permeability estimations [1]. Thus, the geomechanical model proposed in this work is based on the empirical relationships between the petrophysical properties and both, the pressure and temperature, variations.

The overall objective of this work is to study the effect of petrophysical properties modification which occurs during SAGD process using a numerical model based on a meandering fluvial reservoir and particularly focus on the geomechanical productivity index. To focus on this type of reservoirs, we built a 3D reservoir model based on 3D outcrop observations of the Scalby Formation (Middle Jurassic Ravenscar Group) located in the Yorkshire, UK. This reservoir model was studied by [8] to estimate the impact of upscaling but without considering the effect of pressure and temperature changes on petrophysical properties. A 1D model derived from this model was also used by [9] to study the H₂S production during SAGD process, but always without considering the geomechanical effect.

This paper is organized as follows. First, in Section 2, we describe the geological context of the reservoir which is studied here. Then, the details of our geomechanical model and the modification of productivity index are given in Section 3. The conditions of the SAGD process simulation and the results are discussed in Section 4. Finally, in Section 4, we provide some perspectives and conclusions on this work.

2 Reservoir model

Before stating our geomechanical model, we first need to clarify the context by setting the reservoir model. The geological model used for the simulations of the SAGD recovery process is obtained from 3D outcrop observations of the Scalby Formation analogue. It exhibits heterogeneities specific to fluvial and estuarine environments and, therefore, the corresponding reservoir model is a highly complex system. Thus, in order to describe accurately the heterogeneities distribution within the reservoir, the geological model was constructed on the fine grid using a geostatistical approach.

The lithofacies are defined according to the depositional history as described by [8]. In terms of lithofacies distribution and petrophysical properties, the Scalby Formation is very similar to the McMurray Formation of Hangingstone heavy oil field in Athabasca (Alberta, Canada). Therefore, the well-log data from Hangingstone field is used to define the lithofacies of our model. Five lithofacies are distinguished [10]:

• Lithofacies 1 is a clean medium- to coarse-grained sandstone which often corresponds to braided stacked channel deposits;

• Lithofacies 2 is a medium-grained sandstone associated with channel infill fluvial and estuarine sandstones;

• Lithofacies 3 is fine-grained sandstones which occur as fining-upward intervals at the top of channel fills and overbank deposits associated with sandy Inclined Heterolithic Stratification (IHS) in the estuary setting and with sandy tidal flats;

• Lithofacies 4 corresponds to silty shales facies which is associated with the main heterolithic facies associations of the reservoir, represented by tidally influenced point bar facies, estuarine IHS, mud flat and fine shaly overbank deposits;

• Lithofacies 5 is shaly facies related to channel abandonment mud plugs, floodplain or coastal plain shales, muddy IHS or distal bay deposits.

The distribution of the facies is shown in Figure 1 for a half field with a cut along the wells (at Y = 75 m).

In order to accurately account for the heterogeneities distribution within the reservoir, a very fine geological grid is built using geostatistical modeling. Then, the grid is populated with the lithofacies and corresponding initial petrophysical properties such as the porosity, φ, and horizontal and vertical permeability, K_h and K_v, respectively. These properties are summarized in Table 1.

For the sake of simplicity, all facies have the same water-oil and gas-oil relative permeability curves which are shown in Figure 2.

For the reservoir simulation, the initial oil saturation is set to 0.85. The irreducible water saturation, the residual oil saturation to waterflood and to steamflood are set to 0.15, to 0.20 and 0.10, respectively. The oil properties are given in Table 2 [8]. The top of the reservoir is at 250 m of depth. The initial reservoir temperature is uniform, T₀ = 10  ° C, and the reference pressure P₀ = 20 bar is defined at 220 m of depth.

The injector and producer wells are oriented along X-axis and situated in the middle of the reservoir according to Y-axis and at 9.2 m and 15.2-m depth from the top of the reservoir for injector and producer, respectively.

Fig. 1 Lithofacies distribution: the model is cut allong the wells (Y = 75 m).

Fig. 2 Water-oil (a) and gas-oil (b) relative permeabilities.

3 Geomechanical model

We are now in a position to work out the simplified geomechanical model. The basic idea is to use analytical and empirical relationships for the porosity, permeability and the productivity index to take into account the pressure and temperature dependencies of the geomechanical parameters in the reservoir simulations.

The mechanical properties such as Young modulus E and Poisson coefficient ν, for each lithofacies are deduced from the measurements provided by [11] and given in Table 3. For the geomechanical model, we also need the coefficient of the rock compressibility c_p and the rock thermal expansion c_T. Using the hypothesis of the constant stress, c_p is found from $${\phi }_0{c}_p=\frac{3(1-2\nu )}{E}\text{,}$$(1) where φ₀ is the initial porosity.

The thermal expansion coefficient c_T is deduced from [11] for each kind of facies. For the sandstones (lithofacies 1–3), one can use a constant value of c_T corresponding to the sands. For the shales (lithofacies 4–5), it is a temperature dependent function obtained by square fit of the measurements provided by [11]. In Table 3, we can remark that for lithofacies 4–5, c_T becomes negative when the temperature increases. Indeed, when the thermal expansion coefficient is positive, the effective stress increases with temperature whereas it decreases as pore pressure rises. When c_T is negative the effects of heating and steam injection cumulate [12]. These parameters are summarized in Table 3.

High variations in pressure and temperature during SAGD process have a strong impact on the strain state of the reservoir. Therefore, P and T modify the petrophysical properties such as porosity and permeability and it should be taken into account in the well productivity estimation. A direct coupling of the reservoir and geomechanical simulations [13,14] can reveal expensive in terms of computational time, thus, a simplified geomechanical model has been developed. It allows us to take into account the geomechanical effects on the porosity, permeability and, consequently, on the well productivity index via either analytical either empirical relationships. It should be noted that this approach is based on the simplifying assumption of constant total stress and that in this context sand dilation is neglected.

The variations of the porosity are determined by a first derivative of the corresponding stress-strain relationship equations under hypothesis of the small deformations which yields $$\phi ={\phi }_0\exp \left[c_p(P-P_0)+c_T(T-T_0)\right]\text{,}$$(2) where P₀ and T₀ are the initial reservoir pressure and temperature, respectively.

Then, using a multiplier M_i(P, T) for the three directions i = X, Y, Z, the new permeability $K_i^{\text{eff}}$ is obtained by $$K_i^{\text{eff}}\left(P,T\right)=M_i(P,T)K_i\text{,}$$(3) where K_i is the initial permeability.

The multiplier is defined using an appropriate porosity-permeability relationship. For the lithofacies 1–3, the Touhidi-Baghini expression [2] provides M_i(P, T) in the empirical form for i = X, Y, Z as $$M_i(P,T)=\frac{K_i^{\text{eff}}}{K_i}=\exp \left(\frac{c_i{ϵ}_v}{{\phi }_0}\right)\text{,}$$(4) with c_X = c_Y = 2 and c_Z = 5 are the material parameters and ϵ_v the volumetric strain that is computed taking into account the pressure and temperature variations. $${ϵ}_v={\phi }_0{c}_p(P-P_0)+c_T(T-T_0)\mathrm{}\text{.}$$(5)Note that the Touhidi-Baghini expression [2] is obtained for McMurray formation oil sands for the range of porosity φ ≥ 0.3. Therefore, for the shales (lithofacies 4 and 5), we consider that the approximation of Kozeny-Carman [15] is more appropriate. Thus, the multipliers are calculated as $$M_i(P,T)=\frac{K_i^{\text{eff}}}{K_i}=\frac{{\phi }^3{\left(1-{\phi }_0\right)}^2}{{\phi }_0^3{\left(1-\phi \right)}^2}\text{.}$$(6)

Then, the permeability modifications are taken into account for the well injectivity and productivity index calculations by using Peaceman formulation [3,16,17] $$PI=\frac{2\pi \sqrt{K_X^{\text{eff}}{K}_Y^{\text{eff}}}h}{\ln \left(\frac{r_0}{r_w}\right)+s}\text{,}$$(7) where h is an effective thickness of the well, s the skin factor and r_w the well radius.The drainage radius r₀ is obtained by $$r_0=0.28\frac{{\left[{\left(\frac{K_X^{\text{eff}}}{K_Y^{\text{eff}}}\right)}^{1/2}\mathrm{\Delta }{x}^2+{\left(\frac{K_Y^{\text{eff}}}{K_X^{\text{eff}}}\right)}^{1/2}\mathrm{\Delta }{y}^2\right]}^{1/2}}{{\left[\frac{K_X^{\text{eff}}}{K_Y^{\text{eff}}}\right]}^{1/4}+{\left[\frac{K_Y^{\text{eff}}}{K_X^{\text{eff}}}\right]}^{1/4}}\text{.}$$(8)Using (3), the well productivity indices (7,8) are updated according to the pressure and temperature variations as $$PI(P,T)=\frac{2\pi {\left[M_X(P,T)M_Y(P,T)K_X{K}_Y\right]}^{1/2}h}{\ln \left(\frac{r_0}{r_w}\right)+s}\text{,}$$(9)with $$r_0=0.28\frac{{\left[M_Y(P,T)K_Y\mathrm{\Delta }{x}^2+M_X(P,T)K_X\mathrm{\Delta }{y}^2\right]}^{1/2}}{{\left[M_Y{K}_Y\right]}^{1/2}+{\left[M_X{K}_X\right]}^{1/2}}\text{.}$$(10)

4 SAGD simulation

The SAGD process was simulated using the previous geomechanical model in the reservoir simulator used by [8], its details can be found using the link http://www.beicip.com/reservoir-simulation

First, the wells are heated during 120 days at 260 °C. Then, the steam of quality 0.8 is injected with the maximal well bottom hole pressure of 50 bar. The minimal bottom hole pressure of the production well is set to 5 bar. The production rate is controlled by maintaining the difference between production and injection well temperature between 20 and 30 °C. The oil production is simulated during 2000 days.

We are going to compare three simulations, one accounts the effects of pressure and temperature variations on the porosity, permeability and well productivity index, the second one contains only geomechanical effects on the porosity and permeability but not on PI, and the third one uses constant porosity and permeability and PI (no geomechanical effects).

In Figure 3a, b, the temperature is more diffused when geomechanical effects are taken into account. Therefore, due to the pressure and temperature variations on the porosity, permeability and productivity indices, we obtain a better steam propagation within the steam chamber (T ≥ 100  ° C) of the reservoir model. The same effect has been observed by [1]. It was demonstrated that the porosity and permeability increase due to heating accelerates the growth of the steam chamber. Comparing Figure 3c–d, e–f and h–g, the medium is more saturated with gas and water for the case using our geomechanical model and, consequently, more oil has been produced. This is also confirmed by the oil production curve shown for the producer well in Figure 5a. Due to the fracturing of the shales, the steam chamber is extended and thus, more oil can be produced. At the end of the simulation period, the difference in oil production at surface conditions between two cases is about 10%.

Figure 4a shows the effects of the pressure and the temperature on PI. During wells heating, the productivity index significantly increases. Then, at the beginning of the production, the pressure increases and applies small fluctuations of the PI. Once, the pressure and temperature in the well are stabilized, the productivity index becomes constant. However, the resulting value of PI is significantly higher in comparison to the case when pressure and temperature variations are not taken into account. We remark an increasing of 20% or 32% for the examples in Figure 4a.

The heterogeneities distribution also plays an important role on the well productivity index behavior as already showed by [18] which analyse the impact of shale distribution and volume fraction on oil production, and by [19] where it was proposed a set of input attributes for correlating the reservoir heterogeneities to SAGD production performance. For example, at layer X = 0 m (Fig. 4b) the lithofacies are mostly sandstones compared to the layer X = 180 m (Fig. 4c) which is mostly composed of shales on the top of the wells. Therefore, taking into account the formulations (4) and (6), the pressure and temperature effects on PI are more significant for shales than for sands (Fig. 4a). The negative thermal expansion coefficients used for the lithofacies 4–5 in our geomechanical model at high temperature, lead to the fracturing of the shales. Therefore, the porosity, permeability and the productivity index increase within these regions.

In order to underline the role of the PI modifications due to the temperature and pressure variations, the results are also compared with a simulation without geomechanical effects on the productivity index. If we don’t take into account the PI variations, the steam injection flow rate will not be adjusted to the modified porosity and permeability and it will take longer time to fill the pore volume which became larger due to the temperature and pressure effects. Therefore, the oil flow rate at the beginning of the steam injection is even smaller than for the model without any geomechanical effect (Fig. 5b). However, once the available volume is fitted, the oil recovery increases and after simulation of 2000 days it becomes 3% smaller than in the full geomechanical model but about 8% higher than for the case without any geomechanical effect.

Fig. 3 Field properties in a cut at Y = 75 m: T is a temperature, Sg is a gas saturation, Sw is a water saturation and So is an oil saturation; upper black line corresponds to the injector well and bottom one to the producer well. Data corresponds to a steam chamber with minimal temperature of 100 °C.

Fig. 5 Cumulative oil production (a) and oil flow rate (b): (red) − with geomechanical effects, (blue) − without geomechanical effects, (green) − without geomechanical effects for the PI calculations only.

Fig. 4 Well productivity index at X = 0 m (blue) and X = 180 m (black) compared to the constant PI for the case without geomechanical effects (red). (b) and (c) are the field maps corresponding to X = 0 m and X = 180 m, respectively; black points on (b–c) are the producer and injector wells position; the facies colors are the same as in Figure 1.

5 Conclusion

The simplified geomechanical model has been applied to the SAGD process simulation. It allows us to take into account the rock heterogeneities and high pressure and temperature variations for the estimation of the porosity and permeability. Thus, the well productivity can be calculated more accurately according to these modifications of the petrophysical properties.

It was demonstrated that for SAGD process, the simulation without geomechanical effects yields underestimation of the production. The stress induced into the media by heating can change the porosity and permeability. Therefore, it is essential to use an appropriate relationship between these parameters. As a first approach we used the semi-analytical formulations [2] for sands and [15] for the shales. This leads to the increase of the productivity estimation by about 10%.

The results presented here seem to be qualitatively relevant. Nevertheless, performing real coupling of the reservoir simulator with the geomechanical one, we could observe the effects of the plastic strain and rock dilatancy which could appear in such kind of models [1,14,20]. A future work can be to compare our simplified model with geomechanically coupled fluid flow model.

References

All Tables

All Figures

	Fig. 1 Lithofacies distribution: the model is cut allong the wells (Y = 75 m).
In the text

	Fig. 2 Water-oil (a) and gas-oil (b) relative permeabilities.
In the text

	Fig. 3 Field properties in a cut at Y = 75 m: T is a temperature, Sg is a gas saturation, Sw is a water saturation and So is an oil saturation; upper black line corresponds to the injector well and bottom one to the producer well. Data corresponds to a steam chamber with minimal temperature of 100 °C.
In the text

	Fig. 5 Cumulative oil production (a) and oil flow rate (b): (red) − with geomechanical effects, (blue) − without geomechanical effects, (green) − without geomechanical effects for the PI calculations only.
In the text

	Fig. 4 Well productivity index at X = 0 m (blue) and X = 180 m (black) compared to the constant PI for the case without geomechanical effects (red). (b) and (c) are the field maps corresponding to X = 0 m and X = 180 m, respectively; black points on (b–c) are the producer and injector wells position; the facies colors are the same as in Figure 1.
In the text