Numerically Coupled Thermo-Hydro-Mechanical Analyses of Ultra-Heavy Oil Reservoirs during the Micro-Fracturing Stage

Abstract

Water injection in SAGD well pairs has rapidly developed to improve the circulation preheating efficiency and the steam chamber propagation speed. However, a key problem is to evaluate the coupled thermo-geomechanical effects of water injection under complex operations. In this paper, a coupled thermo-hydro-mechanical model considering skeleton shear dilation and phase change of bitumen was established. Major conclusions were drawn that there were large temperature changes only in two separate areas around wells. The pore pressures in the reservoir around wells and the base rock right below the production well increased significantly. The void ratio was improved in the reservoir around wells, especially in the inter-well region. The reservoir was lifted, and the largest uplift was right above the wells. This information can guide engineers in properly evaluating the field operations.

1. Introduction

The steam-assisted gravity drainage (SAGD) technique is one of the efficient recovery methods for the development of heavy oils. However, these burning issues usually occur, such as a long-time circulation preheating period, a low speed for the steam chamber propagation, and a low recovery efficiency for the horizontal section. In this regard, the micro-fracturing technique by water injection in SAGD well pairs has rapidly developed to try to overcome the above-mentioned problems. The key problem is to evaluate the coupled thermo-geomechanical effects of water injection under complex field operations.

The related theories of micro-fracturing in heavy oils have been proven in the oilfield practice [1,2,3]. In the process of water injection into oil sands (the ultra-heavy oils), the transient inhomogeneous diffusion of reservoir pore pressure and temperature leads to the dramatic evolutions of mean effective stress, shear stress, and temperature applied to the oil sand structure, generating tensile micro-fracture, shear fracture, and thermal expansion in the wellbore region. The propagation of fractures significantly increases the porosity and permeability in the wellbore area, which is beneficial to the heat transfer and water flow, aiming at a better circulation preheating effect in a shorter time and a lower economic cost [4,5,6]. At present, the studies on the coupled thermo-hydro-mechanical responses in oil sand reservoirs are mostly focused on the core-scale high-temperature high-pressure triaxial experiments (e.g., the Canadian oil sands and Venezuelan ultra-heavy oils) and engineering-scale numerically coupled thermo-geomechanical investigations [7,8,9,10,11,12]. However, there are a few studies on the laboratory experiments and numerical simulation for the Karamay oil sands in Xinjiang, China [13,14,15], which seriously hinder the understanding of reservoir geomechanics in the process of SAGD production.

Currently, SAGD technology is extensively used in oil sand exploitation in Venezuela, Canada, and China, thus, the simulation must coincide with site operations. Nevertheless, traditional dual-horizontal wells have problems such as low oil–gas ratio and long preheating time. Therefore, it is important to evaluate the establishment of vertically hydraulic and thermal communications through geomechanical approaches. The objective of this paper was to simulate a homogeneous oil sand reservoir under hot water injection and to give the coupled thermo-geomechanical responses.

Micro-fracturing by water injection improves the reservoir permeability and fluid flow capability, which is beneficial to the subsequent preheating and production. We proposed a comprehensive numerical model for conventional dual horizontal wells to calculate the coupled reaction in the process of thermal fluid-structural interaction by injecting hot water within a representative oil sand pay zone, considering structure shear dilatancy, bitumen’s phase transition, the permeability development related to elastoplastic strain, and heat and mass transfer theories.

By comparison with other simulations, this investigation could describe the multiple behaviors containing deformation, temperature change, seepage, and bitumen’s phase transition during micro-fracturing. The numerical computation could emulate and foretell the evolution incorporating heat distribution, rock deformation, and formation pressure for different operations. These findings could be employed to accurately predict the coupled thermo-geomechanical responses and to appropriately assess the vertical rise of caprock and pay zone, oil yield variations, and steam usage efficiency for the site engineers.

2. Coupled Thermo-Hydro-Mechanical Processes

2.1. Physical Model

The bitumen in the pore space of oil sands can be viewed as one of the typical amorphous materials; bitumen has three-phase states, namely, liquid, quasi-solid, and glassy [16,17,18]. The phase state of bitumen is related to temperature, and the phase change occurs during a temperature range. It was said in the literature that the fusion of bitumen takes place in the temperature range of about 35~90 °C [19]. When the temperature exceeds the fusion point, the properties of heavy oil and light oil are similar. However, when the temperature becomes lower than the fusion point, the heavy oil is transformed from a liquid into a quasi-solid with the dramatically increasing viscosity [16,17,18]. The viscosity of bitumen at a low temperature is extremely high, and its mechanical behavior is similar to solid [20]. It is confirmed that there is an obvious contrast in reservoir deformation, fluid seepage, and temperature propagation behaviors between the drained zone (or liquid zone) and the undrained zone (or solid zone) [21,22,23]. For instance, the elastic modulus of rocks in the drained zone is about double that in the undrained zone. The compressibility coefficient, thermal conductivity, and thermal diffusion coefficient are about half of those in the undrained zone [14,24].

