Porothermoelastic response of an oil sand formation subjected to injection and m

Botao Lin · Han Meng · Jingjun Pan · Sen Chen

Abstract Stimulation of unconsolidated formations via horizontal wells has seen its vast implementation in the recent development of heavy oil reservoir to save the time and cost of preheating the reservoir before the steam-assisted gravity drainage (SAGD)process. A mathematical approach was proposed in this research that fully couples the hydraulic, mechanical and thermal responses of unconsolidated sandstone formations and also applies failure criteria for describing either shear dilation or tensile parting mechanism that generates microcracks. The approach was implemented to predict the porothermoelastic response of a pair of SAGD wells subject to injection and subsequent micro-fracturing using hot water. It was found that the predicted bottom hole pressures (BHPs) match closely with the field observed data. An elliptical dilation zone developed around the dual wells with relatively high pore pressure, porosity, permeability and temperature, implying good interwell hydraulic communication between both wells. The activation of microcracks dramatically accelerated the dissipation of pore pressure across the entire formation depth and also facilitated heat convection in between the dual wells, though to a lesser extent. In summary, the approach provides a convenient means to assist field engineers in the optimization of injection efficiency and evaluation of interference among multiple horizontal wells.

Keywords Porothermoelastic response · Unconsolidated sandstone · Injection · Micro-fracturing · Microcracks · Oil sand

1 Introduction

Unconsolidated sandstone reservoirs are widely distributed over the world, bringing in specialized geomechanical treatment or development methods. McLarty et al. (1993) provided an overview of the problems and techniques related to the application of horizontal wells in the drilling and completion of unconsolidated sandstone formations in the Gulf of Mexico. Rabaa Ali et al. (2009a, b) applied rock mechanics methods to analyze well stability during drilling and long-term screen integrity considering in situ stress field for an unconsolidated sandstone reservoir in Saudi Arabia;they also used a multi-disciplinary approach through maximizing reservoir contact to enhance sweep efficiency and increase productivity. Robles et al. (2012) designed slurries containing polymeric admixtures to solve the near wellbore channeling caused by sand production at a heavy oil field in the Neuquina basin, Argentina. Yuan et al. (2011)applied micro-fracturing to the heavy oil sand reservoirs in Alberta, Canada, via dual horizontal wells, so as to achieve dilation of the near wellbore region for fast start-up of a steam-assisted gravity drainage (SAGD) project. Most of these research endeavors dealt with the use or analysis of geomechanics in horizontal wells. It is therefore desirable to learn poromechanics responses of horizontal wellbores in an unconsolidated sandstone formation to assist stimulation or development practices.

Analytical solutions for poroelastic or porothermoelastic responses of a vertical or inclined borehole have been proposed by a number of researchers, including Rice and Cleary(1976), Detournay and Cheng (1988), Cui et al. (1999),Abousleiman and Ekbote (2005), Abousleiman and Chen(2010), Chen and Abousleiman (2016), etc. Abousleiman and Ekbote (2005) formulated porothermoelastic solutions for an inclined borehole in a transversely isotropic medium with the axis of the borehole being orthogonal to the bedding plane of the formation. Ekbote and Abousleiman (2005,2006) proposed poroelastic solutions for an inclined borehole that couple both chemical osmosis and thermal conduction. The above solutions can readily be extended to analyze the porothermoelastic response of horizontal wells, provided that the stress and boundary conditions stay unchanged.Although the non-hydrostatic stress field was considered mainly in these literature works, far- field in situ stresses and in finite boundary conditions were usually assumed. However, they are not able to account for the finite boundary conditions. Also, the permeability has been assumed as constant in deriving analytical solutions (Detournay and Cheng 1988; Abousleiman and Ekbote 2005; Abousleiman and Chen 2010). Such an assumption suffices for analyzing the poroelastic response of boreholes in formations composed of relatively stiff rock such as shale or limestone, but an evolving permeability is essential for describing the dilation behavior of weakly cemented unconsolidated formations.In order to account for finite boundaries and permeability evolution, numerical simulations using either the finite element or finite difference method have been developed to analyze the hydromechanical responses of unconsolidated formations (Pak and Chan 1996; Settari and Walters 2001;Lee and Ghassemi 2010; Ruan et al. 2012). Especially, Xu and Wong (2010, 2013) and Lin et al. (2017a, b) used elastoplastic constitutive models and the finite element method to evaluate the hydromechanical response of an unconsolidated sandstone formation through injection in single or dual horizontal boreholes. It is worth mentioning that for the stimulation of an unconsolidated sandstone formation,the field injection rates are unlikely to trigger plastic deformation in formations buried at a depth of 450 m or deeper(Lin et al. 2017a); in that case, poroelasticity suffices for the evaluation of the hydromechanical behavior of an unconsolidated formation. In the aforementioned numerical efforts,the thermal effect was usually not involved, mostly attributed to the belief that the in fluence of the injected water temperature was negligible on the hydromechanical behavior of an unconsolidated sandstone formation. This may not be true in certain unconsolidated formation such as heavy oil sands,whose bitumen content has its viscosity being extremely sensitive to temperature. Moreover, when the injection pressure rises to certain level, the formation will break down in a similar manner to that of tight formations, generating microcracks instead of linear apertures at the region in proximity to the boreholes (Kry et al. 1992; Collins 2007; Yuan et al.2011; Lin et al. 2016). In literature, some researchers have integrated damage mechanics to the continuum mechanics regime for the simulation of water injection or fracturing in tight formations by taking the permeability as a function of damage variables (Tang et al. 2002; Selvadurai 2004;Lee and Ghassemi 2010). However, no research work has yet been conducted to investigate the numerical simulation about the in fluences of microcracks on the hydromechanical behavior of unconsolidated sandstone formations.

Given such a situation, this study proposes an efficient numerical solution to the porothermoelastic response of a pair of horizontal wellbores subjected to both hot water injection and subsequent micro-fracturing in an unconsolidated sandstone formation, which is sandwiched in between impermeable mudstone barriers. The engineering goal is to bring in controlled dilation of the formation region surrounding the dual horizontal wells, so as to establish interwell hydraulic communication, de fined as the condition that the pressure of one well is quickly responsive to the other well if closing the former and varying the pressure of the latter (Yuan et al. 2011; Dragani and Drover 2016). Consequently, the time for the formation of a steam chamber around the dual wells can be dramatically shortened from more than a year to less than two months, such that not only economic and time costs can be substantially reduced, but also the steam chamber has better conformance along the well length in the later steam circulation of a SAGD project.The solution proposed in the study treats the formation as a porous continuum medium, part of which experiences an abrupt jump in permeability when it is failed either by shear or tensile stresses (with microcracks generated). The solution is primarily featured by the following three characteristics: (1) the thermal, hydraulic and mechanical behavior is fully coupled. (2) The permeability evolves not only with deformation but also temperature changes. (3) The development of micro-fracturing and its in fluences on the porothermoelastic responses of the formation are considered in the simulation.

2 Governing equations for the porothermoelastic responses

For a plane strain injection problem in an in finite porous medium, the governing equations for the deformation of a porothermoelastic isotropic medium under the Cartesian coordinate system can be derived as follows (Biot 1941);Rice and Cleary 1976; Detournay and Cheng 1988; Selvadurai et al. 2017):

whereνandνuare the drained and undrained Poisson’s ratios;uidenotes the displacement along each coordinate axis, m;εvis the increment of fluid content per unit volume;evis the volumetric dilation of the solid matrix per unit volume;αsis the linear thermal expansion coefficient of the solid phase, K−1;βsfis the thermic coefficient of the pore fluid, K−1;Bis Skempton’s pore pressure coefficient;Tis the current temperature, K;∇2denotes Laplacian operator which is given by

The volumetric dilation of the solid matrix per unit volumeevcan be written as:

In Eqs. (2) and (3), the summation rules apply for the repeated indices. For a plane strain problem (on anx–zplane) in an isotropic homogeneous medium, the items containing derivatives with respect toycan be dropped off.Note that the tensile stress and compressive pore pressure are taken as positive, whereas the compressive stress is negative in the present analysis. For a linear thermal expansion without thermal shear strains, the constitutive relations can be expressed as (Abousleiman and Ekbote 2005):

whereGis the shear modulus, Pa;σijdenotes the total stress tensor, Pa;eijdenotes the strain tensor;pis the pore water pressure, Pa;δijis the Kronecker delta function;βsis the thermic coefficients of the solid skeleton, Pa K−1; ΔTis the change from the reference valueT0where no thermal strain exists, K;αandMare the Biot effective stress coefficient and Biot modulus (Pa) given by