As shown in Figure 1, due to the phase change of bitumen, the wellbore region during the water injection process was separated into 3 different zones, viz. the drained, partially drained, and undrained zones [23]. In the undrained zone, the solid bitumen played the role of structure and cement, supporting a part of shear stress. In the pore surrounded by solid bitumen, water was the only seepage phase (gas was not considered), but the heat convection was weak because of its low permeability to water. In the drained zone, with the hot liquid being squeezed into the reservoir, there were both heat convection and heat conduction. The solid bitumen was transformed into a liquid, and the viscosity of bitumen decreased dramatically. The mobile bitumen and water form the oil–water mixed flow, and meanwhile, the cementation among mineral particles did not work. The adequate contact among mineral particles made the Karamay oil sands form an interlocking structure. Therefore, the shear dilation potential was improved significantly, and the shear-induced fracture and tensile parting-induced fracture were generated.

To accurately describe the coupled thermo-hydro-mechanical process for the actual pay zone, the three-dimensional geological geometric model was proposed, considering the shape, size, and buried depth of the actual reservoir, caprock, base rock, boreholes, and downhole string in SAGD wells.

Figure 2 is the formation profile perpendicular to the wellbore extension. The possible phase state of bitumen was displayed under the interference of two wells. During the micro-fracturing stage, inter-well thermal communication (which always happens in the SAGD circulation preheating stage) was not easy to achieve because of the short construction time and low temperature of injection water. Therefore, there were two separate drained zones around the horizontal wells, and the far reservoir was undrained. The area between the undrained and the drained zones is called the partially drained zone, whose outer boundary and inner boundary are called bitumen fusion front and rear interfaces, respectively [19]. The shapes of the partially drained zone depend on the reservoir’s physical properties and construction parameters. The thickness of the Karamay oil sand pay zone is commonly thin, and the impermeable caprock and base rock mudstones maintained a perfect closure for the reservoir in vertical directions. To simplify the three-dimensional geometric model, a regular geometric model was established according to the actual reservoir and closure conditions, as shown in Figure 3.

To sum up, the engineering problem of water injection can be transformed into a three-dimensional coupled thermo-hydro-mechanical problem with two moving phase change boundaries. As the thermal, fluid flow, and mechanical properties in three geo-mechanical zones differed significantly, the heterogeneity for these properties was considered using varying parameters.

2.2. Mathematical Model

2.2.1. Heat Transfer

For the oil sands reservoir incorporating mudstone stringers, cover, and bottom rocks, the temperature propagation could be described as [25]:

$$\frac{\partial }{\partial t}\left(\rho C\Delta T\right)+\left(1-\varphi \right){\beta }_{\mathrm{frT}}\Delta T\frac{\partial {\epsilon }_{\mathrm{v}}}{\partial t}+q_{\mathrm{k}}=-{\partial }_{,j}\left(-{\lambda }_{ij}{T}_{,j}-{\rho }_{\mathrm{L}}{C}_{\mathrm{L}}\Delta T\frac{k_{ij}}{{\mu }_{\mathrm{L}}}{p}_{,j}\right)$$

(1)

where t is the time; ρC is the equivalent specific heat capacity for the skeleton and pore fluid; ΔT refers to the variation of temperature; $\varphi$ means the porosity; β_frT is the thermal expansivity for the skeleton; ε_v refers to the bulk strain; q_k means the equivalent heat source intensity for the skeleton and fluid; λ_ij signifies the equivalent heat conductivity for the skeleton and fluid; T is current temperature; ρ_L refers to pore fluids’ density; C_L means the specific heats for liquids; k_ij refers to fluids’ permeability coefficient; μ_L means the dynamic viscosity for fluids; p signifies formation pressure.

2.2.2. Fluid Flow

The pore pressure propagation within the oil sand pay zone could be shown as [25]:

$${\partial }_{,i}\left(\frac{{\rho }_{\mathrm{L}}{k}_{ij}}{{\mu }_{\mathrm{L}}}{p}_{,j}\right)=-\varphi {\alpha }_{\mathrm{Tf}}\frac{\partial T}{\partial t}+\varphi {\alpha }_{\mathrm{p}}\frac{\partial p}{\partial t}+\frac{1-\varphi }{{\rho }_{\mathrm{fr}}}\frac{\mathrm{d}{\rho }_{\mathrm{fr}}}{\mathrm{d}t}+\left(1-\varphi \right)\frac{\partial {\epsilon }_{\mathrm{v}}}{\partial t}$$

(2)

where α_Tf refers to the bulk thermal expansivity; α_p means liquids’ compressibility; ρ_fr is the densities of rock skeleton. The gravity term in the equilibrium equation and mass balance equations are neglected.