The thermic coefficients of the solid skeleton and the pore fluid,βsandβsf, can be expressed as follows:in whichαfis the volumetric expansion coefficient of the pore fluid, K−1;ϕis the porosity. Following the literature(Abousleiman and Ekbote 2005; Cheng 2016), while adding in the effect of heat convection (Bear and Corapcioglu 1981;Zhao et al. 2014; Cheng 2016), the fluid and heat diffusion equations are given by:

wheretis time, s;cfis the fluid diffusivity, m2s−1;chis the bulk heat diffusivity,is the relative speci fic heat capacity of the fluid to that of the bulk per unit volume.chandare given by

hereκis the mobility, m2Pa−1s−1, and de fined byκ=k/μ, in whichkdenotes the intrinsic permeability, m2andμthe fluid viscosity, Pa s;λis the bulk thermal conductivity, W m−1K−1;ρfandρare the densities of the pore fluid and the bulk,kg m−3;andCvdenote the speci fic heat capacity per unit mass of the pore fluid and of the bulk, J K−1kg−1, respectively.cfandcpare the coupling constants in accordance with

The last term in Eq. (12) denotes the heat transport by fluid convection. BothλandρCvcan be obtained as the weighted average of the properties of the solid and fluid phases (Bear and Corapcioglu 1981; McTigue 1986):

in which the superscripts s and f denote the solid and fluid,respectively. The transport equations follow Darcy’s law and Fourier’s law:

whereqis the fluid flow velocity, m s−1;his the heat flow rate, J s−1. Both flow rates are assumed to be isotropic for simplicity at present. The evolution of permeability for unconsolidated sandstones because of matrix dilation can be achieved using the Kozeny–Poiseuille correlation (Tortike and Farouq Ali 1993; Wang and Xue 2002; Du and Wong 2007):

wherek0is the initial permeability, m2;ϕ0the initial porosity. Because the dynamic viscosity decreases as a function of temperature as:

wherea= 2.417 × 10−5Pa s,b= 247.8 K andc= 140 K provide the best fit for the viscosityμcurve of water with respect to temperatureT. Therefore, the mobility evolution equation can be derived based on Eqs. (21) and (22) as:

As illustrated above, the porothermoelastic solution considers heat transfer due to conduction through the matrix (Eqs. (12) and (20)), convection via the porous fluid(Eq. (12)) and mobility evolution as a function of mechanical volumetric strain, thermal expansion, accompanied by changes in water viscosity because of temperature variation(Eqs. (21)–(23)).

3 Simulation of micro‑fracturing

In the injection phase, the BHPs have to be retained below the breakdown pressure of the unconsolidated formation,opposite to many other reservoir stimulation practices that desire as a high degree of fracture network as possible. A continuing rise in the injection pressure will finally lead to micro-fracturing of the unconsolidated sandstone formation neighboring the boreholes induced by shear dilation or tensile parting (Agar et al. 1987; Oldakowski 1994; Samieh 1995; Wong 1999; de Pater and Dong 2007). Shear or tensile microcracks are generated instead of large linear apertures as in tight formations (Kry et al. 1992; Yuan et al. 2011; Lin et al. 2016). Figure 1 shows the scanning electron micrograph taken on the studied formation sample after shear at a con fining stress of 5.5 MPa and a pore pressure of 5 MPa in a standard compression test. The microcracks and microvoids shown in Fig. 1 are described as microcracks in general in the following discussion.

Fig. 1 Micrograph of Karamay oil sand after triaxial shear

If the microcracks are created by the shear dilation mechanism, the Mohr–Coulomb failure criterion can be used to trace this phenomenon:

whereτfis the shear stress at failure;σnis the normal stress on the failure plane;cis the cohesion; andφis the friction angle. The Coulomb failure criterion can be rewritten in terms of the maximumσ1and minimumσ3effective stresses:

If the microcracks are generated by the tensile parting mechanism, the tensile stress is considered to exceed the inherent tensile strength of the matrix. For the borehole parallel to the orientation ofσh(just as in this case), the fracture initiation pressurepffor a tight formation can be derived as(Fjar et al. 2008):

which was derived from analysis of tensile failure at the horizontal borehole. Here,Stis the tensile strength of the matrix, Pa. However, Eq. (26) is not applicable any more for unconsolidated sandstones, especially an oil sand formation, because microcracks of irregular patterns rather other than clear tensile aperture at a speci fied stress orientation occur in such a scenario. In this respect, a negative effective mean stress with its absolute value exceeding the tensile strength of the matrix is taken to be the stress condition that triggers micro-fracturing (here, the sign convention of rock mechanics is followed, i.e., compressive stress is positive,and tensile stress is negative):

whereσcis the mean con fining stress of the formation, Pa. In Eqs. (25) and (27),F1andF2can be viewed as the shear or tensile failure index. When eitherF1orF2at a location turns to be greater than or equal to zero, microcracks are considered to form. It was found that for the studied Karamay oil sand, the permeabilityksharply increases by a factor of six after microcracks have formed (Lin et al. 2016; Chen et al.2018). Such a permeability increase was implemented in the simulation on the effects of microcracks.

As microcracks propagate from a borehole to the deeper part of the formation, the pore pressure inside the failed region quickly dissipates because of an abrupt increase in permeability, bringing the region back to the porothermoelastic regime. Therefore, the porothermoelastic solution(Eqs. (1)–(23)) and the criteria for microcrack generation(Eqs. (25) and (27)) combine to provide an efficient method to simulate both injection and micro-fracturing in an unconsolidated sandstone formation. If the injection continues,then the pressure will rise again until more microcracks are generated, followed by another pressure drop.

4 Formulation of the problem and finite difference solutions

4.1 Formulation of the problem

In this study, the hot water stimulation in a pair of SAGD wells placed in a heavy oil sand reservoir (a common type of the unconsolidated sandstone reservoir) was proposed as an illustrative example. The geometry of the plane strain problem for analysis is demonstrated in Fig. 2. Two parallel horizontal wells marked as I and P wells are completed in a heavy oil sand formation. The formation is sandwiched in between nearly impermeable mudstone layers. The mechanical, hydraulic and thermal responses of the formation are assumed to be isotropic. The far field in situ stresses are adopted as the initial stress states. Zero normal displacement conditions are imposed on all of the boundaries of the formation except the top surface on which a free vertical displacement is allowed under an overburden stress. Water at a temperature beyond that of the formation is injected via the two horizontally aligned boreholes of a SAGD wellpair. No seepage is allowed at the top and bottom surface to simulate the flow barriers made by the impervious mudstone layers.Also, zero seepage flow condition is imposed on the lateral boundaries far enough from the injection wells such that the pressure and stress disturbances caused by injection would dissipate at that distance.

For an injection project taking place in a heavy oil sand reservoir, the 300–500 m long horizontal SAGD wells were usually completed with screen pipes. In a stimulation project on such a weakly cemented reservoir, the fluid flow through the pay zone and its induced formation dilation are more of a concern, whereas the wellbore stability needs little consideration. Consequently, the boreholes can be simpli fied as two nodes in a rectangular mesh. A finite difference scheme can be adopted to deal with the simple geometry and tackle the top, bottom and lateral boundary conditions as demonstrated in Fig. 3 below.

Fig. 2 Geometry of the dual horizontal wells in an oil sand reservoir(not to scale)

Fig. 3 Schematic diagram of a finite rectangular poroelastic medium subjected to dual well injection

The flow rates at the I and P wells will beqIandqP, whereqI=QI∕LandqP=QP∕L;QIandQPare the fluid flux per unit time through the I and P wells, respectively;Lis the length of the dual wells. The heat fluxes are introduced ashIandhP, which are the flow rates multiplied by the temperature change, the fluid density and the speci fic heat capacity of the fluid. Since the SAGD wells were completed with screen pipes, no restraint is applied to the fluid and heat fluxes via the boreholes.

4.2 Numerical scheme

The mechanical equilibrium equations (Eq. (1)), the mass conservation equation for the pore fluid (Eqs. (10) and (11))and the energy balance equation (Eq. (12)) constitute a system of coupled partial differential equations. The in finite difference method was adopted to solve these coupled equations, in which the central difference approximation was used for the second-order spatial derivatives while a backward difference was used for the first-order derivatives with respect to time, following the finite difference scheme used by Liang et al. (2010), Ruan et al. (2012) and Rahmati et al.(2017). First, the pressure and temperature were derived implicitly using Eqs. (5) and (10)–(12) ( flow module):