2.2.3. Deformation

The water injection-induced deformation could be analyzed from three portions: elastic, plastic, and thermal deformations. Supposing that a variation in the effective stress is (dσ′), the change in the sum of three deformations (dε) becomes:

$$\mathrm{d}\mathit{\epsilon }=\mathrm{d}{\mathit{\epsilon }}^{\mathrm{e}}+\mathrm{d}{\mathit{\epsilon }}^{\mathrm{T}}+\mathrm{d}{\mathit{\epsilon }}^{\mathrm{p}}$$

(3)

where dε^e, dε^p, and dε^T refer to the elastic, plastic, and thermal deformations, severally.

The elastic part could be calculated by the Hook’s law, and the plastic part could be determined by the plastic stress-strain law. The thermo-elastic behaviors could be described by the generalized Hook’s law [25]:

$${{\sigma }^{\prime }}_{ij}=\left(K_{\mathrm{fr}}-\frac{2K_{\mathrm{fr}}{E}_{\mathrm{fr}}}{9K_{\mathrm{fr}}-E_{\mathrm{fr}}}\right){\delta }_{ij}{\mathit{\epsilon }}_{kk}+\frac{6K_{\mathrm{fr}}{E}_{\mathrm{fr}}}{9K_{\mathrm{fr}}-E_{\mathrm{fr}}}{\epsilon }_{ij}-3{\alpha }_{\mathrm{T}}{K}_{\mathrm{fr}}\Delta T{\delta }_{ij}$$

(4)

where ${\sigma }_{ij}$ means the total stress; ${{\sigma }^{\prime }}_{ij}$ refers to the effective stress; α_b refers to the effective stress coefficient; δ_ij means the Kronecker symbol; E_fr refers to the skeleton’s elastic modulus; K_fr and K_s mean the skeleton’s and matrix’s bulk moduli, severally; α_T refers to the linear coefficient of thermal expansion; ε_ij means the strain.

The change in the plastic strain could be described using the flow law:

$$\mathrm{d}{\mathit{\epsilon }}^{\mathrm{p}}=\mathrm{d}{\lambda }_{\mathrm{p}}\frac{\mathrm{d}G}{\mathrm{d}{\mathit{\sigma }}^{\prime }}$$

(5)

where G means the function of plastic potential; dλ_p refers to the plastic proportional multiplier.

In the sections of the shear yield as well as transition, the non-associated flow law can be employed. The function of plastic potential G_t shows:

$$G_{\mathrm{t}}=\sqrt{{\left[\left(p^{\prime }-p_a\right)\tan \beta \left(T\right)\right]}^2+{\left[\frac{q}{1+\alpha -\alpha /\cos \beta \left(T\right)}\right]}^2}$$

(6)

The associated flow law can be adopted to exhibit the function of plastic potential (G_c) in the cap portion:

$$G_{\mathrm{c}}=\sqrt{{\left(p^{\prime }-p_a\right)}^2+{\left[\frac{Rq}{1+\alpha -\alpha /\cos \beta \left(T\right)}\right]}^2}$$

(7)

where p’ is the mean effective stress; β(T) refers to the temperature-related D–P internal friction angles; p_a means the evolved parameters that can control the strain softening or hardening; q is the deviatoric stress; α is the shape factor of the transition zone; R is the eccentricity of the D–P plastic cap.

2.3. Boundary Conditions

2.3.1. Displacement and Stress

According to the in situ stress test data, the initial in situ stress condition is:

$${\left.{\sigma }_{ij}\right|}_{t=0,{\Omega }_{\mathrm{fr}}}={\sigma }_{ij}^0$$

(8)

The boundary condition for stress can be expressed as:

$${\left.X_i\right|}_{t,{\Gamma }_{\mathrm{fr}}}={\sigma }_{ij}\left(t\right)n_j$$

(9)

The boundary condition for displacement is:

$${\left.u_i\right|}_{t,{\Gamma }_{\mathrm{fr}}}={\overline{u}}_i\left(t\right)$$

(10)

where Ω_fr is the space occupied by the rock skeleton; Γ_fr is the boundary of Ω_fr.

2.3.2. Pore Pressure and Flow

According to the formation pressure test data, the initial pore pressure condition is:

$${\left.p\right|}_{t=0,{\Omega }_{\mathrm{l}}}=p_0$$

(11)

The constant flow boundary condition is used:

$${\left.\frac{k}{\mu }\nabla p\cdot n\right|}_{t,{\Gamma }_{\mathrm{l}}}=-\overline{q}\left(t\right)$$

(12)

where Ω_l is the space occupied by pore fluids; Γ_l is the boundary of Ω_l.

2.3.3. Temperature

According to the formation temperature test data, the initial temperature is:

$${\left.T\right|}_{t=0,\Omega }=T_0$$

(13)

The constant temperature boundary condition is used:

$${\left.T\right|}_{t,\Gamma }=T\left(t\right)$$

(14)

where Ω is the space occupied by both rock skeleton and pore fluids; Γ is the boundary of Ω. The specific parameters used for the initial and boundary conditions are shown in the subsequent context.

3. Micro-Fracturing Numerical Simulation

3.1. Model Parameters

The reservoir and operation data of the Well A in the Z1 block in Fengcheng Oilfield was used to give the model geometrical, elastoplastic, fluid flow, heat transfer, and phase change parameters needed for the finite element calculations. The materials used in this study were isotropic in the deformation, fluid flow, and heat transfer models.

3.1.1. Geometrical Parameters

As shown in Figure 3, SAGD well pairs were drilled in Well A in Block Z1, where the upper well was a steam injection well (I well), and the lower well was a production well (P well). According to the previous geological data, the target reservoir owned the thickness of 15 m, the horizontal well length of 380 m, a depth of 372 m for the I well, a vertical distance of 5 m for the SAGD wells, and a vertical distance of 1 m between the P well and the base rock. The caprock and the base rock were both impermeable mudstones. It was assumed that the horizontal minimum principal stress was along with the wellbore extension. The geometric model with a horizontal section of 10 m was used for less calculation.

Therefore, the model geometry is described by the following data. The model sizes in the x, y, and z directions are 30 m, 20 m, and 10 m, respectively. The thickness of the pay zone is 15 m, and the thickness of base rock and caprock is 2.5 m. At first, the well is filled with oil sands. Next, the oil sands in the well are removed, and the screen is installed. The oil sand in direct contact with the screen pipe is assumed. The wellbore diameters of the SAGD production well and steam injection well are both 0.125 m. The inner diameter of the screen is 0.159 m.

3.1.2. Elastoplastic Parameters

According to the high-temperature high-pressure triaxial compression experiments [14], a set of parameters were given for the D–P (Drucker–Prager) elastoplastic constitutive model with a cap for the oil sand under room temperature to 100 °C. The GCTS (Geotechnical Consulting & Testing Systems) RTR 1500 high-temperature high-pressure rock triaxial apparatus was used to obtain the mechanical parameters at different temperatures. As calculated in our published paper [23], when the temperature exceeded 94 °C, the reservoir belonged to the drained zone. These critical elastoplastic parameters used are listed in

Table 1. Elastoplastic parameters of Karamay oil sands at 20~100 °C [14].

Parameters	Undrained Zone			Drained Zone
	20 °C	45 °C	70 °C	100 °C
Young’s modulus (MPa)	663	421	388	724
Bulk modulus (MPa)	806	1000	1217	1125
Poisson’s ratio (-)	0.36	0.43	0.45	0.39
Biot’s coefficient (-)	0.96	0.97	0.99	0.99
D–P friction angle (°)	44.3	29.2	52.5	22.7
D–P cohesion (MPa)	2.26	3.31	0.8	0
Cap parameter pa (MPa)	2.03	1.98	2.2	2.93

.

As shown in Table 1, when the pores were filled with the solid bitumen, with the increase in temperature, the elastic modulus of Karamay oil sands gradually dropped, while the Poisson’s ratio, Biot’s coefficient, and bulk modulus rose. However, when the mobile bitumen was squeezed out of the pores, the trend was the opposite. Table 1 also shows that the D–P friction and D–P cohesion significantly decreased when the bitumen was drained. Especially, the D–P cohesion approached about zero, which meant the Karamay oil sands were similar to the sandy soil. The increase in temperature made the cap yield surface shrink because the cap parameter p_a increased.

Table 2. Young’s modulus and Poisson’s ratio of cap/base rock and screen [2,3].

Elastic Properties	Cap/Base Rock	Screen
Young’s modulus (GPa)	1.667	20
Poisson’s ratio (-)	0.2	0.2

shows Young’s modulus and Poisson’s ratio of cap/base rock and screen.

The eccentricity of the D–P plastic cap R is 0.3. The initial yield surface position ${\epsilon }_{\mathrm{vol}}^{\mathrm{in}}\left|{}_0\right.$ is zero. The shape factor of the transition zone α is 0.03. The ratio of triaxial tensile strength and triaxial compression strength k is 1.

3.1.3. Flow Parameters

The pore pressure diffusion across the reservoir can be described as a combined effect of the thermal expansion of the pore fluids, the compression of the pore fluid, and the volumetric deformation of the matrix skeleton [25,26]. According to the permeability experiments of Karamay oil sands at room temperature to 100 °C, the temperature-related Kozeny–Poiseuille model is used to consider the permeability improvement in the volumetric dilation process [27]. The relation between the permeability and volumetric strain is:

$$\frac{k}{k_0}=\frac{{\left(1+\frac{{\epsilon }_{\mathrm{v}}}{{\varphi }_0}-\frac{{\alpha }_{\mathrm{s}}\Delta T\left(1-{\varphi }_0\right)}{{\varphi }_0}\right)}^3}{1+{\epsilon }_{\mathrm{v}}}$$

(15)

where k is the dilation-induced permeability; k₀ is the initial permeability; ${\varphi }_0$ is the initial porosity; α_s is the thermal expansion coefficient; ΔT is the temperature change.

The reference permeability in Equation (15) at 20 °C and zero volumetric strain is 2.886 mD. The permeability at any temperature and volumetric strain can be calculated by Equation (15). The occurrence of volumetric strain is caused by two reasons: first is the shear dilation of oil sand skeleton under shear stress, and second is the tensile parting expansion induced by the decrease of mean effective stress under water injection. As shown in Figure 4, the dynamic viscosity of Karamay bitumen varies with the temperature.

The water flow capability is quantitatively described by the permeability coefficient:

$$\kappa =\frac{k{\rho }_{\mathrm{L}}g}{{\mu }_{\mathrm{L}}}$$

(16)

where κ is the permeability coefficient; ρ is the liquid density; g is the gravitational acceleration; μ_L is the dynamic viscosity of fluids.

As shown in Figure 5, the permeability coefficient evolved with the volumetric strain or void ratio and temperature in the process of water injection. When the temperature was high, the liquid viscosity became low, therefore, the permeability coefficient was high, which was beneficial to the hot water flowing into the far reservoir. The injected water in the oil field was a by-product of water with a temperature of 70~80 °C [2,3]. The water temperature was not specially heated, but was the natural temperature after treatment. The water used for injection during micro-fracturing was the by-product fluid recycled from previous recovery activities. The treatment mainly contained sedimentation and filtration. The viscosity of water under 70~80 °C was low, so it can maintain a high permeability coefficient. At present, the temperature for water injection is usually between 20~80 °C, and a good effect is proven by field practice.

3.1.4. Heat Transfer Parameters

Table 3. Heat transfer parameters of oil sands, mudstones, screens, and water (room temperature) [25].

Heat Transfer Parameters	Oil Sands	Mudstones	Screens	Water
Thermal conductivity (W/(m·s))	1.0	2.0	16	0.599
Thermal expansion coefficient (°C−1)	1.0 × 10−5	1.0 × 10−5	1.2 × 10−5	2.08 × 10−4
Specific heat capacity (J/(kg·°C))	882.7	2000	460	4200
Density (kg/m3)	1960	2500	7900	1000

shows the heat transfer parameters of oil sands, mudstones, screens, and water at room temperature.

In the partially drained zone, the bitumen started flowing at 56 °C, and is completely mobile at 70 °C. The phase change-induced latent heat is 266 J/g [22]. The average porosity of Karamay oil sands is 0.33, and the oil saturation is 42%, so the equivalent latent heat of the oil sand reservoir is 36.8676 J/g. The initial temperature of the reservoir is 12 °C.

3.2. Simulation Process

The Abaqus element code was used to solve the coupled geomechanical-thermal problem. As shown in Figure 3, the geometric model established consisted of four parts, namely caprock/base rock, reservoir, screen, and the wellbore drilled. The numerical calculation contained three steps, namely, the initial in situ stress equilibrium, drilling, and water injection (micro-fracturing).

The in situ stress test revealed the initial in situ stress conditions: the minimum horizontal effective principal stress on the z-axis is 0.9412 MPa; the maximum horizontal effective principal stress on the x-axis is 3.2218 MPa; the vertical effective principal stress is 3.982 MPa. The overburden pressure applied on the top of the caprock is 7.602 MPa. The boundary conditions were set as follows: the bottom is zero freedom; four sides are restricted in a normal direction; the top is free.

The initial formation pore pressure is 3.62 MPa. The initial void ratios of mudstone and reservoir are 0.05 and 0.2, respectively. The permeability of mudstone is 2 × 10⁻³ mD, and the permeability of the screen is set as 20 D [2,3].

Both the density and viscosity of oil and water dropped with the temperature increase.

Table 4. Critical parameters of water under different temperatures [28].

Critical Parameters	10 °C	80 °C
Density (kg/m3)	999.7	971.8
Viscosity (mPa·s)	1.3077	0.3565
Specific heat capacity (J/(kg·°C))	4178	4195

shows the critical parameter of water under different temperatures. The impact of temperature-induced viscosity change on the permeability coefficient must be considered. The thermal conductivity of water-saturated oil sand at different temperatures is about 1.8 W/(m·s) [24].