Figure 4 reveals that in the flow module, the pressure and temperature were first solved using Eqs. (10)–(12). Secondly, in the geomechanics module, the updated pressure and temperature values were adopted to achieve the displacements of the formation through Eqs. (1)–(5). At each time interval, these two steps were iterated until the convergence was reached to meet the prescribed accuracy. The porothermoelastic solution and the associated numerical approach were integrated into a self-developed computer code, which is veri fied in the following section through comparison on the predicted bottom hole pressures (BHPs) of a SAGD wellpair with the field recorded data. The predicted BHPs of the SAGD wellpair were compared among various grid sizes dx(dx= dz, c.f. Fig. 3) and initial time steps dt(if convergence is not satis fied, then dtis halved) in Fig. 5 to check for the stability of the numerical solution. The details of the input parameters will be presented in the next section. It is demonstrated by Fig. 5 that choices of dx= 0.2 m and dt= 15 s suffice for the numerical solution.

Fig. 4 Numerical scheme for solving the differential equations

Fig. 5 In fluences of grid size and initial time step on the predicted BHPs

5 Applications to a field stimulation project 5.1 Geotechnical and engineering background

The stimulation project took place in a heavy oil reservoir in the Junggar Basin north to the Karamay city of Xinjiang Uygur Autonomous Region, China. The oil sand formation has a thickness of 35 m, the top and bottom of which were measured at depths of 340 and 375 m, respectively, and is sandwiched by nearly impermeable mudstone formations.The SAGD wellpair was drilled to a depth of 366 m for the I well and 371 m for the P well, respectively, with the latter overlying 4 m above the bottom mudstone barrier. The horizontal length of the dual wells was 390 m from heel to toe.The plane strain model was constructed to be 100 m × 35 m stacked with 3.5 × 105square elements having a dimension of 0.1 m × 0.1 m, summing up a total number of 351,351 nodes. The length of the wellpair (L= 390 m) was used to calculate the implemented flow rates and heat fluxes per unit length.

The laboratory-tested material parameters of the selected oil sand formation are presented in Table 1 (data from Li et al. 2015; Lin et al. 2016, 2017b; Chen et al. 2018; Gao et al. 2018). The geostress gradients were measured to be 0.014, 0.017 and 0.021 MPa/m in the minimum horizontal(σh), maximum horizontal (σH) and vertical (σv) directions.Given that the wells were drilled mainly in the direction ofσh, the calculatedσHandσvwere implemented as the initial stress state conditions for the plane strain model depicted in Fig. 3.

Table 1 Material parameters of the oil sand formation

Fig. 6 Comparison between the predicted and the observed BHPs at the boreholes

The project was conducted in October when the local temperature stayed at 6–15 °C, while the formation temperature as monitored by the ten equally spaced thermocouples buried with the P well was about 20 °C before injection. The injected fluid was produced water from previous SAGD projects, which was measured at a temperature of 80 °C. The entire project lasted for 66 h, composing two injection stages and the last stage of micro-fracturing and shut-in. The field engineers applied a net flow rate of 10.7 m3/day for the I well and 12.5 m3/day for the P well during the first 25 h, 18.1 m3/day and 19.6 m3/day for the next 33 h, at the end of which the BHPs became close to (about 0.3–0.5 MPa below) their estimated formation breakdown pressurepf= 6.04 MPa based on Eq. (26), as illustrated by the scattered field data dots in Fig. 6. At that moment, the injection rates were ramped up to be 40.7 m3/day and 41.5 m3/day in the I and P wells in order to microfracture the formation. It was found that the breakdown pressures were not reached at 6.04 MPa as expected.Instead, the micro-fracturing initiated atpf= 6.57 MPa (I well) and 6.70 MPa (P well). The micro-fracturing job lasted only for 4 h to produce microcracks adequate to facilitate interwell (hydraulic) communication, but not run the risk of creating undesirable large linear apertures from the convergence of a large number of microcracks.The large fractures will perform as primary flow channels for steam circulation in the following SAGD process after the stimulation, during which the steam introduced to the I well will rapidly be produced through the P well and circulated to the ground, instead of transferring most of the latent heat to the formation. Noticing that both boreholes had experienced apparent pressure drops, the two wells were quickly shut into prevent excessive fracturing of the formation.

Fig. 7 Pore pressure pro files a after injection and b after micro-fracturing

Fig. 8 Porosity pro files a after injection and b after micro-fracturing

Fig. 9 Permeability pro files a after injection and b after micro-fracturing

5.2 Porothermoelastic responses of the formation