To simulate the coupled thermal-fluid-solid response of water injection, the type of element C3D8RPT, a 3D 8-node thermal-hydro-mechanical coupled reduced integral element, was used for all elements. For a higher calculation accuracy, a better convergence, and a smaller calculation quantity, the dense elements were adopted in the region around the SAGD wells (the smallest step length is set as about 4.55 cm), and the sparse elements were employed in the far reservoir region. The number of nodes is 17039, and the number of elements is 14970. The grid division result is shown in Figure 6.

The field monitoring data revealed that there are three phases in the reservoir reconstruction process, namely, well washing, water injection (micro-fracturing), and connectivity test. The simulation in the study corresponded to the micro-fracturing phase, where the 80 °C hot water was injected into the SAGD wells to improve the bottom hole pressure step by step. According to the wellhead flow and bottom hole pressure data, the micro-fracturing process for Well A could be simplified as two stages, in which 182 m³ water was injected during 58 h (2.5 days). During 0~23 h, the I well and P well were injected at the rates of 47.2 and 54.4 m³/day, respectively. During 24~35 h, the I well and P well were injected at the rates of 113.6 and 71.2 m³/day, respectively.

At present, only the bottom hole pressure (BHP) data in Well A can be monitored, and the reservoir temperature, in situ stress, and deformation can be difficult to acquire because of these reasons, such as equipment and technology. To prove the model reliability, the BHPs predicted by the numerical model were compared with the field monitored data (Figure 7). The comparison revealed that the predicted results agreed well with the monitored data, which means a good model reliability.

4. Numerical Simulation Results and Analysis

4.1. Temperature and Geomechanical Zones

Figure 8 shows the reservoir temperature distribution after water injection. The reservoir in the well-around region is heated by both heat conduction and heat convection effects. However, there were only two narrow separate zones in which temperature changed effectively because of the weak heat convection induced by a low permeability to water (<3 mD) [14], relatively low injection temperature (about 80 °C), and a short injection period (<3 days).

Figure 9 is the bitumen state distribution in the x-y plane around the I well after water injection. It is estimated that the geomechanically drained zone (where bitumen is liquid state) and partially drained zone (where bitumen is the mixture of liquid and solid states) show the widths of 10 cm and 9 cm, respectively, according to the proportion of the geometric model size and the actual formation size. During water injection, the oil sand skeleton underwent a process of effective stress unloading due to the great pore pressure increase and relatively low thermal expansion stress, therefore, the pore space would not collapse even when bitumen was mobile in the drained zone. The term “pore collapse” is used to describe a drastic, irreversible reduction in porosity due to the increase in effective stress on the rock matrix as a consequence of the withdrawal of oil and gas [29]. However, even this narrow drained zone may induce a significant effect on the subsequent cyclic preheating process, because the steam can invade the drained zone and form an initial steam chamber, which quickens the cyclic preheating speed and shortens the preheating period. Duong et al. [30] and Duong [31] observed this phenomenon by field test and proposed a concept named “heat ring” to quantitatively describe the steam chamber area where the perfect heat convection occurs in the preheating stage.

4.2. Pore Pressure

As shown in Figure 10, the BHP increased to 5.6 MPa from 3.62 MPa (the initial reservoir pore pressure). The BHP around the I well was slightly higher than that around the P well, because the injection rate in the I well was higher. It is worth noting that the reservoir pore pressure after operation must be about 0.5 MPa lower than the reservoir cracking pressure to prevent the two wells from inter-well breakdown [32]. The pore pressure of mudstone base rock right under the P well increased significantly, which illustrates that the mudstone interlayers or stringers may be broken through by the water injection operation near the mudstone. Considering the interferences induced by both two wells and mudstone caprock/base rock, the reservoir pore pressure isoline exhibited a shape of semi-ellipse that is symmetric to the mudstone base rock.

4.3. Void Ratio

The void ratio is defined as the ratio of pore volume and solid particle volume. As shown in Figure 11, the void ratio in the reservoir had a small improvement after water injection, especially in the inter-well zone, where the volumetric strain changed from 0 to 0.167%, and the void ratio changed from 0.2 to 0.202. The well-around zone owned a relatively smaller void ratio than the inter-well zone. Generally, the void ratio improvement zone induced by water injection showed a shape of a large quasi semi-ellipse. Compared with the study of Lin et al. [2], the largest volumetric dilation did not occur on the well wall if the transient heat transfer and bitumen state around the wellbore were incorporated. The increase in void ratio meant the increase of both reservoir porosity and permeability, which provided the advantageous flow channel for steam in the subsequent cyclic production and preheating phases.

The next step was to consider the influence of the fracturing pressure along the horizontal wellbore on reservoir dilation, aiming to find a critical fracturing pressure in the heterogeneous reservoir for simultaneous operation, and further provide the guidance for the well selection and construction in the oilfield.

4.4. In Situ Stress