In Fig. 6, the predicted bottom hole pressures (BHPs) by the mathematical approach in this study are plotted as a function of injection time and compared to the observed BHPs in the field. In general, the simulated BHPs match the observed values quite well. Meanwhile, the approach predicts a breakdown pressure ofpf= 6.84 MPa in the P well and 6.81 MPa in the I well, slightly greater than the field recorded values (6.57 and 6.70 MPa). The comparison suggests that the proposed criteria for microcrack generation can properly describe the breakdown pressure in unconsolidated sandstone formations.

To evaluate the porothermoelastic responses of the formation after injection (at the end of the 2nd injection stage as denoted in Fig. 6) and those after micro-fracturing, the corresponding pro files of pore pressure, porosity, permeability and temperature are demonstrated in Figs. 7, 8, 9 and 10 using MATLAB (2016b).

It is shown in Fig. 7 that the pore pressure increases in an elliptical con figuration encompassing the dual wells. The pressure is highest at either the I or P well borehole, dissipating gradually deeper into the formation. A comparison between Fig. 7a and b discloses that the formation of microcracks has dramatically facilitated the pressure dissipation across the region surrounding the dual wells. An elliptical zone with a pressure about 5.0–5.2 MPa has formed around the wells, implying that the interwell hydraulic communication has been well established.

Figure 8a shows that the amount of porosity rises by 0.3% at the boreholes after injection (ϕ0= 15%, Table 1).Upon micro-fracturing, the enhanced unit pore volume at the boreholes decreases down to 15.15% because of pressure drawdown that leads to an increase in effective con fining stress (Fig. 7b). At the same time, the porosity magnitudes distribute more evenly around the dual wells after the microfracturing. It must be noted that these porosity increases are relatively small because only 89.7 m3of hot water was injected into the formation through two 390 m long horizontal wells during the entire stimulation process.

In Fig. 9, the permeability pro files demonstrate con figurations very similar to those of the porosity, except that two small deep-colored spots emerge at the borehole positions.These two spots denote the geometrical extension of the microcracks, exhibiting as circular shapes with a radius of 0.3 m. It is interesting to observe that a relatively narrow zone of microcracks has signi ficantly assisted in the dissipation of pore pressure around a borehole (c.f., Fig. 7b and Fig. 9b).

Fig. 10 Temperature pro files a after injection and b after micro-fracturing

Meanwhile, after injection, the temperature increases in a circular shape around each borehole, maximizing at 29.1 and 29.4 °C at the I and P well boreholes, respectively(Fig. 10a). After micro-fracturing, the two zones of temperature increase expand and merge, while the temperatures at the two boreholes reduce to 26.2 and 26.4 °C (Fig. 10b)because the microcracks have strengthened both hydraulic and thermal dissipation. It must be noted that the temperature displayed by Fig. 10 represents that of the bulk instead of the pore fluid alone. The amount of temperature increase also indicates that injecting hot water will not substantially decrease the viscosity of the bitumen (from 2.6 × 106mPa s at 20 °C to 4.8 × 105mPa s at 30 °C); in other words, the bitumen remains immobile.

5.3 Interwell hydraulic communication

As previously mentioned, the primary goal of this injection example was to establish an interwell communication of the I and P wells such that the subsequent steam circulation via both wells can create a steam chamber within a short period.The interwell communication has been illustrated as the formation of a region with enhanced pore pressure, porosity,permeability and temperature that embraces the dual wells(Figs. 7, 8, 9, 10). Such a communication can be veri fied through either pressure or thermal measurements in the field.After the micro-fracturing and shut-in, the field engineers waited for another three hours during which the BHPs further decreased to 4.94 and 4.78 MPa in the I and P wells,respectively. Afterward, they opened the P well and raised its BHP to 6.4 MPa within 20 min, then quickly stopped pumping and closed the P well again. The BHP of the P well then stayed approximately at 6.1 MPa for 1 h and 18 min. Although remaining closed, the I well recorded a slight BHP decrease of 40 kPa within 20 min, then a steady increase up to 94 kPa afterward (Fig. 11). The subsequent pressure growth indicates that the interwell communication had been accomplished. After two days passed, the engineers reopened the I well and started circulating hot water at around 80 °C to it while keeping the P well closed. The borehole fluid temperature in the P well had decreased from 76 to 65 °C because of thermal dissipation but gradually recovered by an increase of 2.4–3.3 °C during the next six-hour monitoring of temperatures at the ten thermocouples placed in the P well (Fig. 12).In the P well, the ten thermocouples were aligned at an interval of 43.3 m along the 390 m long horizontal well, with the first and the last being located at point A and point B of the horizontal section, respectively. The observation of the temperature evolution implies that the interwell communication has facilitated the heat transfer via the microcracks and expanded pores. In addition, the relatively consistent increase in temperature among the ten thermocouples again supports that the formation is homogeneous because the thermal responses are evenly distributed along the horizontal well length. In subsequence, the porothermoelastic behavior of the formation can be appropriately evaluated using the plane strain model.

Fig. 11 Pressure response of the I well upon hydraulic communication

Fig. 12 Temperature responses of the P well upon hydraulic communication

6 Discussion

6.1 Mechanical responses at the boreholes

It is desirable to learn the evolution of stress states at the two boreholes at three critical states (at the end of 2nd injection stage, at the peak BHP and after micro-fracturing, c.f., Fig. 6), so as to understand when and how the matrix at the borehole location failed to form microcracks.In this regard, the Mohr’s circles representing these stress states are presented in Fig. 13.

The black Mohr’s circles at the two boreholes that represent the original stress state first move to the left (the deep blue ones) because of the increased pore pressure caused by the injection, then further translate leftward and shrink in radii as the BHPs were raised further to initiate micro-fracturing (the red ones). After micro-fracturing,the circles move backward and stay as the purple ones,as the microcracks rapidly dissipated the BHPs. It is also observed that the Mohr’s circles in both boreholes (which are most likely to yield relative to other locations) did not reach the yield surface even when the peak BHPs were encountered. Such a phenomenon indicates that only tensile failure contributed to the microcracks (according to Eq. (27)), whereas the shear dilation mechanism failed to be activated.

Fig. 13 Mohr’s circles representing the stress states at the boreholes

The shear-induced microcracks are permanent because of elastoplastic failure, while the tensile ones are considered to gradually close with a pressure drawdown(Oldakowski 1994; Samieh 1995; Wong 1999). The field experiences also discovered that if excessive time passed after micro-fracturing before starting steam circulation,the efficiency of reservoir stimulation weakened in certain degrees because the expanded pores and the tensile microcracks gradually closed with time if the borehole pressure was not remained relatively high (usually 2 –2.5 MPa beyond the hydrostatic pressure). It is expected that a sharp increase in injection rates may be able to enforce a BHP high enough to induce shear-type microcracks.

6.2 In fluences of injection and micro‑fracturing on the porothermoelastic responses

As aforementioned, the fluid flow and heat convection are not only strengthened by the dilation of pore space because of increased pore pressure, but also facilitated by the microcracks that provide greater transport paths for both fluid and heat fluxes. A deeper look into these mechanisms can be accomplished by investigating the distributions of pressure,porosity, permeability and temperature along the vertical straight line passing both I and P wells, which are displayed in Fig. 14. It must be noted that because the I and P wells arrived at their peak BHPs at different times (Fig. 6), the pro files of the P well whose peak BHP came first are presented in Fig. 14 (the red curves).

Fig. 14 Distributions of a pressure, b porosity, c permeability and d temperature along the vertical line connecting the I and P wells

Figure 14a shows that as the injection proceeded, the pore pressure across the entire formation depth rises in various degrees, with the effect maximizing at the two boreholes.There are two obvious pressure kinks at the borehole locations until micro-fracturing occurred. Upon micro-fracturing, these two kinks disappear as the corresponding pressures decrease drastically, while the pressure at positions abovez= 19 m keeps on increasing. This phenomenon implies that the micro-fracturing renders the pressure dissipate remarkably faster at the boreholes than other parts of the formation, producing a more evenly distributed pressure pro file. The porosity pro files exhibit a similar trend to that of the pressure pro files (c.f., Fig. 14a and b), except that two slight kinks remain at the boreholes. The porosity neighboring the borehole positions kept on increasing until microfracturing happened, then decreased signi ficantly atzbelow 19 m, above which it continued to increase. Such a trend is in response to the tendency of pressure changes as discussed above. It must be noted that the effects of microcracks were implemented as a dramatic increase in permeability, while their contribution to pore volume was not accounted for in the present simulation. Therefore, the porosity as shown represents that of the formation matrix. Figure 14c shows that the permeability pro files are almost identical with the porosity ones. Three circular dots are marked on the pro files denoting the places where the formation has been broken to form microcracks. At those places (circular shapes with a radius of 0.3 m, Fig. 9b), the permeability increased by a factor of six of its value that had already been raised because of pressure and thermal expansion of pores (Eq. (21)). It is worth notifying that no dot exists at the I well position at the peak BHP curve in Fig. 14c. This is because the peak pressure at the I well that activated micro-fracturing occurred 1.5 h later than at the P well (Fig. 6). It is also interesting to observe that the permeability at the P well decreased slightly (from 3927 to 3794 μD) as the pressure went down and the porosity decreased. Figure 14d demonstrates that the temperatures at the I and P wells were 29.1 and 29.4 °C at the end of the 2nd injection stage, increased further to 30.8and 30.9 °C at the peak BHP (of the P well) and dropped down to 26.2 and 26.4 °C after micro-fracturing. The trend of the three curves implies that thermal convection is most signi ficant at the borehole region especially when microcracks were generated. In contrast, the formation is minimally affected by thermal dissipation from the borehole sources atzgreater than 14 m. The region belowz= 14 m other than the borehole positions experienced only a slight increase in temperature. It is expected that a more considerable amount of hot water being injected at a higher rate would induce propagation of microcracks to a further reach of the formation and bring in a larger degree of heat convection.