As shown in Figure 12, the reservoir maximum principal stress changed from 3.2 MPa everywhere to uneven distribution (tension stress is positive). The maximum principal stress in the reservoir near the wellbore was lower (the lowest value was about 2.0 MPa), while the maximum principal stress in the far reservoir stayed unchanged. After drilling the borehole, the redistribution of the stress around the well occurred. The impact of stress redistribution was small, and the drilling-induced maximum principal stress exhibited 3.1~3.4 MPa.

Figure 13 shows the directions of the reservoir’s maximum principal stress after water injection. It could be known that injected water disturbed the direction of the maximum principal stress in the inter-well region and well-around region, from the direction perpendicular to the wellbore to a complex state. Due to the artificial fracture which developed along the direction of the maximum principal stress, the water injection-induced changes in the maximum principal stress magnitude and direction will influence the generation and propagation of the artificial fracture.

As shown in Figure 14, the reservoir vertical stress changed from 3.98 MPa everywhere to uneven distribution. The well-around region possessed lower vertical stress (as low as 2.54 MPa), while the far reservoir had a higher value.

Figure 15 shows the reservoir minimum principal stress changed from 0.94 MPa everywhere to uneven distribution, with the lowest value of 0.46 MPa taking place in the inter-well region and the highest value of 1.6 MPa occurring on the well wall, resulting from water injection-induced coupled thermo-hydro-mechanical process. The relative magnitude of the maximum/minimum/vertical stresses will impact the propagation and pattern of artificial fracture, which is one of the most important issues in the process of reservoir stimulation for engineers.

4.5. Displacement

As shown in Figure 16, the reservoir uplift is caused by water injection. The vertical displacement in the plane perpendicular to the wellbore is symmetrical about the connection between the two wells. The largest displacement reached about 5 cm, and the largest uplift took place in the zone right above the two wells. The displacement along the wellbore was very small, and it can be neglected.

The field monitored data shows that there is a slight uplift on the ground after water injection. To observe the reservoir uplift information at different depths, a vertical line was established through the two wells from the reservoir top to bottom. Figure 17 displays the vertical displacement along the vertical line. It can be observed that the uplift magnitude increased sharply from the reservoir bottom to the middle of the reservoir, then the uplift magnitude went up slowly from the middle of the reservoir. The phenomenon of reservoir uplift verified the mechanism of water injection-induced dilation. The uplift was relatively small because of the short injection period and low water temperature. In the next study, the ground uplift magnitude and caprock stability will be focused on under the condition of long-term high-temperature steam injection.

5. Discussion

In this paper, a case study of the Well A in the Z1 block in Xinjiang Fengcheng Oilfield was numerically simulated to show the geomechanically coupled responses (water injection-induced temperature, pore pressure, stress redistributions, and formation uplift) for the approximately homogeneous reservoir under hot water injection. Water injection is a process of the water being squeezed into the pore pace below the formation, fracturing pressure for 3~4 days. The area where more water was absorbed led to a larger volumetric strain and a higher pore pressure. The pore pressures in the reservoir around wells and the base rock right below the production well were increased significantly, due to the interference effects of SAGD well pairs, mudstone cap, and base rocks [33].

The calculation results induced by the sensitivity of the grid (e.g., density, element type) were acceptable in our research group’s previous study [2,3,22,25]. In order to simulate the coupled thermal-fluid-solid processes, the type of element was chosen as C3D8RPT (the eight-node brick, tri-linear pore pressure, tri-linear displacement, and reduced integration). The grids near the I and P wells were intensive for higher precision and a smaller computation quantity. In reality, when we evaluated the process of water injection inconsiderate to temperature, no difference was found between their simulated results whether using C3D8P or C3D20RP. In the next step, we intend to re-establish the relation between the predicted results (e.g., pore pressure) and the grid number through some calculation examples. The validation of our numerical model can be also verified by comparing the predicted BHP (bottom hole pressure) results with the actual test data, just as shown in Figure 7.

According to the difference in fluid flow, deformation, and temperature propagation properties, the geometric model was partitioned as 5 zones, namely caprock, base rock, pay zone, screen tube, and the wellbore drilled. Only in this way, can we assign different material property parameters for the five meshed zones. Particularly, the reservoir can be further partitioned as 3 zones, namely undrained, partially drained, and drained zones, considering that the geo-mechanical behavior in the process of SAGD is significantly affected by the phase transition of bitumen. The seepage, deformation, and temperature propagation behaviors within the oil sand reservoir exhibit big differences for the geomechanically undrained, partially drained, and drained areas, in which the “phase states” or bitumen viscosity varies a lot [23,25]. We used the temperature-dependent flow and heat transfer parameters to embody this process in the reservoir. In the z-direction, all the zones have the same thickness, which means the geometric model’s division is generated by sweeping along the z-direction.

In our study, good model reliability means that the predicted BHPs are just essentially in agreement with the field recorded BHP values. The predicted results are very difficult to quite match the actual data, because of these factors such as the subjective and objective fluctuations in field operations, the heterogeneity of the reservoir’s material, the assumptions and simplification of the numerical models, and the boundary conditions, and so on. For the field engineers, they desire the predicted results within an acceptable range. As shown in Figure 7, the field recorded BHP values are relatively irregular including some dips and scatters, which may result from the fluctuation of wellhead operation pressure. Besides, the raised BHP may cause fracturing of certain thin lenses, generating preferential flow paths in the reservoir and connecting natural fractures or weak planes [2], which may lead to a sudden drop in BHP. In addition, there are some dips and scatters during 35~58 h because we used a simplified constant injection rate for the simulation. For instance, the I well and P well were injected at the constant average rates of 113.6 and 71.2 m³/day during 35~58 h, respectively. If we divide 35~58 h into several or a dozen periods, the matching degree may be better theoretically, but the workload will increase and it is not necessary for engineers.

At present, however, it is difficult for us to quantitatively give the level of uncertainty. There are several uncertainties concerning the nature of the reservoir rock itself in the prediction of injection. For instance, both compositional and mechanical heterogeneities are significant in the oil sand reservoirs, creating a huge obstacle for accurately assigning proper constitutive properties for the analyzed reservoir. Another noticeable uncertainty lies in the possible existence of caverns and natural fractures in the reservoir. In addition, the formation damage introduced in drilling and completion and by the impurities in the injected water has not yet been studied, whereas its importance in affecting injection is worth further investigation [2].

In our study, it is assumed that the reservoir is completely saturated in the water phase during injection, for which the effective stress concept was adopted [2]. It must be aware that for oil sand reservoirs that are usually in a depth of less than 650 m, the pores may not be fully saturated. In a further study, this challenging task incorporating a multi-phase seepage process can be attempted if we have access to measuring the gas-water two-phase permeability of our oil sand specimens in the laboratory.

In general, laboratory and field experiments should be performed prior to numerical simulations, as experiments can be used to find precise mechanisms, quantitative engineering, and model parameters for subsequent numerical simulations. Since there are no pilot tests, the key simulation parameters used are undetermined in the numerical calculation. Besides, the validity of the simulated results cannot be evaluated. Nevertheless, if sufficient experimental data become available, numerical simulations can be extensively used for the evaluation of the micro-fracturing effects using different site operation parameters. Computational simulations could be employed for the optimization of operating strategies and offer engineers data visualization. In conclusion, experiments are more accurate in the early stages of oil sand development, but numerical simulations are more widely used in subsequent stages [25].

The water injection process may cause water convection and multi-phase seepage near the wells. In the next step, we attempt to combine multi-phase seepage and convection by coupled thermal-fluid-solid methods. Indeed, these simplifications in the simulation are rational and credible for the field issues of oil sands’ exploitation under micro-fracturing. The hypotheses must be redetermined carefully, depending on the specific reservoir and recovery method [25]. For example, cyclic injection types with low reservoir permeability and relatively low injection pressures can make pore fluid flow behavior and its induced thermal convection less critical.

The phase transition of bitumen has an important impact on the geomechanics concerning fluid seepage, pay zone deformation, and temperature propagation behavior during SAGD. The seepage, deformation, and temperature propagation behaviors within the oil sand reservoir exhibited large differences for the geomechanically undrained, partially drained, and drained areas, in which the “phase states” or bitumen viscosity varied a lot [23,25]. Hence, the bitumen’s “phase transition” was taken into account in this paper.

6. Conclusions

This study predicted the underground geomechanical-thermal responses under water injection such as formation pressure, temperature, and in situ stress, as well as displacement, using the coupled thermo-hydro-mechanical element model incorporating bitumen’s state and cap plasticity. Several conclusions can be drawn as follows:

After water injection, the reservoir temperature changed significantly only in two independent areas around the well. The large temperature gradient and thermal expansion stress can lead to the evolution of in situ stress. In the geomechanical drained zone, the bitumen in a small area around the well was completely melted. The width of the partially drained zone was similar to that of the drained zone. The existence of the geomechanical drained zone and partially drained zone expanded the heating ring in the cyclic preheating stage indirectly.

Due to the transient non-uniform distributions of reservoir temperature and pore pressure, the magnitude and orientation of reservoir maximum/minimum principal stress and vertical stress had complex changes, which can influence the fracture propagation trend and pattern. The reservoir pore pressure had a significant improvement, and the isoline exhibited a shape of semi-ellipse that was symmetric about the mudstone base rock.

The reservoir was lifted as a whole, and the maximum uplift area was located in the right upper part of the two wells. The uplift magnitude increased sharply from the reservoir bottom to the middle of the reservoir, then the uplift magnitude went up slowly from the middle of the reservoir. The void ratio in the reservoir had a large improvement after water injection, especially in the inter-well zone.