7 Conclusions

In this study, a mathematical approach that couples thermal,hydraulic and mechanical responses of a porous medium is proposed. The approach is able to simulate both injection and micro-fracturing processes occurring in an unconsolidated sandstone formation. A field project in a heavy oil sand formation was taken as an example to illustrate its implementation. Several conclusions can be drawn as follows.

The approach deals with the formation permeability as a function of porosity dilation and temperature change, while describing the heat conduction and convection to be transiently varying with the proportion of the pore fluid relative to the solid in a unit volume. Moreover, the approach is capable of describing the microcrack generation and propagation as well as the subsequent formation responses,using both shear and tensile failure criteria for continuum media that have previously been validated from laboratory findings in the literature. The predicted BHPs using the approach match closely the field recorded values, which veri fies its applicability.

The injection stages form an elliptical dilated zone wrapping the dual wells, exhibiting highest pore pressure, porosity, permeability and temperature at the borehole positions.The pressure and heat dissipate gradually into the region farther from the boreholes as the injection time passes, both contributing to an expansion of the pore space. After the injection, a semi-elliptical zone with a horizontal span of 22 m and a vertical extent of 16 m develops around the dual wells. The zone has experienced an increase in pore pressure of more than 1 MPa. A circular region of 2.5 m around the dual wells has undergone a temperature increase up to 9.4 °C. After micro-fracturing, a small circular region (with a radius of 30 cm) filled with tensile microcracks forms around each borehole, markedly enhancing fluid transport across the formation, extending the pore pressure semiellipse to a dimension of 38 m in the horizontal direction and 24 m in the vertical direction. The heat convection, on the other hand, is not as signi ficantly in fluenced by the microcracks as the fluid flow. The relatively homogenous dilated zone after micro-fracturing implies good interwell communication, which is supported by the pressure and temperature responses in the closed well while ramping up the BHP or circulating hot water in the other well. The temperature pro file represents that of the oil sand matrix, showing that the injection of hot water only decreases the viscosity of the bitumen to a limited extent.

The mathematical approach proposed in this study provides a robust means for the prediction of porothermoelastic responses in unconsolidated formations, such that the field engineers can readily optimize designs to maximize the stimulation efficiency while maintaining caprock integrity. The approach can further be extended to analyze the injection practices in tight formations, or evaluate the performances of multiple horizontal wells as well as their interferences, as long as the formation deformation stays within the continuum mechanics regime. Further research efforts are dedicated to developing its three-dimensional counterpart to account for the spatial distribution of petrophysical and mechanical heterogeneity. Handling the challenges brought by the hydraulic, mechanical and thermal anisotropy along the horizontal well direction also demands the future investigation.

AcknowledgementsThe outcome of this work was not possible without the financial support from the National Major Science and Technology Projects of China (Grant No. 2017ZX05009-003) and the National Natural Science Foundation of China (No. 51404281). A large amount of precious field data were generously provided by the Xinjiang Oil field Corporation of CNPC. We also owe our deep gratitude to several engineers of the corporation, Ms. Hongjuan You and Mr. Yong Huang,for their wonderful interpretation of the design issues involved in the field project.

Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source,provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons.org/licen ses/by/4.0/.